Can I use a struct in an anonymous function?
Show older comments
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
Accepted Answer
More Answers (1)
Mat
on 4 Oct 2024
0 votes
3 Comments
Bruno Luong
on 7 Oct 2024
"Solution has been posted."
Where?
Bruno Luong
on 10 Oct 2024
This is question not a solution, do I miss anything?
@Steven Lord answer should be accepted.
Categories
Find more on Programming in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!