Can I use a struct in an anonymous function?

조회 수: 8 (최근 30일)
Mat
Mat 2024년 10월 2일
댓글: Bruno Luong 2024년 10월 10일
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

채택된 답변

Steven Lord
Steven Lord 2024년 10월 2일
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
Bruno Luong
Bruno Luong 2024년 10월 7일
편집: Bruno Luong 2024년 10월 7일
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 2024년 10월 7일
편집: Stephen23 2024년 10월 8일
"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.

댓글을 달려면 로그인하십시오.

추가 답변 (1개)

Mat
Mat 2024년 10월 4일
Solution has been posted.
  댓글 수: 3
Bruno Luong
Bruno Luong 2024년 10월 10일
This is question not a solution, do I miss anything?
@Steven Lord answer should be accepted.

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기

제품


릴리스

R2024a

Community Treasure Hunt

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

Start Hunting!

Translated by