Can I use a struct in an anonymous function?
이전 댓글 표시
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
채택된 답변
추가 답변 (1개)
Mat
2024년 10월 4일
0 개 추천
댓글 수: 3
Bruno Luong
2024년 10월 7일
"Solution has been posted."
Where?
Bruno Luong
2024년 10월 10일
This is question not a solution, do I miss anything?
@Steven Lord answer should be accepted.
카테고리
도움말 센터 및 File Exchange에서 Programming에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!