Can I use a struct in an anonymous function?

8 vues (au cours des 30 derniers jours)
Mat
Mat le 2 Oct 2024
Commenté : Bruno Luong le 10 Oct 2024
I am solving a PDE via finite volume methods. As a result, I am left with a set of ODEs which I solve with ode15s. I have a bunch of constants which I use in the function defining the derivative, and I insert these by the inclusion of a struct to the function. When I run the ode15s I obtain an error which pertains to the inclusion of the struct as an input to the function. I don't really want to include the constants within the function as that seems "messy", so I would prefer to do it via a struct.
%This is a code to solve the isothermal sintering equations.
%%---Physical parameters---
N = 251; %This is the number of cells-1
L = 1; %This is the initial length of the rod
%%---Initial conditions
%Length
X=linspace(0,L,N); %This is the partition of the initial length
%Density
rho_0_int=0.5+0.1*sin(6*pi*X);
rho_0=0.5*(rho_0_int(1:N-1)+rho_0_int(2:N));
%Amount of mass
h=cumtrapz(X,rho_0_int);
%Initial velocity
u0=zeros(N-1,1);
% Define the system of equations as a function
y0 = [1./rho_0'; u0];
% Finite Volume Method solution loop
tspan = [0 3]; % You can adjust this based on the time range you're interested in
[t, y] = ode15s(@ode_system, tspan, y0);
nu=y(:,1:N-1);
u=y(:,N:2*N-2);
Tspan=length(y(:,1));
for i=1:Tspan
plot(u(i,:))
pause(0.05);
end
function dydt = ode_system(t, y)
g = 9.81; %Acceleration due to gravity.
K = 1e-2; %This is the Laplace constant from the sintering potential.
eta_0 = 0.005; %The shear stress of the skeleton
T = 2.7; %This is the time of sintering.
rho_m = 1; %This is the density of the individual metallic particles.
N = 251; %This is the number of cells-1
L = 1; %Initial length of rod
X=linspace(0,L,N); %This is the partition of the initial length
%Density
rho_0_int=0.5+0.1*sin(6*pi*X);
rho_0=0.5*(rho_0_int(1:N-1)+rho_0_int(2:N));
%Amount of mass
h=cumtrapz(X,rho_0_int);
M=h(N);
U=M/(rho_m*T);
a_1=T/(U*M);
a_2=T*rho_m/M^2;
a_3=g*T/U;
dh=diff(h);
nu = y(1:N-1)'; % u(t)
u = y(N:2*N-2)'; % v(t)
%Compute the left and right specific volumes:
nu_left=15*nu(1)/8-5*nu(2)/4+3*nu(3)/8;
nu_right=15*nu(N-3)/8-5*nu(N-2)/5+3*nu(N-1)/8;
%Compute the fluxes on the interfaces
nu_half=NaN(1,N);
nu_half(2:N-1)=0.5*(nu(2:N-1)+nu(1:N-2));
nu_half(1)=nu_left;
nu_half(N)=nu_right;
%Fluxes for nu equation
nu_flux=zeros(N,1);
nu_flux(2:N-1)=0.5*(u(1:N-2)+u(2:N-1));
du_dh_left=-P_L(K,nu_left)*nu_left/zeta(nu_left);
du_dh_right=-P_L(K,nu_right)*nu_right/zeta(nu_right);
nu_flux(1)=-3*du_dh_left*dh(1)/8+9*u(1)/8-u(2)/8;
nu_flux(N)=3*du_dh_right*dh(N-1)/8-9*u(N-1)/8+u(N-2)/8;
%Fluxes for u equation
u_grad=2*(u(2:N-1)-u(1:N-2))./(dh(2:N-1)+dh(1:N-2));
u_flux=zeros(1,N);
u_flux(2:N-1)=a_1*P_L(K,nu_half(2:N-1))+(a_2*zeta(nu_half(2:N-1))./nu_half(2:N-1)).*u_grad;
% The system of equationsY
dnu_dt = nu_flux(2:N)-nu_flux(1:N-1);
du_dt =u_flux(2:N)'-u_flux(1:N-1)';
% Return the derivatives
dydt = [dnu_dt; du_dt];
end
function y=P_L(K,nu)
y=K./nu;
end
function y=zeta(nu)
y=(2/3)*(nu).^(-2)./(nu-1);
end

Réponse acceptée

Steven Lord
Steven Lord le 2 Oct 2024
When I run the ode15s I obtain an error which pertains to the inclusion of the struct as an input to the function.
Please show us the full and exact text of that error message, including all the text displayed in red in the Command Window (and if there are any warning messages displayed in orange, please show us those too.) The exact text may be useful and/or necessary to determine what's going on and how to avoid the warning and/or error.
At a quick glance, though, you're not actually including the struct as an input to your ODE function.
[t, y] = ode15s(@ode_system, tspan, y0);
function dydt = ode_system(t, y, S)
As you're calling it, ode15s will call your ode_system function with two input arguments, t and y. You have not told it that it should pass the struct S into that function; use one of the techniques shown on this documentation page to do so. [I suspect you expected that MATLAB would automatically "figure out" it should pass the struct from the caller's workspace into ode_system based on the name. It does not.]
  17 commentaires
Bruno Luong
Bruno Luong le 7 Oct 2024
Modifié(e) : Bruno Luong le 7 Oct 2024
Put the parameter variables in a single extra structure, then dispatche it using input structure within your inner function. This is less messy and more efficient than calling the anonymous function with a long list of parameters.
Advantage: When you want to add or remove parameters you only need to change the structure, not the calling interface.
You could do similar with a specific own object class and put parameters in properties when you instantiate the object, but if you are not comfortable with structure no need to work with OOP.
Stephen23
Stephen23 le 7 Oct 2024
Modifié(e) : Stephen23 le 8 Oct 2024
"Turns out there was."
Turns out there is not.
"The issue is ode15s, where it uses 2 inputs rather than three as I origianally thought...."
As everyone here keeps advising you, that is exactly what function parameterization is for:
Show us your current code and we will show you how.

Connectez-vous pour commenter.

Plus de réponses (1)

Mat
Mat le 4 Oct 2024
Solution has been posted.
  3 commentaires
Bruno Luong
Bruno Luong le 10 Oct 2024
This is question not a solution, do I miss anything?
@Steven Lord answer should be accepted.

Connectez-vous pour commenter.

Produits


Version

R2024a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by