I need help to solve differential equations??

First M file:

us=1;
Mm=29;
pt=1.0e5;
rhob=1300;
pb0=0.211e5;
Cp=0.992;
H=1285409;
U=0.096;
dt=2.54e-2;
Tr=625.0;
R=8314.5;
z0=0;
zend=3;
m=100;
t0=0;
tend=10;
n=9;
pin=0.015e5;
Tin=625;
par = [us Mm pt rhob pb0 Cp H U dt Tr R m pin Tin];
p0=zeros(1,m);
p0(1)=pin;
T0=625*ones(1,m);
z=[z0:(zend-z0)/(m-1):zend];
t= [t0:(tend-t0)/(n-1):tend];
init=[p0 T0];
M=eye(2*m,2*m);
M(1,1)=0;
M(m+1, m+1)=0;
M=spares(M);
options=odeset('Mass',M,'Stats', 'on');
[t,x]=ode15s(@derivt,t,init,options,par,z);
p_out=x(:,1:m);
T_out=x(:,(m+1):(m*2));
l=2;
n=1;
nr=1;
subplot(l,n,nr);
plot(z.p_out(2,:));
title('Partial pressure profile');
xlabel('z [m]');
ylabel('p [Pa]');
nr=nr+1;
subplot(l,n,nr);
plot(z.T_out(2,:));
title('Temperature profile');
xlabel('z [m]');
ylabel('T [K]');

second M file:

function derivt_alt=deriv(t,var,par,z)
us=par(1);
Mm=par(2);
pt=par(3);
rhob=par(4);
pb0=par(5);
Cp=par(6);
H=par(7);
U=par(8);
dt=par(9);
Tr=par(10);
R=par(11);
m=par(12);
pin=par(13);
Tin=par(14);
p=var(1:m);
T=var((m+1):2*m);
dpdz=dss020(z(1).z(m),m,p,1);
dTdz=dss020(z(1),z(m),m,T,1);
p_i=p(2:m);
T_i=T(2:m);
rhog=(Mm*pt)./(R*T_i);
k=1.0e-10*exp(19.837-1363.0./T_i)/3600;
r=k.*p_i*pb0;
dTdt=1.0./(rhog*cp).*(H*rhob*r-(4.0*U/dt)*(T_i-Tr))-us*dTdz(2:m);
dpdt=p_i./T_i.*dTdt-us*dpdz(2:m)+us*p_i./T_i.*dTdz(2:m)-(Mm*pt*rhob./rhog.*r);
resp=pin-p(1);
resT=Tin-T(1);
derivt_alt=[resp,dpdt,resT,dTdt]';

[Star Strider]

Please use the [{}Code] button to format your code so it is readable.

Answer 1

[Star Strider]

I didn’t run your code, but there is one obvious problem and one somewhat less than obvious problem:

[t,x]=ode15s(@(t,var)deriv(t,var,par,z),t,init,options,par,z);

The obvious problem is that you spelled your ‘deriv’ function ‘derivt’ in your ode15s call. The less than obvious one is that the ODE solvers only want the first two arguments, so only pass those to ode15s. This syntax: @(t,var)deriv(t,var,par,z) does that. See the ODE solver documentation for details.

Change that line and see if your code gives the desired result.

[shahid hussain ansari]

I did but in vain,

us=1;
Mm=29;
pt=1.0e5;
rhob=1300;
pb0=0.211e5;
Cp=0.992;
H=1285409;
U=0.096;
dt=2.54e-2;
Tr=625.0;
R=8314.5;
z0=0;
zend=3;
m=100;
t0=0;
tend=10;
n=9;
pin=0.015e5;
Tin=625;
par = [us Mm pt rhob pb0 Cp H U dt Tr R m pin Tin];
p0=zeros(1,m);
p0(1)=pin;
T0=625*ones(1,m);
z=[z0:(zend-z0)/(m-1):zend];
t= [t0:(tend-t0)/(n-1):tend];
init=[p0 T0];
M=eye(2*m,2*m);
M(1,1)=0;
M(m+1, m+1)=0;
M=sparse(M);
options=odeset('Mass',M,'Stats', 'on');
[t,x]=ode15s(@deriv,t,init,options,par,z);
p_out=x(:,1:m);
T_out=x(:,(m+1):(m*2));
l=2;
n=1;
nr=1;
subplot(l,n,nr);
plot(z.p_out(2,:));
title('Partial pressure profile');
xlabel('z [m]');
ylabel('p [Pa]');
nr=nr+1;
subplot(l,n,nr);
plot(z.T_out(2,:));
title('Temperature profile');
xlabel('z [m]');
ylabel('T [K]');

[Star Strider]

First, please format your code.

Second, what was the problem? What did your code do that you did not want it to, or what did it not do that you wanted it to?