Category: Matlab
Category Archives: Matlab
How to fit kinetic model to estimate kinetic parameter from experimental data
I am getting errors in the output. Please help be solve the problem. I have attached the Excel file of the experimental data.
Thank you.
%Fermentation data
Xdata = readmatrix(‘Batch1.xlsx’,’Sheet’,’sheet1′,’Range’,’e2:e13′);
Gdata = readmatrix(‘Batch1.xlsx’,’Sheet’,’sheet1′,’Range’,’b2:b13′);
Bdata = readmatrix(‘Batch1.xlsx’,’Sheet’,’sheet1′,’Range’,’c2:c13′);
Edata = readmatrix(‘Batch1.xlsx’,’Sheet’,’sheet1′,’Range’,’d2:d13′);
% Group all data into yariable y
yinv =[Xdata’; Gdata’; Bdata’; Edata’];
%Data for time
timeex = readmatrix(‘batch1.xlsx’,’Sheet’,’sheet1′,’Range’,’a2:a13′);
%Set ODE solver
% Initial cell biomass concentration (Initial condition)
y0=[1.04 0.78 0.00 0.00];
%Fermentation time
tspan = linspace(0,72);
%Set optimization problem
%Let b = matrix of paramters to be determined
% b= [um ks k1 K1G K1B k2 K2G YXG m k3]
b = optimvar(‘b’,10,"LowerBound",0,"UpperBound",72);
%Set functions for ODE solver for solving ODE
function solferment = toODE(b,tspan,y0)
sol = ode45(@(t,y)batchferment(t,y,b),tspan,y0);
solferment = deval(sol,tspan);
end
%Convert function for ODE solving to optimization expression
%To use RtoODE in an objective function, convert the function to an
%optimization expression by using fcn2optimexpr.
myfcn = fcn2optimexpr(@toODE,b,timeex,y0);
%Express the objective function as the sum of squared differences between
%the ODE solution and the solution with true parameters.
SSE = sum(sum((myfcn-yinv).^2));
%Create an optimization problem with the objective function SSE.
prob = optimproblem ("Description","Fit ODE parameters",’ObjectiveSense’,’min’);
%Objective function (to be minimized)
prob.Objective = SSE;
%Show structure of problem
show(prob)
%Solve Problem
%To find the best-fitting parameters x, give an initial guess
%x0 for the solver and call solve.
% Set initial guesses of parameters
initialGuess.b = [0.18 1.0 0.61 0.18 5.85 3.20 16.25 0.11 3.40 3.02];
%Solve optimization problem
[sol,optval] = solve(prob,initialGuess);
%Extract the results
%Fitted Parameters
bfinal =sol.b;
%Sum of square Error
SSEfinal = optval;
%Plot the simulated data and experimental data
%Call ODE to solve an equation using Final Fitted Parameters (bfinal)
solysim = ode45(@(t,y)batchferment(t,y,bfinal),tspan,y0);
%Evaluate the simulated results at specified tspan
ysim = deval(solysim,tspan);
%Plot graphs
figure;
plot(timeex,Xdata,’*r’,tspan,ysim(1,:),’-b’)
legend(‘exp’,’sim’)
xlabel(‘Time (h)’)
ylabel(‘X (g/L)’)
figure;
plot(timeex,Gdata’,’*r’,tspan,ysim(2,:),’-b’)
legend(‘exp’,’sim’)
xlabel(‘Time (h)’)
ylabel(‘G (g/L)’)
figure;
plot(timeex,Bdata’,’*r’,tspan,ysim(3,:),’-b’)
legend(‘exp’,’sim’)
xlabel(‘Time (h)’)
ylabel(‘B (g/L)’)
figure;
plot(timeex,Edata’,’*r’,tspan,ysim(4,:),’-b’)
legend(‘exp’,’sim’)
xlabel(‘Time (h)’)
ylabel(‘E (g/L)’)
%Equations for Batch
%y(1) = X = Biomass Concentration (g/l)
%y(2) = G = Glucose Concentration (g/l)
%y(3) = B = Cellobiose Concentration (g/L)
%y(4) = E = Ethanol Concentration (g/l)
function dydt = batchferment(t,y,b)
y(1) = X;
y(2) = G;
y(3) = B;
y(4) = E;
%%%Growth equations
%dx/dt = (um*X*G)/(KG+G)
u = (b(1)*y(1)*y(2))/(b(2)+y(2));
%u = (b(1)*X*G)/(b(2)+G);
%parameters
b(1) = um (1/h);
b(2) = ks (g/L);
%%%Cellobiose equations
C = b(11)-(0.9*y(2))-(0.947*y(3))-(0.9*y(4)/0.511)-(1.137*(y(1)-1.04));
r1 = (b(3)*C)/(1+(y(2)/b(4))+(y(3)/b(5)));
r2 = (b(6)*y(3))/(1+(y(2)/b(7)));
b(3) = k1;
b(4) = K1G;
b(5) = K1B;
b(6) = k2;
b(7) = K2G;
%C = Cellulose concentration (g/l)
%C = C0-(0.9*G)-(0.947*B)-(0.9*E/0.511)-(1.137*(X-X0))
%r1 = (k1*C)/(1+(G/K1G)+(B/K1B))
%r2 = (k2*B)/(1+(G/K2G)
%C = C0-(0.9*G)-(0.947*B)-(0.9*E/0.511)-(1.137*(X-X0))
%C = 83.65-(0.9*y(2))-(0.947*y(3))-(0.9*y(4)/0.511)-(1.137*(y(1)-1.04))
%%% Glucose equations
rG = ((1/b(8))*(dx/dt))-(b(9)*y(1));
b(8) = YXG;
b(9) = m;
%rG = ((1/YXG)*(dX/dt))-(m*X)
%%% Ethanol Concentration
Et = k3*(u*X)*(1/YXG);
b(10) = k3;
%Material balance equations
dXdt = u*X;
dGdt = (r2/0.95)-rG;
dBdt = (r1/0.947)-r2;
dEdt = Et;
dydt = [dXdt;dGdt;dBdt;dEdt];
end
This is the output after running:
>> Batch1
Unrecognized function or variable ‘X’.
Error in Batch1>batchferment (line 110)
y(1) = X;
Error in Batch1>@(t,y)batchferment(t,y,b) (line 30)
sol = ode45(@(t,y)batchferment(t,y,b),tspan,y0);
Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 104)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
Error in Batch1>toODE (line 30)
sol = ode45(@(t,y)batchferment(t,y,b),tspan,y0);
Error in optim.problemdef.fcn2optimexpr
Error in optim.problemdef.fcn2optimexpr
Error in fcn2optimexpr
Error in Batch1 (line 38)
myfcn = fcn2optimexpr(@toODE,b,timeex,y0);
Caused by:
Function evaluation failed while attempting to determine output size. The function might contain an error, or might not be well-defined at
the automatically-chosen point. To specify output size without function evaluation, use ‘OutputSize’.
>>I am getting errors in the output. Please help be solve the problem. I have attached the Excel file of the experimental data.
Thank you.
%Fermentation data
Xdata = readmatrix(‘Batch1.xlsx’,’Sheet’,’sheet1′,’Range’,’e2:e13′);
Gdata = readmatrix(‘Batch1.xlsx’,’Sheet’,’sheet1′,’Range’,’b2:b13′);
Bdata = readmatrix(‘Batch1.xlsx’,’Sheet’,’sheet1′,’Range’,’c2:c13′);
Edata = readmatrix(‘Batch1.xlsx’,’Sheet’,’sheet1′,’Range’,’d2:d13′);
% Group all data into yariable y
yinv =[Xdata’; Gdata’; Bdata’; Edata’];
%Data for time
timeex = readmatrix(‘batch1.xlsx’,’Sheet’,’sheet1′,’Range’,’a2:a13′);
%Set ODE solver
% Initial cell biomass concentration (Initial condition)
y0=[1.04 0.78 0.00 0.00];
%Fermentation time
tspan = linspace(0,72);
%Set optimization problem
%Let b = matrix of paramters to be determined
% b= [um ks k1 K1G K1B k2 K2G YXG m k3]
b = optimvar(‘b’,10,"LowerBound",0,"UpperBound",72);
%Set functions for ODE solver for solving ODE
function solferment = toODE(b,tspan,y0)
sol = ode45(@(t,y)batchferment(t,y,b),tspan,y0);
solferment = deval(sol,tspan);
end
%Convert function for ODE solving to optimization expression
%To use RtoODE in an objective function, convert the function to an
%optimization expression by using fcn2optimexpr.
myfcn = fcn2optimexpr(@toODE,b,timeex,y0);
%Express the objective function as the sum of squared differences between
%the ODE solution and the solution with true parameters.
SSE = sum(sum((myfcn-yinv).^2));
%Create an optimization problem with the objective function SSE.
prob = optimproblem ("Description","Fit ODE parameters",’ObjectiveSense’,’min’);
%Objective function (to be minimized)
prob.Objective = SSE;
%Show structure of problem
show(prob)
%Solve Problem
%To find the best-fitting parameters x, give an initial guess
%x0 for the solver and call solve.
% Set initial guesses of parameters
initialGuess.b = [0.18 1.0 0.61 0.18 5.85 3.20 16.25 0.11 3.40 3.02];
%Solve optimization problem
[sol,optval] = solve(prob,initialGuess);
%Extract the results
%Fitted Parameters
bfinal =sol.b;
%Sum of square Error
SSEfinal = optval;
%Plot the simulated data and experimental data
%Call ODE to solve an equation using Final Fitted Parameters (bfinal)
solysim = ode45(@(t,y)batchferment(t,y,bfinal),tspan,y0);
%Evaluate the simulated results at specified tspan
ysim = deval(solysim,tspan);
%Plot graphs
figure;
plot(timeex,Xdata,’*r’,tspan,ysim(1,:),’-b’)
legend(‘exp’,’sim’)
xlabel(‘Time (h)’)
ylabel(‘X (g/L)’)
figure;
plot(timeex,Gdata’,’*r’,tspan,ysim(2,:),’-b’)
legend(‘exp’,’sim’)
xlabel(‘Time (h)’)
ylabel(‘G (g/L)’)
figure;
plot(timeex,Bdata’,’*r’,tspan,ysim(3,:),’-b’)
legend(‘exp’,’sim’)
xlabel(‘Time (h)’)
ylabel(‘B (g/L)’)
figure;
plot(timeex,Edata’,’*r’,tspan,ysim(4,:),’-b’)
legend(‘exp’,’sim’)
xlabel(‘Time (h)’)
ylabel(‘E (g/L)’)
%Equations for Batch
%y(1) = X = Biomass Concentration (g/l)
%y(2) = G = Glucose Concentration (g/l)
%y(3) = B = Cellobiose Concentration (g/L)
%y(4) = E = Ethanol Concentration (g/l)
function dydt = batchferment(t,y,b)
y(1) = X;
y(2) = G;
y(3) = B;
y(4) = E;
%%%Growth equations
%dx/dt = (um*X*G)/(KG+G)
u = (b(1)*y(1)*y(2))/(b(2)+y(2));
%u = (b(1)*X*G)/(b(2)+G);
%parameters
b(1) = um (1/h);
b(2) = ks (g/L);
%%%Cellobiose equations
C = b(11)-(0.9*y(2))-(0.947*y(3))-(0.9*y(4)/0.511)-(1.137*(y(1)-1.04));
r1 = (b(3)*C)/(1+(y(2)/b(4))+(y(3)/b(5)));
r2 = (b(6)*y(3))/(1+(y(2)/b(7)));
b(3) = k1;
b(4) = K1G;
b(5) = K1B;
b(6) = k2;
b(7) = K2G;
%C = Cellulose concentration (g/l)
%C = C0-(0.9*G)-(0.947*B)-(0.9*E/0.511)-(1.137*(X-X0))
%r1 = (k1*C)/(1+(G/K1G)+(B/K1B))
%r2 = (k2*B)/(1+(G/K2G)
%C = C0-(0.9*G)-(0.947*B)-(0.9*E/0.511)-(1.137*(X-X0))
%C = 83.65-(0.9*y(2))-(0.947*y(3))-(0.9*y(4)/0.511)-(1.137*(y(1)-1.04))
%%% Glucose equations
rG = ((1/b(8))*(dx/dt))-(b(9)*y(1));
b(8) = YXG;
b(9) = m;
%rG = ((1/YXG)*(dX/dt))-(m*X)
%%% Ethanol Concentration
Et = k3*(u*X)*(1/YXG);
b(10) = k3;
%Material balance equations
dXdt = u*X;
dGdt = (r2/0.95)-rG;
dBdt = (r1/0.947)-r2;
dEdt = Et;
dydt = [dXdt;dGdt;dBdt;dEdt];
end
This is the output after running:
>> Batch1
Unrecognized function or variable ‘X’.
Error in Batch1>batchferment (line 110)
y(1) = X;
Error in Batch1>@(t,y)batchferment(t,y,b) (line 30)
sol = ode45(@(t,y)batchferment(t,y,b),tspan,y0);
Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 104)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
Error in Batch1>toODE (line 30)
sol = ode45(@(t,y)batchferment(t,y,b),tspan,y0);
Error in optim.problemdef.fcn2optimexpr
Error in optim.problemdef.fcn2optimexpr
Error in fcn2optimexpr
Error in Batch1 (line 38)
myfcn = fcn2optimexpr(@toODE,b,timeex,y0);
Caused by:
Function evaluation failed while attempting to determine output size. The function might contain an error, or might not be well-defined at
the automatically-chosen point. To specify output size without function evaluation, use ‘OutputSize’.
>> I am getting errors in the output. Please help be solve the problem. I have attached the Excel file of the experimental data.
Thank you.
%Fermentation data
Xdata = readmatrix(‘Batch1.xlsx’,’Sheet’,’sheet1′,’Range’,’e2:e13′);
Gdata = readmatrix(‘Batch1.xlsx’,’Sheet’,’sheet1′,’Range’,’b2:b13′);
Bdata = readmatrix(‘Batch1.xlsx’,’Sheet’,’sheet1′,’Range’,’c2:c13′);
Edata = readmatrix(‘Batch1.xlsx’,’Sheet’,’sheet1′,’Range’,’d2:d13′);
% Group all data into yariable y
yinv =[Xdata’; Gdata’; Bdata’; Edata’];
%Data for time
timeex = readmatrix(‘batch1.xlsx’,’Sheet’,’sheet1′,’Range’,’a2:a13′);
%Set ODE solver
% Initial cell biomass concentration (Initial condition)
y0=[1.04 0.78 0.00 0.00];
%Fermentation time
tspan = linspace(0,72);
%Set optimization problem
%Let b = matrix of paramters to be determined
% b= [um ks k1 K1G K1B k2 K2G YXG m k3]
b = optimvar(‘b’,10,"LowerBound",0,"UpperBound",72);
%Set functions for ODE solver for solving ODE
function solferment = toODE(b,tspan,y0)
sol = ode45(@(t,y)batchferment(t,y,b),tspan,y0);
solferment = deval(sol,tspan);
end
%Convert function for ODE solving to optimization expression
%To use RtoODE in an objective function, convert the function to an
%optimization expression by using fcn2optimexpr.
myfcn = fcn2optimexpr(@toODE,b,timeex,y0);
%Express the objective function as the sum of squared differences between
%the ODE solution and the solution with true parameters.
SSE = sum(sum((myfcn-yinv).^2));
%Create an optimization problem with the objective function SSE.
prob = optimproblem ("Description","Fit ODE parameters",’ObjectiveSense’,’min’);
%Objective function (to be minimized)
prob.Objective = SSE;
%Show structure of problem
show(prob)
%Solve Problem
%To find the best-fitting parameters x, give an initial guess
%x0 for the solver and call solve.
% Set initial guesses of parameters
initialGuess.b = [0.18 1.0 0.61 0.18 5.85 3.20 16.25 0.11 3.40 3.02];
%Solve optimization problem
[sol,optval] = solve(prob,initialGuess);
%Extract the results
%Fitted Parameters
bfinal =sol.b;
%Sum of square Error
SSEfinal = optval;
%Plot the simulated data and experimental data
%Call ODE to solve an equation using Final Fitted Parameters (bfinal)
solysim = ode45(@(t,y)batchferment(t,y,bfinal),tspan,y0);
%Evaluate the simulated results at specified tspan
ysim = deval(solysim,tspan);
%Plot graphs
figure;
plot(timeex,Xdata,’*r’,tspan,ysim(1,:),’-b’)
legend(‘exp’,’sim’)
xlabel(‘Time (h)’)
ylabel(‘X (g/L)’)
figure;
plot(timeex,Gdata’,’*r’,tspan,ysim(2,:),’-b’)
legend(‘exp’,’sim’)
xlabel(‘Time (h)’)
ylabel(‘G (g/L)’)
figure;
plot(timeex,Bdata’,’*r’,tspan,ysim(3,:),’-b’)
legend(‘exp’,’sim’)
xlabel(‘Time (h)’)
ylabel(‘B (g/L)’)
figure;
plot(timeex,Edata’,’*r’,tspan,ysim(4,:),’-b’)
legend(‘exp’,’sim’)
xlabel(‘Time (h)’)
ylabel(‘E (g/L)’)
%Equations for Batch
%y(1) = X = Biomass Concentration (g/l)
%y(2) = G = Glucose Concentration (g/l)
%y(3) = B = Cellobiose Concentration (g/L)
%y(4) = E = Ethanol Concentration (g/l)
function dydt = batchferment(t,y,b)
y(1) = X;
y(2) = G;
y(3) = B;
y(4) = E;
%%%Growth equations
%dx/dt = (um*X*G)/(KG+G)
u = (b(1)*y(1)*y(2))/(b(2)+y(2));
%u = (b(1)*X*G)/(b(2)+G);
%parameters
b(1) = um (1/h);
b(2) = ks (g/L);
%%%Cellobiose equations
C = b(11)-(0.9*y(2))-(0.947*y(3))-(0.9*y(4)/0.511)-(1.137*(y(1)-1.04));
r1 = (b(3)*C)/(1+(y(2)/b(4))+(y(3)/b(5)));
r2 = (b(6)*y(3))/(1+(y(2)/b(7)));
b(3) = k1;
b(4) = K1G;
b(5) = K1B;
b(6) = k2;
b(7) = K2G;
%C = Cellulose concentration (g/l)
%C = C0-(0.9*G)-(0.947*B)-(0.9*E/0.511)-(1.137*(X-X0))
%r1 = (k1*C)/(1+(G/K1G)+(B/K1B))
%r2 = (k2*B)/(1+(G/K2G)
%C = C0-(0.9*G)-(0.947*B)-(0.9*E/0.511)-(1.137*(X-X0))
%C = 83.65-(0.9*y(2))-(0.947*y(3))-(0.9*y(4)/0.511)-(1.137*(y(1)-1.04))
%%% Glucose equations
rG = ((1/b(8))*(dx/dt))-(b(9)*y(1));
b(8) = YXG;
b(9) = m;
%rG = ((1/YXG)*(dX/dt))-(m*X)
%%% Ethanol Concentration
Et = k3*(u*X)*(1/YXG);
b(10) = k3;
%Material balance equations
dXdt = u*X;
dGdt = (r2/0.95)-rG;
dBdt = (r1/0.947)-r2;
dEdt = Et;
dydt = [dXdt;dGdt;dBdt;dEdt];
end
This is the output after running:
>> Batch1
Unrecognized function or variable ‘X’.
Error in Batch1>batchferment (line 110)
y(1) = X;
Error in Batch1>@(t,y)batchferment(t,y,b) (line 30)
sol = ode45(@(t,y)batchferment(t,y,b),tspan,y0);
Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 104)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
Error in Batch1>toODE (line 30)
sol = ode45(@(t,y)batchferment(t,y,b),tspan,y0);
Error in optim.problemdef.fcn2optimexpr
Error in optim.problemdef.fcn2optimexpr
Error in fcn2optimexpr
Error in Batch1 (line 38)
myfcn = fcn2optimexpr(@toODE,b,timeex,y0);
Caused by:
Function evaluation failed while attempting to determine output size. The function might contain an error, or might not be well-defined at
the automatically-chosen point. To specify output size without function evaluation, use ‘OutputSize’.
>> data fitting, parameter estimation, optimization, isqnonlin method MATLAB Answers — New Questions
How to set the font of the Diagnostic Viewer in Simulink?
Hi All,
I met some problem in the Diagnostic Viewer in Simulink, the font of the text has became extremely which I can hardly see. Is there anybody know how to fix this problem?
Regards,
KeranHi All,
I met some problem in the Diagnostic Viewer in Simulink, the font of the text has became extremely which I can hardly see. Is there anybody know how to fix this problem?
Regards,
Keran Hi All,
I met some problem in the Diagnostic Viewer in Simulink, the font of the text has became extremely which I can hardly see. Is there anybody know how to fix this problem?
Regards,
Keran simulink, font MATLAB Answers — New Questions
Where to find which toolbox a MATLAB example belongs to, or which toolboxes I need to run the example?
I am currently working through a MathWorks example using the "fit" function:
https://www.mathworks.com/help/curvefit/fit-exponential-model-to-census-data.html
But the example documentation does not specify which toolbox this "fit" function belongs to. Where can I find which toolboxes I need for running a MATLAB example?I am currently working through a MathWorks example using the "fit" function:
https://www.mathworks.com/help/curvefit/fit-exponential-model-to-census-data.html
But the example documentation does not specify which toolbox this "fit" function belongs to. Where can I find which toolboxes I need for running a MATLAB example? I am currently working through a MathWorks example using the "fit" function:
https://www.mathworks.com/help/curvefit/fit-exponential-model-to-census-data.html
But the example documentation does not specify which toolbox this "fit" function belongs to. Where can I find which toolboxes I need for running a MATLAB example? toolboxes, documentation MATLAB Answers — New Questions
How to implement custom function in Flight Log Analyzer
I am trying to create a custom function within the Flight Log Analyzer app in the UAV Toolbox, and I am getting errors no matter what I try. I am simply trying to plot a Fast-Fourier Transform of Accel and Gyro data. Below is an example of a function I wrote that works on its own, but I can not get it to work within the app. I am wondering what the function sturcture is that the app is looking for. I see what is on the pop up window when you go to create the function, but even when I try creating more arguments for the function I get index mismatch errors. An example of a function that can be passed into the app would be very helpful.
(in this example, "signalData" would be ACC_0 from the workspace file)
function afft = ArdupilotFFT(singalData)
Fs = 2000;
L = length(singalData);
L2 = floor(L/2);
f0 = Fs/L*(0:L2-1);
X0 = fft(singalData(:,5));
X0 = 10*log10(abs(X0(1:L2)));
Y0 = fft(singalData(:,6));
Y0 = 10*log10(abs(Y0(1:L2)));
Z0 = fft(singalData(:,7));
Z0 = 10*log10(abs(Z0(1:L2)));
figure(‘Name’,’FFT’);
%subplot(1,3,1);
semilogx(f0,X0,f0,Y0,f0,Z0)
grid on; xlabel(‘Frequency (Hz)’); xlim([1,1000]);
endI am trying to create a custom function within the Flight Log Analyzer app in the UAV Toolbox, and I am getting errors no matter what I try. I am simply trying to plot a Fast-Fourier Transform of Accel and Gyro data. Below is an example of a function I wrote that works on its own, but I can not get it to work within the app. I am wondering what the function sturcture is that the app is looking for. I see what is on the pop up window when you go to create the function, but even when I try creating more arguments for the function I get index mismatch errors. An example of a function that can be passed into the app would be very helpful.
(in this example, "signalData" would be ACC_0 from the workspace file)
function afft = ArdupilotFFT(singalData)
Fs = 2000;
L = length(singalData);
L2 = floor(L/2);
f0 = Fs/L*(0:L2-1);
X0 = fft(singalData(:,5));
X0 = 10*log10(abs(X0(1:L2)));
Y0 = fft(singalData(:,6));
Y0 = 10*log10(abs(Y0(1:L2)));
Z0 = fft(singalData(:,7));
Z0 = 10*log10(abs(Z0(1:L2)));
figure(‘Name’,’FFT’);
%subplot(1,3,1);
semilogx(f0,X0,f0,Y0,f0,Z0)
grid on; xlabel(‘Frequency (Hz)’); xlim([1,1000]);
end I am trying to create a custom function within the Flight Log Analyzer app in the UAV Toolbox, and I am getting errors no matter what I try. I am simply trying to plot a Fast-Fourier Transform of Accel and Gyro data. Below is an example of a function I wrote that works on its own, but I can not get it to work within the app. I am wondering what the function sturcture is that the app is looking for. I see what is on the pop up window when you go to create the function, but even when I try creating more arguments for the function I get index mismatch errors. An example of a function that can be passed into the app would be very helpful.
(in this example, "signalData" would be ACC_0 from the workspace file)
function afft = ArdupilotFFT(singalData)
Fs = 2000;
L = length(singalData);
L2 = floor(L/2);
f0 = Fs/L*(0:L2-1);
X0 = fft(singalData(:,5));
X0 = 10*log10(abs(X0(1:L2)));
Y0 = fft(singalData(:,6));
Y0 = 10*log10(abs(Y0(1:L2)));
Z0 = fft(singalData(:,7));
Z0 = 10*log10(abs(Z0(1:L2)));
figure(‘Name’,’FFT’);
%subplot(1,3,1);
semilogx(f0,X0,f0,Y0,f0,Z0)
grid on; xlabel(‘Frequency (Hz)’); xlim([1,1000]);
end uav, flightloganalyzer, fft, function MATLAB Answers — New Questions
How do I determine if a function is user-defined or belongs to a MATLAB toolbox?
How do I determine if a function is user-defined or belongs to a MATLAB toolbox?How do I determine if a function is user-defined or belongs to a MATLAB toolbox? How do I determine if a function is user-defined or belongs to a MATLAB toolbox? locate, function, name, source MATLAB Answers — New Questions
Why does my Simulink Coder build, Rapid Accelerator build or FMU Export fail with “fatal error C1083: Cannot open include file: ‘xxx.h’: No such file or directory” when using Visual Studio C++ compiler, referencing C Standard Library headers?
I encountered errors while running a Simulink Coder build, Rapid Accelerator build, or FMU Export using the Visual Studio C++ compiler on Windows. The errors were related to missing C Standard Library headers like "limits.h", "string.h", "stdlib.h", or "stddef.h":
fatal error C1083: Cannot open include file: ‘xxx.h’: No such file or directory
The same build works with the MinGW64 compiler. I verified that my Visual Studio setup is correct by successfully compiling a "Hello World" C program and creating a MEX file from MATLAB.I encountered errors while running a Simulink Coder build, Rapid Accelerator build, or FMU Export using the Visual Studio C++ compiler on Windows. The errors were related to missing C Standard Library headers like "limits.h", "string.h", "stdlib.h", or "stddef.h":
fatal error C1083: Cannot open include file: ‘xxx.h’: No such file or directory
The same build works with the MinGW64 compiler. I verified that my Visual Studio setup is correct by successfully compiling a "Hello World" C program and creating a MEX file from MATLAB. I encountered errors while running a Simulink Coder build, Rapid Accelerator build, or FMU Export using the Visual Studio C++ compiler on Windows. The errors were related to missing C Standard Library headers like "limits.h", "string.h", "stdlib.h", or "stddef.h":
fatal error C1083: Cannot open include file: ‘xxx.h’: No such file or directory
The same build works with the MinGW64 compiler. I verified that my Visual Studio setup is correct by successfully compiling a "Hello World" C program and creating a MEX file from MATLAB. raccel, rapid, accel, stdlib, cstdlib, cstd MATLAB Answers — New Questions
How to input a novel boundary condition for a coupled PDE system
Hi all,
I am trying to simulate a coupled PDE system with a non typical boundary conditions.
My coupled PDE system is as such.
The system is solved for two time periods
For first time period 0<tau_M1<tau
For s at the L.H.S we have
For s at the R.H.S we have
For species m at the L.H.S we have the to set m value to a constant
For species m at the R.H.S we have the no flux condition
For second time period tau<tau_M1<tau_2
For s at the L.H.S and R.H.S we have the same such as
For m we have
For the species we have the L.H.S boundary condition for timr period (tau<tau_M1<tau_2). It is this boundary condtiuon that I am seeking help with. How can I go about this.
function [c25, c2, vode] = pde1dm_PS_OCP_v1()
% Piecewise PDE Implementation
% Pulsing Simulation
% Constants Used in the Nernst Equation
n = 2; % No. of Electrons Involved
F = 96485.3329; % C/mol
R = 8.314; % J/(mol.K) or CV/(mol.K)
Temp = 273.15+20; % K
N = 21;
m = 0;
kappa=1;
eta=1;
gamma=1;
mu=1;
tau_M1=1;
tau_M2=2;
l=1;
D_r=1;
Area = 0.0011; % Geometric Area (cm^2)
m_layer = 4.13E-05; % m_layer (mol/cm^3)
D_m = 8.14E-08; % cm^2/s (De of polymer)
k = n*F*Area*l*m_layer;
% DIMENSIONLESS PARAMETERS
chi = linspace(0, 1, N);
tau = linspace(0, tau_M1, N); % Dimensionless Time
% FIRST POTENTIAL STEP
E0 = 0.22;
E_set = .450;
epsilon1 = ((n*F)/(R*Temp))*(E_set-E0); % Dimensionless Potential
c_M_ox = 1/(1+exp(-epsilon1)); % Mox BC
ic_arg_1 = {@(x)ones(size(N)) ; @(x)zeros(size(N))};
IC = @(x)PDE_ND_PS_EK_IC(x, ic_arg_1);
BC = @(xl, yl, xr, yr, t)PDE_PS_EK_BC(xl, yl, xr, yr, c_M_ox);
sol1 = pde1dm(m, @(x, t, c, DcDx)PDE_ND_PS_EK(chi, tau, c, DcDx, kappa, eta, gamma, mu, D_r), …
IC, BC, chi, tau);
c1 = sol1(:, :, 1); % Substrate Conc.
c2 = sol1(:, :, 2); % Mox Conc.
% OPEN CIRCUIT POTENTIAL
tau2 = linspace(tau_M1, tau_M2, N); % Dimensionless Time
ic_arg_1 = {@(x)interp1(chi, sol1(N, :, 1), x) ; @(x)interp1(chi, sol1(N, :, 2), x)};
IC2 = @(x)PDE_ND_PS_EK_IC(x, ic_arg_1);
BC2 = @(xl, yl, xr, yr, t, v)PDE_PSw_EK_BC_2(xl, yl, xr, yr, v);
ode_IC = @() ode_ICFun(tau_M1);
opts.vectorized=’on’;
% sol2 = pde1dm(m, @(x, t, c, DcDx)PDE_ND_PS_EK(x, t, c, DcDx, kappa, eta, gamma, mu, D_r), …
% IC2, BC2, chi, tau2);
[sol2, vode] = pde1dm(m, @(x, t, c, DcDx)PDE_ND_PS_EK(x, t, c, DcDx, kappa, eta, gamma, mu, D_r), …
IC2, BC2, chi, tau2, @odeFunc, ode_IC, .99);
% Concentration Profiles c(i, j, k)(Solutions)
c14 = [c1; sol2(:, :, 1)]; % Substrate Conc.
c25 = [c2; sol2(:, :, 2)]; % Mox Conc.
c36 = 1-c25; % Mred Conc.
mox = abs(sol2(:, 1, 2));
mred = 1-mox;
E_PS = E_set.*ones(1,N);
for counter2 = 1:N
E_OCP(counter2) = E0 + (((R*Temp)/(n*F).*(log(mox(counter2,:)./mred(counter2,:)))));
end
E_array = [E_PS, E_OCP];
% This is the case finder function
% Video(kappa, eta, gamma, mu, tau_M1, tau_M2, l, D_m, N, chi, c14, c25, c36, E_array)
function [cc, ff, ss] = PDE_ND_PS_EK(chi, tau, c, DcDx, kappa, eta, gamma, mu, D_r)
% S; Mox;
cc = [ D_r; 1];
ff = [ 1; 1].*DcDx;
S_kin = ((gamma/eta)*(kappa^2)*c(1)*c(2))./…
(gamma.*c(2).*(1+(mu.*c(1)))+c(1));
M_kin = ((kappa^2).*c(1).*c(2))./…
(gamma.*c(2).*(1+(mu.*c(1)))+c(1));
ss = [-S_kin; -M_kin];
end
function c0 = PDE_ND_PS_EK_IC(x, ic_arg_1)
% Substrate; Mox;
c0 = [ic_arg_1{1}(x).’; ic_arg_1{2}(x).’];
end
function [pl, ql, pr, qr] = PDE_PS_EK_BC(xl, yl, xr, yr, c_M_ox)
% ———————
% | | |
% | | |
% | | |
% | | |
% | | |
% | | |
% | | |
% ——-|————-
% pl pr
% Substrate; Mediator;
pl = [0 ; yl(2)-c_M_ox];
ql = [1 ; 0];
pr = [yr(1)-1 ; 0];
qr = [0 ; 1];
end
function [pl, ql, pr, qr] = PDE_PSw_EK_BC_2(xl, yl, xr, yr, t, v, vdot)
pl = [0 ; 0];
ql = [1 ; 1];
pr = [yr(1)-1 ; 0];
qr = [0 ; 1];
end
function f = odeFunc(x, t, c, DcDx, v, vdot)
f = ((1E-03./(v(1)*(1-v(1))))*(DcDx))-vdot(1);
end
function v0=ode_ICFun(tau_0)
v0 = ones(tau_0).*.9;
end
function myVideo = Video(kappa, eta, gamma, mu, tau_M1, tau_M2, l, D_m, N, chi, c14, c25, c36, E_array)
% Initialize Video
G = figure(1);
myVideo = VideoWriter(sprintf(‘κ%.2f γ%.2f η%.2f µ%.2f’, kappa, gamma, eta, mu));
myVideo.FrameRate = 10; % can adjust this, 5 – 10 works well for me
myVideo.Quality = 100;
open(myVideo);
color = ‘red’;
u = uicontrol(G, ‘Style’,’slider’,’Position’,[0 40 10 360],…
‘Min’,0,’Max’,N*2,’Value’,0);
tspan1 = linspace(0, round(((tau_M1*(l^2))/D_m)), N);
tspan2 = linspace(round(((tau_M1*(l^2))/D_m)), round(((tau_M2*l^2)/D_m),2), N);
tspan = [tspan1,tspan2];
for ii = 2:(N*2)
% Plot Species Conc. in the Layer
subplot(2,1,1);
yyaxis left
plot(chi,c25(ii,:),chi,c36(ii,:), ‘LineWidth’, 2);
ylabel(‘M’)
ylim([0,1]);
yyaxis right
plot(chi, c14(ii,:), ‘LineWidth’, 2);
ylabel(‘S’)
ylim([0,1]);
xlabel(‘chi’);
legend(‘M_{ox}’,’M_{red}’, ‘S’);
title(‘Normalized Concentration’);
subplot(2,1,2);
addpoints(animatedline(tspan,E_array,’marker’, ‘.’, ‘markersize’, 6, ‘color’, color,…
‘linestyle’, ‘–‘, ‘MaximumNumPoints’, 1),tspan(ii),E_array(ii));
ylim([0,0.5]);
xlim([0, tspan(end)]);
hold on;
title([‘Potential = ‘, num2str(E_array(ii))]);
% text(8, 8,{”,”})
drawnow
xlabel(‘t / s’); ylabel(‘E / V’);
u.Value = ii;
uicontrol(‘Style’,’Edit’,’Position’,[0,00,40,40], …
‘String’,num2str(tspan(ii),3));
pause(0.001)
M = getframe(G);
writeVideo(myVideo, M);
end
end
endHi all,
I am trying to simulate a coupled PDE system with a non typical boundary conditions.
My coupled PDE system is as such.
The system is solved for two time periods
For first time period 0<tau_M1<tau
For s at the L.H.S we have
For s at the R.H.S we have
For species m at the L.H.S we have the to set m value to a constant
For species m at the R.H.S we have the no flux condition
For second time period tau<tau_M1<tau_2
For s at the L.H.S and R.H.S we have the same such as
For m we have
For the species we have the L.H.S boundary condition for timr period (tau<tau_M1<tau_2). It is this boundary condtiuon that I am seeking help with. How can I go about this.
function [c25, c2, vode] = pde1dm_PS_OCP_v1()
% Piecewise PDE Implementation
% Pulsing Simulation
% Constants Used in the Nernst Equation
n = 2; % No. of Electrons Involved
F = 96485.3329; % C/mol
R = 8.314; % J/(mol.K) or CV/(mol.K)
Temp = 273.15+20; % K
N = 21;
m = 0;
kappa=1;
eta=1;
gamma=1;
mu=1;
tau_M1=1;
tau_M2=2;
l=1;
D_r=1;
Area = 0.0011; % Geometric Area (cm^2)
m_layer = 4.13E-05; % m_layer (mol/cm^3)
D_m = 8.14E-08; % cm^2/s (De of polymer)
k = n*F*Area*l*m_layer;
% DIMENSIONLESS PARAMETERS
chi = linspace(0, 1, N);
tau = linspace(0, tau_M1, N); % Dimensionless Time
% FIRST POTENTIAL STEP
E0 = 0.22;
E_set = .450;
epsilon1 = ((n*F)/(R*Temp))*(E_set-E0); % Dimensionless Potential
c_M_ox = 1/(1+exp(-epsilon1)); % Mox BC
ic_arg_1 = {@(x)ones(size(N)) ; @(x)zeros(size(N))};
IC = @(x)PDE_ND_PS_EK_IC(x, ic_arg_1);
BC = @(xl, yl, xr, yr, t)PDE_PS_EK_BC(xl, yl, xr, yr, c_M_ox);
sol1 = pde1dm(m, @(x, t, c, DcDx)PDE_ND_PS_EK(chi, tau, c, DcDx, kappa, eta, gamma, mu, D_r), …
IC, BC, chi, tau);
c1 = sol1(:, :, 1); % Substrate Conc.
c2 = sol1(:, :, 2); % Mox Conc.
% OPEN CIRCUIT POTENTIAL
tau2 = linspace(tau_M1, tau_M2, N); % Dimensionless Time
ic_arg_1 = {@(x)interp1(chi, sol1(N, :, 1), x) ; @(x)interp1(chi, sol1(N, :, 2), x)};
IC2 = @(x)PDE_ND_PS_EK_IC(x, ic_arg_1);
BC2 = @(xl, yl, xr, yr, t, v)PDE_PSw_EK_BC_2(xl, yl, xr, yr, v);
ode_IC = @() ode_ICFun(tau_M1);
opts.vectorized=’on’;
% sol2 = pde1dm(m, @(x, t, c, DcDx)PDE_ND_PS_EK(x, t, c, DcDx, kappa, eta, gamma, mu, D_r), …
% IC2, BC2, chi, tau2);
[sol2, vode] = pde1dm(m, @(x, t, c, DcDx)PDE_ND_PS_EK(x, t, c, DcDx, kappa, eta, gamma, mu, D_r), …
IC2, BC2, chi, tau2, @odeFunc, ode_IC, .99);
% Concentration Profiles c(i, j, k)(Solutions)
c14 = [c1; sol2(:, :, 1)]; % Substrate Conc.
c25 = [c2; sol2(:, :, 2)]; % Mox Conc.
c36 = 1-c25; % Mred Conc.
mox = abs(sol2(:, 1, 2));
mred = 1-mox;
E_PS = E_set.*ones(1,N);
for counter2 = 1:N
E_OCP(counter2) = E0 + (((R*Temp)/(n*F).*(log(mox(counter2,:)./mred(counter2,:)))));
end
E_array = [E_PS, E_OCP];
% This is the case finder function
% Video(kappa, eta, gamma, mu, tau_M1, tau_M2, l, D_m, N, chi, c14, c25, c36, E_array)
function [cc, ff, ss] = PDE_ND_PS_EK(chi, tau, c, DcDx, kappa, eta, gamma, mu, D_r)
% S; Mox;
cc = [ D_r; 1];
ff = [ 1; 1].*DcDx;
S_kin = ((gamma/eta)*(kappa^2)*c(1)*c(2))./…
(gamma.*c(2).*(1+(mu.*c(1)))+c(1));
M_kin = ((kappa^2).*c(1).*c(2))./…
(gamma.*c(2).*(1+(mu.*c(1)))+c(1));
ss = [-S_kin; -M_kin];
end
function c0 = PDE_ND_PS_EK_IC(x, ic_arg_1)
% Substrate; Mox;
c0 = [ic_arg_1{1}(x).’; ic_arg_1{2}(x).’];
end
function [pl, ql, pr, qr] = PDE_PS_EK_BC(xl, yl, xr, yr, c_M_ox)
% ———————
% | | |
% | | |
% | | |
% | | |
% | | |
% | | |
% | | |
% ——-|————-
% pl pr
% Substrate; Mediator;
pl = [0 ; yl(2)-c_M_ox];
ql = [1 ; 0];
pr = [yr(1)-1 ; 0];
qr = [0 ; 1];
end
function [pl, ql, pr, qr] = PDE_PSw_EK_BC_2(xl, yl, xr, yr, t, v, vdot)
pl = [0 ; 0];
ql = [1 ; 1];
pr = [yr(1)-1 ; 0];
qr = [0 ; 1];
end
function f = odeFunc(x, t, c, DcDx, v, vdot)
f = ((1E-03./(v(1)*(1-v(1))))*(DcDx))-vdot(1);
end
function v0=ode_ICFun(tau_0)
v0 = ones(tau_0).*.9;
end
function myVideo = Video(kappa, eta, gamma, mu, tau_M1, tau_M2, l, D_m, N, chi, c14, c25, c36, E_array)
% Initialize Video
G = figure(1);
myVideo = VideoWriter(sprintf(‘κ%.2f γ%.2f η%.2f µ%.2f’, kappa, gamma, eta, mu));
myVideo.FrameRate = 10; % can adjust this, 5 – 10 works well for me
myVideo.Quality = 100;
open(myVideo);
color = ‘red’;
u = uicontrol(G, ‘Style’,’slider’,’Position’,[0 40 10 360],…
‘Min’,0,’Max’,N*2,’Value’,0);
tspan1 = linspace(0, round(((tau_M1*(l^2))/D_m)), N);
tspan2 = linspace(round(((tau_M1*(l^2))/D_m)), round(((tau_M2*l^2)/D_m),2), N);
tspan = [tspan1,tspan2];
for ii = 2:(N*2)
% Plot Species Conc. in the Layer
subplot(2,1,1);
yyaxis left
plot(chi,c25(ii,:),chi,c36(ii,:), ‘LineWidth’, 2);
ylabel(‘M’)
ylim([0,1]);
yyaxis right
plot(chi, c14(ii,:), ‘LineWidth’, 2);
ylabel(‘S’)
ylim([0,1]);
xlabel(‘chi’);
legend(‘M_{ox}’,’M_{red}’, ‘S’);
title(‘Normalized Concentration’);
subplot(2,1,2);
addpoints(animatedline(tspan,E_array,’marker’, ‘.’, ‘markersize’, 6, ‘color’, color,…
‘linestyle’, ‘–‘, ‘MaximumNumPoints’, 1),tspan(ii),E_array(ii));
ylim([0,0.5]);
xlim([0, tspan(end)]);
hold on;
title([‘Potential = ‘, num2str(E_array(ii))]);
% text(8, 8,{”,”})
drawnow
xlabel(‘t / s’); ylabel(‘E / V’);
u.Value = ii;
uicontrol(‘Style’,’Edit’,’Position’,[0,00,40,40], …
‘String’,num2str(tspan(ii),3));
pause(0.001)
M = getframe(G);
writeVideo(myVideo, M);
end
end
end Hi all,
I am trying to simulate a coupled PDE system with a non typical boundary conditions.
My coupled PDE system is as such.
The system is solved for two time periods
For first time period 0<tau_M1<tau
For s at the L.H.S we have
For s at the R.H.S we have
For species m at the L.H.S we have the to set m value to a constant
For species m at the R.H.S we have the no flux condition
For second time period tau<tau_M1<tau_2
For s at the L.H.S and R.H.S we have the same such as
For m we have
For the species we have the L.H.S boundary condition for timr period (tau<tau_M1<tau_2). It is this boundary condtiuon that I am seeking help with. How can I go about this.
function [c25, c2, vode] = pde1dm_PS_OCP_v1()
% Piecewise PDE Implementation
% Pulsing Simulation
% Constants Used in the Nernst Equation
n = 2; % No. of Electrons Involved
F = 96485.3329; % C/mol
R = 8.314; % J/(mol.K) or CV/(mol.K)
Temp = 273.15+20; % K
N = 21;
m = 0;
kappa=1;
eta=1;
gamma=1;
mu=1;
tau_M1=1;
tau_M2=2;
l=1;
D_r=1;
Area = 0.0011; % Geometric Area (cm^2)
m_layer = 4.13E-05; % m_layer (mol/cm^3)
D_m = 8.14E-08; % cm^2/s (De of polymer)
k = n*F*Area*l*m_layer;
% DIMENSIONLESS PARAMETERS
chi = linspace(0, 1, N);
tau = linspace(0, tau_M1, N); % Dimensionless Time
% FIRST POTENTIAL STEP
E0 = 0.22;
E_set = .450;
epsilon1 = ((n*F)/(R*Temp))*(E_set-E0); % Dimensionless Potential
c_M_ox = 1/(1+exp(-epsilon1)); % Mox BC
ic_arg_1 = {@(x)ones(size(N)) ; @(x)zeros(size(N))};
IC = @(x)PDE_ND_PS_EK_IC(x, ic_arg_1);
BC = @(xl, yl, xr, yr, t)PDE_PS_EK_BC(xl, yl, xr, yr, c_M_ox);
sol1 = pde1dm(m, @(x, t, c, DcDx)PDE_ND_PS_EK(chi, tau, c, DcDx, kappa, eta, gamma, mu, D_r), …
IC, BC, chi, tau);
c1 = sol1(:, :, 1); % Substrate Conc.
c2 = sol1(:, :, 2); % Mox Conc.
% OPEN CIRCUIT POTENTIAL
tau2 = linspace(tau_M1, tau_M2, N); % Dimensionless Time
ic_arg_1 = {@(x)interp1(chi, sol1(N, :, 1), x) ; @(x)interp1(chi, sol1(N, :, 2), x)};
IC2 = @(x)PDE_ND_PS_EK_IC(x, ic_arg_1);
BC2 = @(xl, yl, xr, yr, t, v)PDE_PSw_EK_BC_2(xl, yl, xr, yr, v);
ode_IC = @() ode_ICFun(tau_M1);
opts.vectorized=’on’;
% sol2 = pde1dm(m, @(x, t, c, DcDx)PDE_ND_PS_EK(x, t, c, DcDx, kappa, eta, gamma, mu, D_r), …
% IC2, BC2, chi, tau2);
[sol2, vode] = pde1dm(m, @(x, t, c, DcDx)PDE_ND_PS_EK(x, t, c, DcDx, kappa, eta, gamma, mu, D_r), …
IC2, BC2, chi, tau2, @odeFunc, ode_IC, .99);
% Concentration Profiles c(i, j, k)(Solutions)
c14 = [c1; sol2(:, :, 1)]; % Substrate Conc.
c25 = [c2; sol2(:, :, 2)]; % Mox Conc.
c36 = 1-c25; % Mred Conc.
mox = abs(sol2(:, 1, 2));
mred = 1-mox;
E_PS = E_set.*ones(1,N);
for counter2 = 1:N
E_OCP(counter2) = E0 + (((R*Temp)/(n*F).*(log(mox(counter2,:)./mred(counter2,:)))));
end
E_array = [E_PS, E_OCP];
% This is the case finder function
% Video(kappa, eta, gamma, mu, tau_M1, tau_M2, l, D_m, N, chi, c14, c25, c36, E_array)
function [cc, ff, ss] = PDE_ND_PS_EK(chi, tau, c, DcDx, kappa, eta, gamma, mu, D_r)
% S; Mox;
cc = [ D_r; 1];
ff = [ 1; 1].*DcDx;
S_kin = ((gamma/eta)*(kappa^2)*c(1)*c(2))./…
(gamma.*c(2).*(1+(mu.*c(1)))+c(1));
M_kin = ((kappa^2).*c(1).*c(2))./…
(gamma.*c(2).*(1+(mu.*c(1)))+c(1));
ss = [-S_kin; -M_kin];
end
function c0 = PDE_ND_PS_EK_IC(x, ic_arg_1)
% Substrate; Mox;
c0 = [ic_arg_1{1}(x).’; ic_arg_1{2}(x).’];
end
function [pl, ql, pr, qr] = PDE_PS_EK_BC(xl, yl, xr, yr, c_M_ox)
% ———————
% | | |
% | | |
% | | |
% | | |
% | | |
% | | |
% | | |
% ——-|————-
% pl pr
% Substrate; Mediator;
pl = [0 ; yl(2)-c_M_ox];
ql = [1 ; 0];
pr = [yr(1)-1 ; 0];
qr = [0 ; 1];
end
function [pl, ql, pr, qr] = PDE_PSw_EK_BC_2(xl, yl, xr, yr, t, v, vdot)
pl = [0 ; 0];
ql = [1 ; 1];
pr = [yr(1)-1 ; 0];
qr = [0 ; 1];
end
function f = odeFunc(x, t, c, DcDx, v, vdot)
f = ((1E-03./(v(1)*(1-v(1))))*(DcDx))-vdot(1);
end
function v0=ode_ICFun(tau_0)
v0 = ones(tau_0).*.9;
end
function myVideo = Video(kappa, eta, gamma, mu, tau_M1, tau_M2, l, D_m, N, chi, c14, c25, c36, E_array)
% Initialize Video
G = figure(1);
myVideo = VideoWriter(sprintf(‘κ%.2f γ%.2f η%.2f µ%.2f’, kappa, gamma, eta, mu));
myVideo.FrameRate = 10; % can adjust this, 5 – 10 works well for me
myVideo.Quality = 100;
open(myVideo);
color = ‘red’;
u = uicontrol(G, ‘Style’,’slider’,’Position’,[0 40 10 360],…
‘Min’,0,’Max’,N*2,’Value’,0);
tspan1 = linspace(0, round(((tau_M1*(l^2))/D_m)), N);
tspan2 = linspace(round(((tau_M1*(l^2))/D_m)), round(((tau_M2*l^2)/D_m),2), N);
tspan = [tspan1,tspan2];
for ii = 2:(N*2)
% Plot Species Conc. in the Layer
subplot(2,1,1);
yyaxis left
plot(chi,c25(ii,:),chi,c36(ii,:), ‘LineWidth’, 2);
ylabel(‘M’)
ylim([0,1]);
yyaxis right
plot(chi, c14(ii,:), ‘LineWidth’, 2);
ylabel(‘S’)
ylim([0,1]);
xlabel(‘chi’);
legend(‘M_{ox}’,’M_{red}’, ‘S’);
title(‘Normalized Concentration’);
subplot(2,1,2);
addpoints(animatedline(tspan,E_array,’marker’, ‘.’, ‘markersize’, 6, ‘color’, color,…
‘linestyle’, ‘–‘, ‘MaximumNumPoints’, 1),tspan(ii),E_array(ii));
ylim([0,0.5]);
xlim([0, tspan(end)]);
hold on;
title([‘Potential = ‘, num2str(E_array(ii))]);
% text(8, 8,{”,”})
drawnow
xlabel(‘t / s’); ylabel(‘E / V’);
u.Value = ii;
uicontrol(‘Style’,’Edit’,’Position’,[0,00,40,40], …
‘String’,num2str(tspan(ii),3));
pause(0.001)
M = getframe(G);
writeVideo(myVideo, M);
end
end
end pde MATLAB Answers — New Questions
How to add a legend for a boxplot that indicates how the boxplot was created (summary statistics)?
Hi folks, I have a simple boxplot and I can’t figure out how to make a legend like the one shown in the photograph below. Ideally, the symbols and line specs would all match the associated text.
Perhaps doing it using the annotation or note tool? Was wondering if anyone has done this before. This type of formatting is a requirement for a journal paper.
For example (see example.png) I’ve gotten this far:
data = [1 2 3 4 4 5 5 6 6 7 8 9 13]
figure; boxplot(data);
a = get(get(gca,’children’),’children’); % Get the handles of all the objects
legend([a(1) a(2) a(3) a(4)],{‘Outliers’,’Median’,’25-75%’,’+/-1.5 IQR’})
But am wondering if there are alternative or better ways, and perhaps a way to show the blue bounding box? Just wanted to hear y’alls thoughts. Cheers.Hi folks, I have a simple boxplot and I can’t figure out how to make a legend like the one shown in the photograph below. Ideally, the symbols and line specs would all match the associated text.
Perhaps doing it using the annotation or note tool? Was wondering if anyone has done this before. This type of formatting is a requirement for a journal paper.
For example (see example.png) I’ve gotten this far:
data = [1 2 3 4 4 5 5 6 6 7 8 9 13]
figure; boxplot(data);
a = get(get(gca,’children’),’children’); % Get the handles of all the objects
legend([a(1) a(2) a(3) a(4)],{‘Outliers’,’Median’,’25-75%’,’+/-1.5 IQR’})
But am wondering if there are alternative or better ways, and perhaps a way to show the blue bounding box? Just wanted to hear y’alls thoughts. Cheers. Hi folks, I have a simple boxplot and I can’t figure out how to make a legend like the one shown in the photograph below. Ideally, the symbols and line specs would all match the associated text.
Perhaps doing it using the annotation or note tool? Was wondering if anyone has done this before. This type of formatting is a requirement for a journal paper.
For example (see example.png) I’ve gotten this far:
data = [1 2 3 4 4 5 5 6 6 7 8 9 13]
figure; boxplot(data);
a = get(get(gca,’children’),’children’); % Get the handles of all the objects
legend([a(1) a(2) a(3) a(4)],{‘Outliers’,’Median’,’25-75%’,’+/-1.5 IQR’})
But am wondering if there are alternative or better ways, and perhaps a way to show the blue bounding box? Just wanted to hear y’alls thoughts. Cheers. boxplot, legend, figure, labels MATLAB Answers — New Questions
Unable to save trainingstats in MATLAB 2024a
I am runnning a training script based on a simulink environment. At the end when I try to save the training results, it always throws an error. I am wondering if it has something to do with the .mat format. This issue does not occur in Matlab 2023a.
%% Train the agent
trainingStats = train(agent, env, trainOpts);
%% Save Agent
results_dir = strcat(‘./Results/’,num2str(time(1)),’0′,num2str(time(2)),num2str(time(3)));
mkdir(results_dir)
save(strcat(results_dir,’/TrainingStats.mat’),"trainingStats")
generatePolicyBlock(agent)
save_system(‘untitled’,’Agent’)
%% Simulate
%addpath(‘D:MastersHiWiTikz’)
simEpisodes = 1;
simOpts = rlSimulationOptions("MaxSteps",1250,…
"NumSimulations", simEpisodes);
experience = sim(env,agent,simOpts);
save(strcat(results_dir,’/Experience.mat’),"experience")
The error message is usually:
>> load(‘C:Usersvasu3DocumentsWorksharma_ma_rlonh2dfacados_implementationrlsavedAgentsPPO20240814Agent5852.mat’, ‘savedAgentResult’)
Error using load
Cannot read file
C:Usersvasu3DocumentsWorksharma_ma_rlonh2dfacados_implementationrlsavedAgentsPPO20240814Agent5852.mat.
Warning: Directory already exists.
Error using save
Unable to save file
‘C:Usersvasu3DocumentsWorksharma_ma_rlonh2dfacados_implementationrlResults20240814TrainingStats.mat’.
The file could not be closed, and might now be corrupt.
Error in run_H2DFRL_GRU (line 195)
save(strcat(results_dir,’/TrainingStats.mat’),"trainingStats")
Warning: Unable to read some of the variables due to unknown MAT-file error.
> In matfinfo (line 9)
In finfo>fetchDescriptions (line 278)
In finfo (line 82)
In uiimport/gatherFilePreviewData (line 419)
In uiimport (line 260)
Error using load
Unable to read file
‘C:Usersvasu3DocumentsWorksharma_ma_rlonh2dfacados_implementationrlResults20240814TrainingStats.mat’. Input must be a
MAT-file or an ASCII file containing numeric data with same number of columns in each row.
Error in uiimport/runImportdata (line 470)
datastruct = load(‘-ascii’, fileAbsolutePath);
Error in uiimport/gatherFilePreviewData (line 438)
[datastruct, textDelimiter, headerLines]= runImportdata(fileAbsolutePath, type);
Error in uiimport (line 260)
gatherFilePreviewData(fileAbsolutePath);
Could someone please help outI am runnning a training script based on a simulink environment. At the end when I try to save the training results, it always throws an error. I am wondering if it has something to do with the .mat format. This issue does not occur in Matlab 2023a.
%% Train the agent
trainingStats = train(agent, env, trainOpts);
%% Save Agent
results_dir = strcat(‘./Results/’,num2str(time(1)),’0′,num2str(time(2)),num2str(time(3)));
mkdir(results_dir)
save(strcat(results_dir,’/TrainingStats.mat’),"trainingStats")
generatePolicyBlock(agent)
save_system(‘untitled’,’Agent’)
%% Simulate
%addpath(‘D:MastersHiWiTikz’)
simEpisodes = 1;
simOpts = rlSimulationOptions("MaxSteps",1250,…
"NumSimulations", simEpisodes);
experience = sim(env,agent,simOpts);
save(strcat(results_dir,’/Experience.mat’),"experience")
The error message is usually:
>> load(‘C:Usersvasu3DocumentsWorksharma_ma_rlonh2dfacados_implementationrlsavedAgentsPPO20240814Agent5852.mat’, ‘savedAgentResult’)
Error using load
Cannot read file
C:Usersvasu3DocumentsWorksharma_ma_rlonh2dfacados_implementationrlsavedAgentsPPO20240814Agent5852.mat.
Warning: Directory already exists.
Error using save
Unable to save file
‘C:Usersvasu3DocumentsWorksharma_ma_rlonh2dfacados_implementationrlResults20240814TrainingStats.mat’.
The file could not be closed, and might now be corrupt.
Error in run_H2DFRL_GRU (line 195)
save(strcat(results_dir,’/TrainingStats.mat’),"trainingStats")
Warning: Unable to read some of the variables due to unknown MAT-file error.
> In matfinfo (line 9)
In finfo>fetchDescriptions (line 278)
In finfo (line 82)
In uiimport/gatherFilePreviewData (line 419)
In uiimport (line 260)
Error using load
Unable to read file
‘C:Usersvasu3DocumentsWorksharma_ma_rlonh2dfacados_implementationrlResults20240814TrainingStats.mat’. Input must be a
MAT-file or an ASCII file containing numeric data with same number of columns in each row.
Error in uiimport/runImportdata (line 470)
datastruct = load(‘-ascii’, fileAbsolutePath);
Error in uiimport/gatherFilePreviewData (line 438)
[datastruct, textDelimiter, headerLines]= runImportdata(fileAbsolutePath, type);
Error in uiimport (line 260)
gatherFilePreviewData(fileAbsolutePath);
Could someone please help out I am runnning a training script based on a simulink environment. At the end when I try to save the training results, it always throws an error. I am wondering if it has something to do with the .mat format. This issue does not occur in Matlab 2023a.
%% Train the agent
trainingStats = train(agent, env, trainOpts);
%% Save Agent
results_dir = strcat(‘./Results/’,num2str(time(1)),’0′,num2str(time(2)),num2str(time(3)));
mkdir(results_dir)
save(strcat(results_dir,’/TrainingStats.mat’),"trainingStats")
generatePolicyBlock(agent)
save_system(‘untitled’,’Agent’)
%% Simulate
%addpath(‘D:MastersHiWiTikz’)
simEpisodes = 1;
simOpts = rlSimulationOptions("MaxSteps",1250,…
"NumSimulations", simEpisodes);
experience = sim(env,agent,simOpts);
save(strcat(results_dir,’/Experience.mat’),"experience")
The error message is usually:
>> load(‘C:Usersvasu3DocumentsWorksharma_ma_rlonh2dfacados_implementationrlsavedAgentsPPO20240814Agent5852.mat’, ‘savedAgentResult’)
Error using load
Cannot read file
C:Usersvasu3DocumentsWorksharma_ma_rlonh2dfacados_implementationrlsavedAgentsPPO20240814Agent5852.mat.
Warning: Directory already exists.
Error using save
Unable to save file
‘C:Usersvasu3DocumentsWorksharma_ma_rlonh2dfacados_implementationrlResults20240814TrainingStats.mat’.
The file could not be closed, and might now be corrupt.
Error in run_H2DFRL_GRU (line 195)
save(strcat(results_dir,’/TrainingStats.mat’),"trainingStats")
Warning: Unable to read some of the variables due to unknown MAT-file error.
> In matfinfo (line 9)
In finfo>fetchDescriptions (line 278)
In finfo (line 82)
In uiimport/gatherFilePreviewData (line 419)
In uiimport (line 260)
Error using load
Unable to read file
‘C:Usersvasu3DocumentsWorksharma_ma_rlonh2dfacados_implementationrlResults20240814TrainingStats.mat’. Input must be a
MAT-file or an ASCII file containing numeric data with same number of columns in each row.
Error in uiimport/runImportdata (line 470)
datastruct = load(‘-ascii’, fileAbsolutePath);
Error in uiimport/gatherFilePreviewData (line 438)
[datastruct, textDelimiter, headerLines]= runImportdata(fileAbsolutePath, type);
Error in uiimport (line 260)
gatherFilePreviewData(fileAbsolutePath);
Could someone please help out reinforcement learning., deep learning, machine learning, bug, matlab MATLAB Answers — New Questions
Cannot compile app with Segment Anything Model SAM due to license?
Hi,
I could not compile my app with SAM due to the following error.
Where can I find the specific LICENSE that states why and how SAM cannot be packaged?
Or am I misunderstanding anything about specifying parameters for packaging? I’m very new to compiling app in Matlab.
Thanks,
Lam
mcc -o v9 -W ‘WinMain:v9,version=1.0’ -T link:exe -d ‘D:OneDrive – ZeonCorporation1.WORK24_MatlabSAMworkedv9for_testing’ -v ‘D:OneDrive – ZeonCorporation1.WORK24_MatlabSAMworkedv9.m’ -a ‘D:OneDrive – ZeonCorporation1.WORK24_MatlabSAMworkedLICENSE’ -a C:ProgramDataMATLABSupportPackagesR2024atoolboximagessupportpackagessamdatapreTrainedSAM.mat -a C:ProgramDataMATLABSupportPackagesR2024atoolboximagessupportpackagessam+matlabshared+supportpkg+internal+sppkglegacySAM.m -r ‘C:Program FilesMATLABR2024atoolboxcompilerpackagingResourcesdefault_icon.ico’ -Z ‘Deep Learning Toolbox Converter for ONNX Model Format’
Compiler version: 24.1 (R2024a)
Analyzing file dependencies.
Warning: Removed http proxy service credentials from preference settings.
Warning: In "C:ProgramDataMATLABSupportPackagesR2024atoolboximagessupportpackagessam+matlabshared+supportpkg+internal+sppkglegacySAM.m", "matlabshared.supportpkg.internal.sppkglegacy.SupportPackageRegistryPluginBase" are excluded from packaging for the MATLAB Runtime environment according to the MATLAB Compiler license. Either remove the file or function from your code, or use the MATLAB function "isdeployed" to ensure the function is not invoked in the deployed component.
foundation::storage::vfs::Exception
mcc failed.Hi,
I could not compile my app with SAM due to the following error.
Where can I find the specific LICENSE that states why and how SAM cannot be packaged?
Or am I misunderstanding anything about specifying parameters for packaging? I’m very new to compiling app in Matlab.
Thanks,
Lam
mcc -o v9 -W ‘WinMain:v9,version=1.0’ -T link:exe -d ‘D:OneDrive – ZeonCorporation1.WORK24_MatlabSAMworkedv9for_testing’ -v ‘D:OneDrive – ZeonCorporation1.WORK24_MatlabSAMworkedv9.m’ -a ‘D:OneDrive – ZeonCorporation1.WORK24_MatlabSAMworkedLICENSE’ -a C:ProgramDataMATLABSupportPackagesR2024atoolboximagessupportpackagessamdatapreTrainedSAM.mat -a C:ProgramDataMATLABSupportPackagesR2024atoolboximagessupportpackagessam+matlabshared+supportpkg+internal+sppkglegacySAM.m -r ‘C:Program FilesMATLABR2024atoolboxcompilerpackagingResourcesdefault_icon.ico’ -Z ‘Deep Learning Toolbox Converter for ONNX Model Format’
Compiler version: 24.1 (R2024a)
Analyzing file dependencies.
Warning: Removed http proxy service credentials from preference settings.
Warning: In "C:ProgramDataMATLABSupportPackagesR2024atoolboximagessupportpackagessam+matlabshared+supportpkg+internal+sppkglegacySAM.m", "matlabshared.supportpkg.internal.sppkglegacy.SupportPackageRegistryPluginBase" are excluded from packaging for the MATLAB Runtime environment according to the MATLAB Compiler license. Either remove the file or function from your code, or use the MATLAB function "isdeployed" to ensure the function is not invoked in the deployed component.
foundation::storage::vfs::Exception
mcc failed. Hi,
I could not compile my app with SAM due to the following error.
Where can I find the specific LICENSE that states why and how SAM cannot be packaged?
Or am I misunderstanding anything about specifying parameters for packaging? I’m very new to compiling app in Matlab.
Thanks,
Lam
mcc -o v9 -W ‘WinMain:v9,version=1.0’ -T link:exe -d ‘D:OneDrive – ZeonCorporation1.WORK24_MatlabSAMworkedv9for_testing’ -v ‘D:OneDrive – ZeonCorporation1.WORK24_MatlabSAMworkedv9.m’ -a ‘D:OneDrive – ZeonCorporation1.WORK24_MatlabSAMworkedLICENSE’ -a C:ProgramDataMATLABSupportPackagesR2024atoolboximagessupportpackagessamdatapreTrainedSAM.mat -a C:ProgramDataMATLABSupportPackagesR2024atoolboximagessupportpackagessam+matlabshared+supportpkg+internal+sppkglegacySAM.m -r ‘C:Program FilesMATLABR2024atoolboxcompilerpackagingResourcesdefault_icon.ico’ -Z ‘Deep Learning Toolbox Converter for ONNX Model Format’
Compiler version: 24.1 (R2024a)
Analyzing file dependencies.
Warning: Removed http proxy service credentials from preference settings.
Warning: In "C:ProgramDataMATLABSupportPackagesR2024atoolboximagessupportpackagessam+matlabshared+supportpkg+internal+sppkglegacySAM.m", "matlabshared.supportpkg.internal.sppkglegacy.SupportPackageRegistryPluginBase" are excluded from packaging for the MATLAB Runtime environment according to the MATLAB Compiler license. Either remove the file or function from your code, or use the MATLAB function "isdeployed" to ensure the function is not invoked in the deployed component.
foundation::storage::vfs::Exception
mcc failed. segment anything model, sam, compiler, matlab compiler MATLAB Answers — New Questions
Parameter influence estimation using monte carlo
Hi,
I have 14 parameters of groundwater, Ca, Mg, pH, Cl, NO3, TD, EC.. etc. I want to estimate the influence of each parameters for groundwater quality using monte carlo but had no idea about the coding. I try to look out for some coding but really didnt get it.
Please helpHi,
I have 14 parameters of groundwater, Ca, Mg, pH, Cl, NO3, TD, EC.. etc. I want to estimate the influence of each parameters for groundwater quality using monte carlo but had no idea about the coding. I try to look out for some coding but really didnt get it.
Please help Hi,
I have 14 parameters of groundwater, Ca, Mg, pH, Cl, NO3, TD, EC.. etc. I want to estimate the influence of each parameters for groundwater quality using monte carlo but had no idea about the coding. I try to look out for some coding but really didnt get it.
Please help monte carlo MATLAB Answers — New Questions
How to setup system path to vivado path?
While using hdlsetuptoolpath the system path is being prepended twice.
First the correct path and second time a non existent path.
Due to this I can’t validate the model it says Xilinx Vivado not found.While using hdlsetuptoolpath the system path is being prepended twice.
First the correct path and second time a non existent path.
Due to this I can’t validate the model it says Xilinx Vivado not found. While using hdlsetuptoolpath the system path is being prepended twice.
First the correct path and second time a non existent path.
Due to this I can’t validate the model it says Xilinx Vivado not found. hdlsetuptoolpath, soc builder MATLAB Answers — New Questions
sym to double data type:
I want to obtain the y_values in double but they are in symbolic form. Anyway I can get around this? Thank you!
L_e = 0.7199;
Y_poly = @(x) a4*x.^4 + a3*x.^3 + a2*x.^2 + a1*x + a0;
x_values = linspace(0, L_e, 100); % 100 points from 0 to L_e
y_values = Y_poly(x_values);
The values for the coefficients are:
a4 = -0.324785
a3 = 0.472768
a2 = -0.011102
a1 = 0.000000
a0 = 1.000000I want to obtain the y_values in double but they are in symbolic form. Anyway I can get around this? Thank you!
L_e = 0.7199;
Y_poly = @(x) a4*x.^4 + a3*x.^3 + a2*x.^2 + a1*x + a0;
x_values = linspace(0, L_e, 100); % 100 points from 0 to L_e
y_values = Y_poly(x_values);
The values for the coefficients are:
a4 = -0.324785
a3 = 0.472768
a2 = -0.011102
a1 = 0.000000
a0 = 1.000000 I want to obtain the y_values in double but they are in symbolic form. Anyway I can get around this? Thank you!
L_e = 0.7199;
Y_poly = @(x) a4*x.^4 + a3*x.^3 + a2*x.^2 + a1*x + a0;
x_values = linspace(0, L_e, 100); % 100 points from 0 to L_e
y_values = Y_poly(x_values);
The values for the coefficients are:
a4 = -0.324785
a3 = 0.472768
a2 = -0.011102
a1 = 0.000000
a0 = 1.000000 #sym, #double, #polynomial MATLAB Answers — New Questions
Issue with YOLO4
Hi,
I tried to run the yolo code , but i get the error below :
‘validateInputData’ is used in the following examples:
Object Detection Using YOLO v3 Deep Learning
Object Detection Using YOLO v4 Deep Learning
Prune Filters in a Detection Network Using Taylor Scores
Error in yolo4Test (line 45)
validateInputData(trainingData);
Could you help please ?
Iam trying to run the Mathworks code for yolo4 ?Hi,
I tried to run the yolo code , but i get the error below :
‘validateInputData’ is used in the following examples:
Object Detection Using YOLO v3 Deep Learning
Object Detection Using YOLO v4 Deep Learning
Prune Filters in a Detection Network Using Taylor Scores
Error in yolo4Test (line 45)
validateInputData(trainingData);
Could you help please ?
Iam trying to run the Mathworks code for yolo4 ? Hi,
I tried to run the yolo code , but i get the error below :
‘validateInputData’ is used in the following examples:
Object Detection Using YOLO v3 Deep Learning
Object Detection Using YOLO v4 Deep Learning
Prune Filters in a Detection Network Using Taylor Scores
Error in yolo4Test (line 45)
validateInputData(trainingData);
Could you help please ?
Iam trying to run the Mathworks code for yolo4 ? deep learning, yolo4 MATLAB Answers — New Questions
Error: variable-size matrix type is not supported for HDL code
I use the fixed-point tool to fixed-point the subsystem and then generate Verilog, but the following error occurs.
The model contains constructs that are unsupported for HDL code generation. HDL Coder ‘c’ : Error: variable-size matrix type is not supported for HDL code generation. Function ’eml_fixpt_times’ (#33554529.1887.1910), line 65, column 5 Function ‘times’ (#33554530.5290.5318), line 146, column 27 Function ‘mtimes’ (#33554528.2252.2264), line 62, column 9 Function ‘forSubsystem/MATLAB Function1/MATLAB Function1_FixPt’ (#369.531.542), line 23, column 13
In #369.531.542, the error is related to miu*x*en(i), i tried to use sss=miu*x*en(i) instead, but it ended up showing Subscripted assignment dimension mismatch: [1] ~= [5]. Error in ‘testfixed/testforSubsystem/MATLAB Function’ (line 23) sss = miu*x*en(i);
I fixed the above problem,
sss= zeros(1,5);
sss = miu*x*en(i);
wn(:)=wn+sss;
but still failed to generate Verilog,
The model contains constructs that are unsupported for HDL code generation. HDL Coder ‘c’ : Error: variable-size matrix type is not supported for HDL code generation. Function ’eml_fixpt_times’ (#33554529.1887.1910), line 65, column 5 Function ‘times’ (#33554530.5290.5318), line 146, column 27 Function ‘mtimes’ (#33554528.2252.2264), line 62, column 9 Function ‘forsubsystem/MATLAB Function’ (#998.580.591), line 23, column 1
The key to the error is on en(i), why does it cause this error?I use the fixed-point tool to fixed-point the subsystem and then generate Verilog, but the following error occurs.
The model contains constructs that are unsupported for HDL code generation. HDL Coder ‘c’ : Error: variable-size matrix type is not supported for HDL code generation. Function ’eml_fixpt_times’ (#33554529.1887.1910), line 65, column 5 Function ‘times’ (#33554530.5290.5318), line 146, column 27 Function ‘mtimes’ (#33554528.2252.2264), line 62, column 9 Function ‘forSubsystem/MATLAB Function1/MATLAB Function1_FixPt’ (#369.531.542), line 23, column 13
In #369.531.542, the error is related to miu*x*en(i), i tried to use sss=miu*x*en(i) instead, but it ended up showing Subscripted assignment dimension mismatch: [1] ~= [5]. Error in ‘testfixed/testforSubsystem/MATLAB Function’ (line 23) sss = miu*x*en(i);
I fixed the above problem,
sss= zeros(1,5);
sss = miu*x*en(i);
wn(:)=wn+sss;
but still failed to generate Verilog,
The model contains constructs that are unsupported for HDL code generation. HDL Coder ‘c’ : Error: variable-size matrix type is not supported for HDL code generation. Function ’eml_fixpt_times’ (#33554529.1887.1910), line 65, column 5 Function ‘times’ (#33554530.5290.5318), line 146, column 27 Function ‘mtimes’ (#33554528.2252.2264), line 62, column 9 Function ‘forsubsystem/MATLAB Function’ (#998.580.591), line 23, column 1
The key to the error is on en(i), why does it cause this error? I use the fixed-point tool to fixed-point the subsystem and then generate Verilog, but the following error occurs.
The model contains constructs that are unsupported for HDL code generation. HDL Coder ‘c’ : Error: variable-size matrix type is not supported for HDL code generation. Function ’eml_fixpt_times’ (#33554529.1887.1910), line 65, column 5 Function ‘times’ (#33554530.5290.5318), line 146, column 27 Function ‘mtimes’ (#33554528.2252.2264), line 62, column 9 Function ‘forSubsystem/MATLAB Function1/MATLAB Function1_FixPt’ (#369.531.542), line 23, column 13
In #369.531.542, the error is related to miu*x*en(i), i tried to use sss=miu*x*en(i) instead, but it ended up showing Subscripted assignment dimension mismatch: [1] ~= [5]. Error in ‘testfixed/testforSubsystem/MATLAB Function’ (line 23) sss = miu*x*en(i);
I fixed the above problem,
sss= zeros(1,5);
sss = miu*x*en(i);
wn(:)=wn+sss;
but still failed to generate Verilog,
The model contains constructs that are unsupported for HDL code generation. HDL Coder ‘c’ : Error: variable-size matrix type is not supported for HDL code generation. Function ’eml_fixpt_times’ (#33554529.1887.1910), line 65, column 5 Function ‘times’ (#33554530.5290.5318), line 146, column 27 Function ‘mtimes’ (#33554528.2252.2264), line 62, column 9 Function ‘forsubsystem/MATLAB Function’ (#998.580.591), line 23, column 1
The key to the error is on en(i), why does it cause this error? simulink model, variable-size, generate code MATLAB Answers — New Questions
quadprog output: this problem is non-convex
I am trying to solve a quadratic optimization problem but quadprog keeps telling me that my problem is non-convex.
After several experiments, I found that the problem comes from the equation constraints matrix A, which is a 57250*57441 matrix.
For the following code,
[m, n] = size(A);
assert(m < n);
options = optimoptions(‘quadprog’,’Display’,’off’);
[Pwp,fval,exitflag,output] = quadprog(speye(n), zeros(n,1), [], [], A, zeros(m, 1), [], [], [], options);
obviously the solution should be the all-zero vector. But the output still said that this is a nonconvex problem.I am trying to solve a quadratic optimization problem but quadprog keeps telling me that my problem is non-convex.
After several experiments, I found that the problem comes from the equation constraints matrix A, which is a 57250*57441 matrix.
For the following code,
[m, n] = size(A);
assert(m < n);
options = optimoptions(‘quadprog’,’Display’,’off’);
[Pwp,fval,exitflag,output] = quadprog(speye(n), zeros(n,1), [], [], A, zeros(m, 1), [], [], [], options);
obviously the solution should be the all-zero vector. But the output still said that this is a nonconvex problem. I am trying to solve a quadratic optimization problem but quadprog keeps telling me that my problem is non-convex.
After several experiments, I found that the problem comes from the equation constraints matrix A, which is a 57250*57441 matrix.
For the following code,
[m, n] = size(A);
assert(m < n);
options = optimoptions(‘quadprog’,’Display’,’off’);
[Pwp,fval,exitflag,output] = quadprog(speye(n), zeros(n,1), [], [], A, zeros(m, 1), [], [], [], options);
obviously the solution should be the all-zero vector. But the output still said that this is a nonconvex problem. optimization MATLAB Answers — New Questions
How to change in marker size in the global legend?
Hi,
I’m making a plot containing a few subplots using the function tiledlayout, and I created a global legend using the code
leg = legend({‘A’,’B’,’C’})
leg.Layout.Tile = ‘North’
However with this I cann’t use the previous method to change the marker size in the legend, because it requirs two outputs from the legend, and it will override the previous code.
[~,icons] = legend({‘A’,’B’,’C’})
icons1=findobj(icons,’type’,’patch’);
set(icons1,’MarkerSize’,15,’Linewidth’,1.5);
Anyone know the workaround of this? many thanks!Hi,
I’m making a plot containing a few subplots using the function tiledlayout, and I created a global legend using the code
leg = legend({‘A’,’B’,’C’})
leg.Layout.Tile = ‘North’
However with this I cann’t use the previous method to change the marker size in the legend, because it requirs two outputs from the legend, and it will override the previous code.
[~,icons] = legend({‘A’,’B’,’C’})
icons1=findobj(icons,’type’,’patch’);
set(icons1,’MarkerSize’,15,’Linewidth’,1.5);
Anyone know the workaround of this? many thanks! Hi,
I’m making a plot containing a few subplots using the function tiledlayout, and I created a global legend using the code
leg = legend({‘A’,’B’,’C’})
leg.Layout.Tile = ‘North’
However with this I cann’t use the previous method to change the marker size in the legend, because it requirs two outputs from the legend, and it will override the previous code.
[~,icons] = legend({‘A’,’B’,’C’})
icons1=findobj(icons,’type’,’patch’);
set(icons1,’MarkerSize’,15,’Linewidth’,1.5);
Anyone know the workaround of this? many thanks! tiledlayout, makersize, global legend MATLAB Answers — New Questions
Computing gradient of mean curvature on a mesh using gp toolbox, unexpected error. What am I missing?
I am trying to compute this quantity here H grad(H) where H is the mean curvture
The code breaks at the last line because there is a dimension mismatch between H and grad(H). size(grad(H))= [960 3] and size(V)=[162 3] and size(mean_curvture)=[162 3]. Why is the grad(H) has that dimension? where is the mistake?
% Load a mesh
[V, F] = load_mesh(‘sphere.off’);
% Compute Mean Curvature and Normal Vector
% 2. Compute Cotangent Laplacian (L) and Mass Matrix (M)
L = cotmatrix(V, F);
M = massmatrix(V, F, ‘barycentric’);
% 3. Compute Mean Curvature Vector (H)
H = -inv(M) * (L * V);
% 4. Compute the Magnitude of Mean Curvature (mean curvature at each vertex)
mean_curvature = sqrt(sum(H.^2, 2));
% Compute the Normal Vector
normals = per_vertex_normals(V, F);
% Compute Gaussian Curvature
gaussian_curvature = discrete_gaussian_curvature(V,F);
% Compute Laplacian of Mean Curvature
laplacian_H = L * mean_curvature;
% Compute the Gradient of Mean Curvature
% Use gptoolbox’s gradient operator for scalar fields
grad_H = grad(V, F) * mean_curvature;
% Compute the new term: newthing
%H_squared = mean_curvature .^ 2;
%newthing = laplacian_H + ((H_squared – 2 * gaussian_curvature) .* mean_curvature) .* normals + mean_curvature .* grad_H;
% Compute the new term: newthing
H_squared = mean_curvature .^ 2;
newthing = laplacian_H + ((H_squared – 2 * gaussian_curvature) .* mean_curvature) .* normals + mean_curvature .* grad_H;
% Define small arc equation
h = 0.1; % Small parameter h
vertex = V; % V contains the vertices
% Compute the small arc
f_h = vertex + (mean_curvature .* normals) * h + newthing * (h^2 / 2);
% Visualization of small arc
figure;
trisurf(F, V(:,1), V(:,2), V(:,3), ‘FaceColor’, ‘cyan’, ‘EdgeColor’, ‘none’);
hold on;
plot3(f_h(:,1), f_h(:,2), f_h(:,3), ‘r.’, ‘MarkerSize’, 10);
axis equal;
title(‘Small Arc Visualization’);
xlabel(‘X’);
ylabel(‘Y’);
zlabel(‘Z’);
the grad(H) script
% Step 1: Load a Mesh
[V, F] = load_mesh(‘sphere.off’);
% Step 2: Compute Laplace-Beltrami Operator
L = cotmatrix(V, F); % cotangent Laplace-Beltrami operator
% Step 3: Compute the Mean Curvature Vector (H)
H = -L * V; % Mean curvature normal vector
% Step 4: Compute the Mean Curvature Magnitude
mean_curvature = sqrt(sum(H.^2, 2));
% Step 5: Compute the Gradient of Mean Curvature
% Use gptoolbox’s gradient operator for scalar fields
grad_H = grad(V, F) * mean_curvature;
% Step 6: Visualization or Further Processing
figure;
trisurf(F, V(:,1), V(:,2), V(:,3), mean_curvature, ‘EdgeColor’, ‘none’);
axis equal;
lighting gouraud;
camlight;
colorbar;
title(‘Mean Curvature Gradient’);
xlabel(‘X’);
ylabel(‘Y’);
zlabel(‘Z’);
grad(V, F) from gp toolbox
function [G] = grad(V,F)
% GRAD Compute the numerical gradient operator for triangle or tet meshes
%
% G = grad(V,F)
%
% Inputs:
% V #vertices by dim list of mesh vertex positions
% F #faces by simplex-size list of mesh face indices
% Outputs:
% G #faces*dim by #V Gradient operator
%
% Example:
% L = cotmatrix(V,F);
% G = grad(V,F);
% dblA = doublearea(V,F);
% GMG = -G’*repdiag(diag(sparse(dblA)/2),size(V,2))*G;
%
% % Columns of W are scalar fields
% G = grad(V,F);
% % Compute gradient magnitude for each column in W
% GM = squeeze(sqrt(sum(reshape(G*W,size(F,1),size(V,2),size(W,2)).^2,2)));
%
dim = size(V,2);
ss = size(F,2);
switch ss
case 2
% Edge lengths
len = normrow(V(F(:,2),:)-V(F(:,1),:));
% Gradient is just staggered grid finite difference
G = sparse(repmat(1:size(F,1),2,1)’,F,[1 -1]./len, size(F,1),size(V,1));
case 3
% append with 0s for convenience
if size(V,2) == 2
V = [V zeros(size(V,1),1)];
end
% Gradient of a scalar function defined on piecewise linear elements (mesh)
% is constant on each triangle i,j,k:
% grad(Xijk) = (Xj-Xi) * (Vi – Vk)^R90 / 2A + (Xk-Xi) * (Vj – Vi)^R90 / 2A
% grad(Xijk) = Xj * (Vi – Vk)^R90 / 2A + Xk * (Vj – Vi)^R90 / 2A +
% -Xi * (Vi – Vk)^R90 / 2A – Xi * (Vj – Vi)^R90 / 2A
% where Xi is the scalar value at vertex i, Vi is the 3D position of vertex
% i, and A is the area of triangle (i,j,k). ^R90 represent a rotation of
% 90 degrees
%
% renaming indices of vertices of triangles for convenience
i1 = F(:,1); i2 = F(:,2); i3 = F(:,3);
% #F x 3 matrices of triangle edge vectors, named after opposite vertices
v32 = V(i3,:) – V(i2,:); v13 = V(i1,:) – V(i3,:); v21 = V(i2,:) – V(i1,:);
% area of parallelogram is twice area of triangle
% area of parallelogram is || v1 x v2 ||
n = cross(v32,v13,2);
% This does correct l2 norm of rows, so that it contains #F list of twice
% triangle areas
dblA = normrow(n);
% now normalize normals to get unit normals
u = normalizerow(n);
% rotate each vector 90 degrees around normal
%eperp21 = bsxfun(@times,normalizerow(cross(u,v21)),normrow(v21)./dblA);
%eperp13 = bsxfun(@times,normalizerow(cross(u,v13)),normrow(v13)./dblA);
eperp21 = bsxfun(@times,cross(u,v21),1./dblA);
eperp13 = bsxfun(@times,cross(u,v13),1./dblA);
%g = …
% ( …
% repmat(X(F(:,2)) – X(F(:,1)),1,3).*eperp13 + …
% repmat(X(F(:,3)) – X(F(:,1)),1,3).*eperp21 …
% );
GI = …
[0*size(F,1)+repmat(1:size(F,1),1,4) …
1*size(F,1)+repmat(1:size(F,1),1,4) …
2*size(F,1)+repmat(1:size(F,1),1,4)]’;
GJ = repmat([F(:,2);F(:,1);F(:,3);F(:,1)],3,1);
GV = [eperp13(:,1);-eperp13(:,1);eperp21(:,1);-eperp21(:,1); …
eperp13(:,2);-eperp13(:,2);eperp21(:,2);-eperp21(:,2); …
eperp13(:,3);-eperp13(:,3);eperp21(:,3);-eperp21(:,3)];
G = sparse(GI,GJ,GV, 3*size(F,1), size(V,1));
%% Alternatively
%%
%% f(x) is piecewise-linear function:
%%
%% f(x) = ∑ φi(x) fi, f(x ∈ T) = φi(x) fi + φj(x) fj + φk(x) fk
%% ∇f(x) = … = ∇φi(x) fi + ∇φj(x) fj + ∇φk(x) fk
%% = ∇φi fi + ∇φj fj + ∇φk) fk
%%
%% ∇φi = 1/hjk ((Vj-Vk)/||Vj-Vk||)^perp =
%% = ||Vj-Vk|| /(2 Aijk) * ((Vj-Vk)/||Vj-Vk||)^perp
%% = 1/(2 Aijk) * (Vj-Vk)^perp
%%
%m = size(F,1);
%eperp32 = bsxfun(@times,cross(u,v32),1./dblA);
%G = sparse( …
% [0*m + repmat(1:m,1,3) …
% 1*m + repmat(1:m,1,3) …
% 2*m + repmat(1:m,1,3)]’, …
% repmat([F(:,1);F(:,2);F(:,3)],3,1), …
% [eperp32(:,1);eperp13(:,1);eperp21(:,1); …
% eperp32(:,2);eperp13(:,2);eperp21(:,2); …
% eperp32(:,3);eperp13(:,3);eperp21(:,3)], …
% 3*m,size(V,1));
if dim == 2
G = G(1:(size(F,1)*dim),:);
end
% Should be the same as:
% g = …
% bsxfun(@times,X(F(:,1)),cross(u,v32)) + …
% bsxfun(@times,X(F(:,2)),cross(u,v13)) + …
% bsxfun(@times,X(F(:,3)),cross(u,v21));
% g = bsxfun(@rdivide,g,dblA);
case 4
% really dealing with tets
T = F;
% number of dimensions
assert(dim == 3);
% number of vertices
n = size(V,1);
% number of elements
m = size(T,1);
% simplex size
assert(size(T,2) == 4);
if m == 1 && ~isnumeric(V)
simple_volume = @(ad,r) -sum(ad.*r,2)./6;
simple_volume = @(ad,bd,cd) simple_volume(ad, …
[bd(:,2).*cd(:,3)-bd(:,3).*cd(:,2), …
bd(:,3).*cd(:,1)-bd(:,1).*cd(:,3), …
bd(:,1).*cd(:,2)-bd(:,2).*cd(:,1)]);
P = sym(‘P’,[1 3]);
V1 = V(T(:,1),:);
V2 = V(T(:,2),:);
V3 = V(T(:,3),:);
V4 = V(T(:,4),:);
V1P = V1-P;
V2P = V2-P;
V3P = V3-P;
V4P = V4-P;
A1 = simple_volume(V2P,V4P,V3P);
A2 = simple_volume(V1P,V3P,V4P);
A3 = simple_volume(V1P,V4P,V2P);
A4 = simple_volume(V1P,V2P,V3P);
B = [A1 A2 A3 A4]./(A1+A2+A3+A4);
G = simplify(jacobian(B,P).’);
return
end
% f(x) is piecewise-linear function:
%
% f(x) = ∑ φi(x) fi, f(x ∈ T) = φi(x) fi + φj(x) fj + φk(x) fk + φl(x) fl
% ∇f(x) = … = ∇φi(x) fi + ∇φj(x) fj + ∇φk(x) fk + ∇φl(x) fl
% = ∇φi fi + ∇φj fj + ∇φk fk + ∇φl fl
%
% ∇φi = 1/hjk = Ajkl / 3V * (Facejkl)^perp
% = Ajkl / 3V * (Vj-Vk)x(Vl-Vk)
% = Ajkl / 3V * Njkl / ||Njkl||
%
% get all faces
F = [ …
T(:,1) T(:,2) T(:,3); …
T(:,1) T(:,3) T(:,4); …
T(:,1) T(:,4) T(:,2); …
T(:,2) T(:,4) T(:,3)];
% compute areas of each face
A = doublearea(V,F)/2;
N = normalizerow(normals(V,F));
% compute volume of each tet
vol = volume(V,T);
GI = …
[0*m + repmat(1:m,1,4) …
1*m + repmat(1:m,1,4) …
2*m + repmat(1:m,1,4)];
GJ = repmat([T(:,4);T(:,2);T(:,3);T(:,1)],3,1);
GV = repmat(A./(3*repmat(vol,4,1)),3,1).*N(:);
G = sparse(GI,GJ,GV, 3*m,n);
end
end
the sphere.off file (mesh data)I am trying to compute this quantity here H grad(H) where H is the mean curvture
The code breaks at the last line because there is a dimension mismatch between H and grad(H). size(grad(H))= [960 3] and size(V)=[162 3] and size(mean_curvture)=[162 3]. Why is the grad(H) has that dimension? where is the mistake?
% Load a mesh
[V, F] = load_mesh(‘sphere.off’);
% Compute Mean Curvature and Normal Vector
% 2. Compute Cotangent Laplacian (L) and Mass Matrix (M)
L = cotmatrix(V, F);
M = massmatrix(V, F, ‘barycentric’);
% 3. Compute Mean Curvature Vector (H)
H = -inv(M) * (L * V);
% 4. Compute the Magnitude of Mean Curvature (mean curvature at each vertex)
mean_curvature = sqrt(sum(H.^2, 2));
% Compute the Normal Vector
normals = per_vertex_normals(V, F);
% Compute Gaussian Curvature
gaussian_curvature = discrete_gaussian_curvature(V,F);
% Compute Laplacian of Mean Curvature
laplacian_H = L * mean_curvature;
% Compute the Gradient of Mean Curvature
% Use gptoolbox’s gradient operator for scalar fields
grad_H = grad(V, F) * mean_curvature;
% Compute the new term: newthing
%H_squared = mean_curvature .^ 2;
%newthing = laplacian_H + ((H_squared – 2 * gaussian_curvature) .* mean_curvature) .* normals + mean_curvature .* grad_H;
% Compute the new term: newthing
H_squared = mean_curvature .^ 2;
newthing = laplacian_H + ((H_squared – 2 * gaussian_curvature) .* mean_curvature) .* normals + mean_curvature .* grad_H;
% Define small arc equation
h = 0.1; % Small parameter h
vertex = V; % V contains the vertices
% Compute the small arc
f_h = vertex + (mean_curvature .* normals) * h + newthing * (h^2 / 2);
% Visualization of small arc
figure;
trisurf(F, V(:,1), V(:,2), V(:,3), ‘FaceColor’, ‘cyan’, ‘EdgeColor’, ‘none’);
hold on;
plot3(f_h(:,1), f_h(:,2), f_h(:,3), ‘r.’, ‘MarkerSize’, 10);
axis equal;
title(‘Small Arc Visualization’);
xlabel(‘X’);
ylabel(‘Y’);
zlabel(‘Z’);
the grad(H) script
% Step 1: Load a Mesh
[V, F] = load_mesh(‘sphere.off’);
% Step 2: Compute Laplace-Beltrami Operator
L = cotmatrix(V, F); % cotangent Laplace-Beltrami operator
% Step 3: Compute the Mean Curvature Vector (H)
H = -L * V; % Mean curvature normal vector
% Step 4: Compute the Mean Curvature Magnitude
mean_curvature = sqrt(sum(H.^2, 2));
% Step 5: Compute the Gradient of Mean Curvature
% Use gptoolbox’s gradient operator for scalar fields
grad_H = grad(V, F) * mean_curvature;
% Step 6: Visualization or Further Processing
figure;
trisurf(F, V(:,1), V(:,2), V(:,3), mean_curvature, ‘EdgeColor’, ‘none’);
axis equal;
lighting gouraud;
camlight;
colorbar;
title(‘Mean Curvature Gradient’);
xlabel(‘X’);
ylabel(‘Y’);
zlabel(‘Z’);
grad(V, F) from gp toolbox
function [G] = grad(V,F)
% GRAD Compute the numerical gradient operator for triangle or tet meshes
%
% G = grad(V,F)
%
% Inputs:
% V #vertices by dim list of mesh vertex positions
% F #faces by simplex-size list of mesh face indices
% Outputs:
% G #faces*dim by #V Gradient operator
%
% Example:
% L = cotmatrix(V,F);
% G = grad(V,F);
% dblA = doublearea(V,F);
% GMG = -G’*repdiag(diag(sparse(dblA)/2),size(V,2))*G;
%
% % Columns of W are scalar fields
% G = grad(V,F);
% % Compute gradient magnitude for each column in W
% GM = squeeze(sqrt(sum(reshape(G*W,size(F,1),size(V,2),size(W,2)).^2,2)));
%
dim = size(V,2);
ss = size(F,2);
switch ss
case 2
% Edge lengths
len = normrow(V(F(:,2),:)-V(F(:,1),:));
% Gradient is just staggered grid finite difference
G = sparse(repmat(1:size(F,1),2,1)’,F,[1 -1]./len, size(F,1),size(V,1));
case 3
% append with 0s for convenience
if size(V,2) == 2
V = [V zeros(size(V,1),1)];
end
% Gradient of a scalar function defined on piecewise linear elements (mesh)
% is constant on each triangle i,j,k:
% grad(Xijk) = (Xj-Xi) * (Vi – Vk)^R90 / 2A + (Xk-Xi) * (Vj – Vi)^R90 / 2A
% grad(Xijk) = Xj * (Vi – Vk)^R90 / 2A + Xk * (Vj – Vi)^R90 / 2A +
% -Xi * (Vi – Vk)^R90 / 2A – Xi * (Vj – Vi)^R90 / 2A
% where Xi is the scalar value at vertex i, Vi is the 3D position of vertex
% i, and A is the area of triangle (i,j,k). ^R90 represent a rotation of
% 90 degrees
%
% renaming indices of vertices of triangles for convenience
i1 = F(:,1); i2 = F(:,2); i3 = F(:,3);
% #F x 3 matrices of triangle edge vectors, named after opposite vertices
v32 = V(i3,:) – V(i2,:); v13 = V(i1,:) – V(i3,:); v21 = V(i2,:) – V(i1,:);
% area of parallelogram is twice area of triangle
% area of parallelogram is || v1 x v2 ||
n = cross(v32,v13,2);
% This does correct l2 norm of rows, so that it contains #F list of twice
% triangle areas
dblA = normrow(n);
% now normalize normals to get unit normals
u = normalizerow(n);
% rotate each vector 90 degrees around normal
%eperp21 = bsxfun(@times,normalizerow(cross(u,v21)),normrow(v21)./dblA);
%eperp13 = bsxfun(@times,normalizerow(cross(u,v13)),normrow(v13)./dblA);
eperp21 = bsxfun(@times,cross(u,v21),1./dblA);
eperp13 = bsxfun(@times,cross(u,v13),1./dblA);
%g = …
% ( …
% repmat(X(F(:,2)) – X(F(:,1)),1,3).*eperp13 + …
% repmat(X(F(:,3)) – X(F(:,1)),1,3).*eperp21 …
% );
GI = …
[0*size(F,1)+repmat(1:size(F,1),1,4) …
1*size(F,1)+repmat(1:size(F,1),1,4) …
2*size(F,1)+repmat(1:size(F,1),1,4)]’;
GJ = repmat([F(:,2);F(:,1);F(:,3);F(:,1)],3,1);
GV = [eperp13(:,1);-eperp13(:,1);eperp21(:,1);-eperp21(:,1); …
eperp13(:,2);-eperp13(:,2);eperp21(:,2);-eperp21(:,2); …
eperp13(:,3);-eperp13(:,3);eperp21(:,3);-eperp21(:,3)];
G = sparse(GI,GJ,GV, 3*size(F,1), size(V,1));
%% Alternatively
%%
%% f(x) is piecewise-linear function:
%%
%% f(x) = ∑ φi(x) fi, f(x ∈ T) = φi(x) fi + φj(x) fj + φk(x) fk
%% ∇f(x) = … = ∇φi(x) fi + ∇φj(x) fj + ∇φk(x) fk
%% = ∇φi fi + ∇φj fj + ∇φk) fk
%%
%% ∇φi = 1/hjk ((Vj-Vk)/||Vj-Vk||)^perp =
%% = ||Vj-Vk|| /(2 Aijk) * ((Vj-Vk)/||Vj-Vk||)^perp
%% = 1/(2 Aijk) * (Vj-Vk)^perp
%%
%m = size(F,1);
%eperp32 = bsxfun(@times,cross(u,v32),1./dblA);
%G = sparse( …
% [0*m + repmat(1:m,1,3) …
% 1*m + repmat(1:m,1,3) …
% 2*m + repmat(1:m,1,3)]’, …
% repmat([F(:,1);F(:,2);F(:,3)],3,1), …
% [eperp32(:,1);eperp13(:,1);eperp21(:,1); …
% eperp32(:,2);eperp13(:,2);eperp21(:,2); …
% eperp32(:,3);eperp13(:,3);eperp21(:,3)], …
% 3*m,size(V,1));
if dim == 2
G = G(1:(size(F,1)*dim),:);
end
% Should be the same as:
% g = …
% bsxfun(@times,X(F(:,1)),cross(u,v32)) + …
% bsxfun(@times,X(F(:,2)),cross(u,v13)) + …
% bsxfun(@times,X(F(:,3)),cross(u,v21));
% g = bsxfun(@rdivide,g,dblA);
case 4
% really dealing with tets
T = F;
% number of dimensions
assert(dim == 3);
% number of vertices
n = size(V,1);
% number of elements
m = size(T,1);
% simplex size
assert(size(T,2) == 4);
if m == 1 && ~isnumeric(V)
simple_volume = @(ad,r) -sum(ad.*r,2)./6;
simple_volume = @(ad,bd,cd) simple_volume(ad, …
[bd(:,2).*cd(:,3)-bd(:,3).*cd(:,2), …
bd(:,3).*cd(:,1)-bd(:,1).*cd(:,3), …
bd(:,1).*cd(:,2)-bd(:,2).*cd(:,1)]);
P = sym(‘P’,[1 3]);
V1 = V(T(:,1),:);
V2 = V(T(:,2),:);
V3 = V(T(:,3),:);
V4 = V(T(:,4),:);
V1P = V1-P;
V2P = V2-P;
V3P = V3-P;
V4P = V4-P;
A1 = simple_volume(V2P,V4P,V3P);
A2 = simple_volume(V1P,V3P,V4P);
A3 = simple_volume(V1P,V4P,V2P);
A4 = simple_volume(V1P,V2P,V3P);
B = [A1 A2 A3 A4]./(A1+A2+A3+A4);
G = simplify(jacobian(B,P).’);
return
end
% f(x) is piecewise-linear function:
%
% f(x) = ∑ φi(x) fi, f(x ∈ T) = φi(x) fi + φj(x) fj + φk(x) fk + φl(x) fl
% ∇f(x) = … = ∇φi(x) fi + ∇φj(x) fj + ∇φk(x) fk + ∇φl(x) fl
% = ∇φi fi + ∇φj fj + ∇φk fk + ∇φl fl
%
% ∇φi = 1/hjk = Ajkl / 3V * (Facejkl)^perp
% = Ajkl / 3V * (Vj-Vk)x(Vl-Vk)
% = Ajkl / 3V * Njkl / ||Njkl||
%
% get all faces
F = [ …
T(:,1) T(:,2) T(:,3); …
T(:,1) T(:,3) T(:,4); …
T(:,1) T(:,4) T(:,2); …
T(:,2) T(:,4) T(:,3)];
% compute areas of each face
A = doublearea(V,F)/2;
N = normalizerow(normals(V,F));
% compute volume of each tet
vol = volume(V,T);
GI = …
[0*m + repmat(1:m,1,4) …
1*m + repmat(1:m,1,4) …
2*m + repmat(1:m,1,4)];
GJ = repmat([T(:,4);T(:,2);T(:,3);T(:,1)],3,1);
GV = repmat(A./(3*repmat(vol,4,1)),3,1).*N(:);
G = sparse(GI,GJ,GV, 3*m,n);
end
end
the sphere.off file (mesh data) I am trying to compute this quantity here H grad(H) where H is the mean curvture
The code breaks at the last line because there is a dimension mismatch between H and grad(H). size(grad(H))= [960 3] and size(V)=[162 3] and size(mean_curvture)=[162 3]. Why is the grad(H) has that dimension? where is the mistake?
% Load a mesh
[V, F] = load_mesh(‘sphere.off’);
% Compute Mean Curvature and Normal Vector
% 2. Compute Cotangent Laplacian (L) and Mass Matrix (M)
L = cotmatrix(V, F);
M = massmatrix(V, F, ‘barycentric’);
% 3. Compute Mean Curvature Vector (H)
H = -inv(M) * (L * V);
% 4. Compute the Magnitude of Mean Curvature (mean curvature at each vertex)
mean_curvature = sqrt(sum(H.^2, 2));
% Compute the Normal Vector
normals = per_vertex_normals(V, F);
% Compute Gaussian Curvature
gaussian_curvature = discrete_gaussian_curvature(V,F);
% Compute Laplacian of Mean Curvature
laplacian_H = L * mean_curvature;
% Compute the Gradient of Mean Curvature
% Use gptoolbox’s gradient operator for scalar fields
grad_H = grad(V, F) * mean_curvature;
% Compute the new term: newthing
%H_squared = mean_curvature .^ 2;
%newthing = laplacian_H + ((H_squared – 2 * gaussian_curvature) .* mean_curvature) .* normals + mean_curvature .* grad_H;
% Compute the new term: newthing
H_squared = mean_curvature .^ 2;
newthing = laplacian_H + ((H_squared – 2 * gaussian_curvature) .* mean_curvature) .* normals + mean_curvature .* grad_H;
% Define small arc equation
h = 0.1; % Small parameter h
vertex = V; % V contains the vertices
% Compute the small arc
f_h = vertex + (mean_curvature .* normals) * h + newthing * (h^2 / 2);
% Visualization of small arc
figure;
trisurf(F, V(:,1), V(:,2), V(:,3), ‘FaceColor’, ‘cyan’, ‘EdgeColor’, ‘none’);
hold on;
plot3(f_h(:,1), f_h(:,2), f_h(:,3), ‘r.’, ‘MarkerSize’, 10);
axis equal;
title(‘Small Arc Visualization’);
xlabel(‘X’);
ylabel(‘Y’);
zlabel(‘Z’);
the grad(H) script
% Step 1: Load a Mesh
[V, F] = load_mesh(‘sphere.off’);
% Step 2: Compute Laplace-Beltrami Operator
L = cotmatrix(V, F); % cotangent Laplace-Beltrami operator
% Step 3: Compute the Mean Curvature Vector (H)
H = -L * V; % Mean curvature normal vector
% Step 4: Compute the Mean Curvature Magnitude
mean_curvature = sqrt(sum(H.^2, 2));
% Step 5: Compute the Gradient of Mean Curvature
% Use gptoolbox’s gradient operator for scalar fields
grad_H = grad(V, F) * mean_curvature;
% Step 6: Visualization or Further Processing
figure;
trisurf(F, V(:,1), V(:,2), V(:,3), mean_curvature, ‘EdgeColor’, ‘none’);
axis equal;
lighting gouraud;
camlight;
colorbar;
title(‘Mean Curvature Gradient’);
xlabel(‘X’);
ylabel(‘Y’);
zlabel(‘Z’);
grad(V, F) from gp toolbox
function [G] = grad(V,F)
% GRAD Compute the numerical gradient operator for triangle or tet meshes
%
% G = grad(V,F)
%
% Inputs:
% V #vertices by dim list of mesh vertex positions
% F #faces by simplex-size list of mesh face indices
% Outputs:
% G #faces*dim by #V Gradient operator
%
% Example:
% L = cotmatrix(V,F);
% G = grad(V,F);
% dblA = doublearea(V,F);
% GMG = -G’*repdiag(diag(sparse(dblA)/2),size(V,2))*G;
%
% % Columns of W are scalar fields
% G = grad(V,F);
% % Compute gradient magnitude for each column in W
% GM = squeeze(sqrt(sum(reshape(G*W,size(F,1),size(V,2),size(W,2)).^2,2)));
%
dim = size(V,2);
ss = size(F,2);
switch ss
case 2
% Edge lengths
len = normrow(V(F(:,2),:)-V(F(:,1),:));
% Gradient is just staggered grid finite difference
G = sparse(repmat(1:size(F,1),2,1)’,F,[1 -1]./len, size(F,1),size(V,1));
case 3
% append with 0s for convenience
if size(V,2) == 2
V = [V zeros(size(V,1),1)];
end
% Gradient of a scalar function defined on piecewise linear elements (mesh)
% is constant on each triangle i,j,k:
% grad(Xijk) = (Xj-Xi) * (Vi – Vk)^R90 / 2A + (Xk-Xi) * (Vj – Vi)^R90 / 2A
% grad(Xijk) = Xj * (Vi – Vk)^R90 / 2A + Xk * (Vj – Vi)^R90 / 2A +
% -Xi * (Vi – Vk)^R90 / 2A – Xi * (Vj – Vi)^R90 / 2A
% where Xi is the scalar value at vertex i, Vi is the 3D position of vertex
% i, and A is the area of triangle (i,j,k). ^R90 represent a rotation of
% 90 degrees
%
% renaming indices of vertices of triangles for convenience
i1 = F(:,1); i2 = F(:,2); i3 = F(:,3);
% #F x 3 matrices of triangle edge vectors, named after opposite vertices
v32 = V(i3,:) – V(i2,:); v13 = V(i1,:) – V(i3,:); v21 = V(i2,:) – V(i1,:);
% area of parallelogram is twice area of triangle
% area of parallelogram is || v1 x v2 ||
n = cross(v32,v13,2);
% This does correct l2 norm of rows, so that it contains #F list of twice
% triangle areas
dblA = normrow(n);
% now normalize normals to get unit normals
u = normalizerow(n);
% rotate each vector 90 degrees around normal
%eperp21 = bsxfun(@times,normalizerow(cross(u,v21)),normrow(v21)./dblA);
%eperp13 = bsxfun(@times,normalizerow(cross(u,v13)),normrow(v13)./dblA);
eperp21 = bsxfun(@times,cross(u,v21),1./dblA);
eperp13 = bsxfun(@times,cross(u,v13),1./dblA);
%g = …
% ( …
% repmat(X(F(:,2)) – X(F(:,1)),1,3).*eperp13 + …
% repmat(X(F(:,3)) – X(F(:,1)),1,3).*eperp21 …
% );
GI = …
[0*size(F,1)+repmat(1:size(F,1),1,4) …
1*size(F,1)+repmat(1:size(F,1),1,4) …
2*size(F,1)+repmat(1:size(F,1),1,4)]’;
GJ = repmat([F(:,2);F(:,1);F(:,3);F(:,1)],3,1);
GV = [eperp13(:,1);-eperp13(:,1);eperp21(:,1);-eperp21(:,1); …
eperp13(:,2);-eperp13(:,2);eperp21(:,2);-eperp21(:,2); …
eperp13(:,3);-eperp13(:,3);eperp21(:,3);-eperp21(:,3)];
G = sparse(GI,GJ,GV, 3*size(F,1), size(V,1));
%% Alternatively
%%
%% f(x) is piecewise-linear function:
%%
%% f(x) = ∑ φi(x) fi, f(x ∈ T) = φi(x) fi + φj(x) fj + φk(x) fk
%% ∇f(x) = … = ∇φi(x) fi + ∇φj(x) fj + ∇φk(x) fk
%% = ∇φi fi + ∇φj fj + ∇φk) fk
%%
%% ∇φi = 1/hjk ((Vj-Vk)/||Vj-Vk||)^perp =
%% = ||Vj-Vk|| /(2 Aijk) * ((Vj-Vk)/||Vj-Vk||)^perp
%% = 1/(2 Aijk) * (Vj-Vk)^perp
%%
%m = size(F,1);
%eperp32 = bsxfun(@times,cross(u,v32),1./dblA);
%G = sparse( …
% [0*m + repmat(1:m,1,3) …
% 1*m + repmat(1:m,1,3) …
% 2*m + repmat(1:m,1,3)]’, …
% repmat([F(:,1);F(:,2);F(:,3)],3,1), …
% [eperp32(:,1);eperp13(:,1);eperp21(:,1); …
% eperp32(:,2);eperp13(:,2);eperp21(:,2); …
% eperp32(:,3);eperp13(:,3);eperp21(:,3)], …
% 3*m,size(V,1));
if dim == 2
G = G(1:(size(F,1)*dim),:);
end
% Should be the same as:
% g = …
% bsxfun(@times,X(F(:,1)),cross(u,v32)) + …
% bsxfun(@times,X(F(:,2)),cross(u,v13)) + …
% bsxfun(@times,X(F(:,3)),cross(u,v21));
% g = bsxfun(@rdivide,g,dblA);
case 4
% really dealing with tets
T = F;
% number of dimensions
assert(dim == 3);
% number of vertices
n = size(V,1);
% number of elements
m = size(T,1);
% simplex size
assert(size(T,2) == 4);
if m == 1 && ~isnumeric(V)
simple_volume = @(ad,r) -sum(ad.*r,2)./6;
simple_volume = @(ad,bd,cd) simple_volume(ad, …
[bd(:,2).*cd(:,3)-bd(:,3).*cd(:,2), …
bd(:,3).*cd(:,1)-bd(:,1).*cd(:,3), …
bd(:,1).*cd(:,2)-bd(:,2).*cd(:,1)]);
P = sym(‘P’,[1 3]);
V1 = V(T(:,1),:);
V2 = V(T(:,2),:);
V3 = V(T(:,3),:);
V4 = V(T(:,4),:);
V1P = V1-P;
V2P = V2-P;
V3P = V3-P;
V4P = V4-P;
A1 = simple_volume(V2P,V4P,V3P);
A2 = simple_volume(V1P,V3P,V4P);
A3 = simple_volume(V1P,V4P,V2P);
A4 = simple_volume(V1P,V2P,V3P);
B = [A1 A2 A3 A4]./(A1+A2+A3+A4);
G = simplify(jacobian(B,P).’);
return
end
% f(x) is piecewise-linear function:
%
% f(x) = ∑ φi(x) fi, f(x ∈ T) = φi(x) fi + φj(x) fj + φk(x) fk + φl(x) fl
% ∇f(x) = … = ∇φi(x) fi + ∇φj(x) fj + ∇φk(x) fk + ∇φl(x) fl
% = ∇φi fi + ∇φj fj + ∇φk fk + ∇φl fl
%
% ∇φi = 1/hjk = Ajkl / 3V * (Facejkl)^perp
% = Ajkl / 3V * (Vj-Vk)x(Vl-Vk)
% = Ajkl / 3V * Njkl / ||Njkl||
%
% get all faces
F = [ …
T(:,1) T(:,2) T(:,3); …
T(:,1) T(:,3) T(:,4); …
T(:,1) T(:,4) T(:,2); …
T(:,2) T(:,4) T(:,3)];
% compute areas of each face
A = doublearea(V,F)/2;
N = normalizerow(normals(V,F));
% compute volume of each tet
vol = volume(V,T);
GI = …
[0*m + repmat(1:m,1,4) …
1*m + repmat(1:m,1,4) …
2*m + repmat(1:m,1,4)];
GJ = repmat([T(:,4);T(:,2);T(:,3);T(:,1)],3,1);
GV = repmat(A./(3*repmat(vol,4,1)),3,1).*N(:);
G = sparse(GI,GJ,GV, 3*m,n);
end
end
the sphere.off file (mesh data) image processing, image segmentation, digital image processing, geometry processing, toolbox, matrix MATLAB Answers — New Questions
How to compute the Shapley value of BP neural network
Hello, I trained a BP nerual network using newff function, and wanted to obtain its Shapley value. But error occurs like this:
How can I deal with it?Hello, I trained a BP nerual network using newff function, and wanted to obtain its Shapley value. But error occurs like this:
How can I deal with it? Hello, I trained a BP nerual network using newff function, and wanted to obtain its Shapley value. But error occurs like this:
How can I deal with it? interpret machine learning models MATLAB Answers — New Questions
Find what toolboxes a script uses
Hello to all,
We need to buy a Matlab license for a project and we already have a working collection of scripts. However, there are so many functions used that it is very difficult to find out what toolboxes we need to buy (it will take a lot of time to check each function’s origin). Is there a tool created by Mathworks which will allow us to find out what toolboxes or Matlab packages are used for a given script?
Best regards,
JeanHello to all,
We need to buy a Matlab license for a project and we already have a working collection of scripts. However, there are so many functions used that it is very difficult to find out what toolboxes we need to buy (it will take a lot of time to check each function’s origin). Is there a tool created by Mathworks which will allow us to find out what toolboxes or Matlab packages are used for a given script?
Best regards,
Jean Hello to all,
We need to buy a Matlab license for a project and we already have a working collection of scripts. However, there are so many functions used that it is very difficult to find out what toolboxes we need to buy (it will take a lot of time to check each function’s origin). Is there a tool created by Mathworks which will allow us to find out what toolboxes or Matlab packages are used for a given script?
Best regards,
Jean toolbox, package, matlab, script, functions MATLAB Answers — New Questions