another function for solving differential equation other than dsolve in MATLAB

Question

Leila Dehdari 2020-10-5

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/604957-another-function-for-solving-differential-equation-other-than-dsolve-in-matlab

评论： Leila Dehdari 2020-10-12

I am using dsolve function in MATLAB to solve a differential equation. It is so slow and can not end up with a result. Is there a way to speed it up, or is there another way for solving a differential equation? I don't want to solve it numerically.

The equations are as follow:

all the derivatives in above equations are replaced with equivalent expressions, so all are only function of y and T only.

I need to solve this equation:

Cg, Cs and qi are constants. The initial conditions for positive root is (y=0.01, T=0) and for negative root is (y=0.99, T=0).

Here is my code:

R=1.9872;                          %cal/mol/K
eps=0.763;
rhos=0.484;                        %gr/cm3
T0=300;                            %theta K
y0=0.01;
ys=0.99;
p=1;                              %atm
Upg=0.68e-3;                      %mol/cm2/s
L=100;                            %cm   column length
Cg=1.04*239.006/1000*28.0134;      % kj/kg/k  to  cal/mol/k   N2
Cs=0.22;                          %cal/gr/k
%  m(mol/gr)     b(1/atm)        q(cal/mol)
cA=[3.65e-3      2.8e-4            4900];  %  co2 BPL AC
cB=[3.65e-3      13.5e-4           2500];  %  N2 BPL AC 
i=1;
IPA=cA(i,:); IPB=cB(i,:);
IPA1=IPA(1);IPA2=IPA(2);IPA3=IPA(3);
IPB1=IPB(1);IPB2=IPB(2);IPB3=IPB(3);
syms T(y)
DN1T= @(y,T)    (IPA1*IPA2*p*y*exp(IPA3/(R*(T + T0)))*((IPA2*IPA3*p*y*exp(IPA3/(R*(T + T0))))/(R*(T + T0)^2) - (IPB2*IPB3*p*exp(IPB3/(R*(T + T0)))...
    *(y - 1))/(R*(T + T0)^2)))/(IPA2*p*y*exp(IPA3/(R*(T + T0))) - IPB2*p*exp(IPB3/(R*(T + T0)))*(y - 1) + 1)^2 -...
    (IPA1*IPA2*IPA3*p*y*exp(IPA3/(R*(T + T0))))/(R*(T + T0)^2*(IPA2*p*y*exp(IPA3/(R*(T + T0))) - ...
    IPB2*p*exp(IPB3/(R*(T + T0)))*(y - 1) + 1));
DN2T =@(y,T)     (IPB1*IPB2*IPB3*p*exp(IPB3/(R*(T + T0)))*(y - 1))/(R*(T + T0)^2*(IPA2*p*y*exp(IPA3/(R*(T + T0))) - ...
    IPB2*p*exp(IPB3/(R*(T + T0)))*(y - 1) + 1)) - (IPB1*IPB2*p*exp(IPB3/(R*(T + T0)))*(y - 1)*((IPA2*IPA3*p*y*...
    exp(IPA3/(R*(T + T0))))/(R*(T + T0)^2) - (IPB2*IPB3*p*exp(IPB3/(R*(T + T0)))*(y - 1))/(R*(T + T0)^2)))...
    /(IPA2*p*y*exp(IPA3/(R*(T + T0))) - IPB2*p*exp(IPB3/(R*(T + T0)))*(y - 1) + 1)^2;
 
 
DN1Y =@(y,T)      (IPA1*IPA2*p*exp(IPA3/(R*(T + T0))))/(IPA2*p*y*exp(IPA3/(R*(T + T0))) - IPB2*p*exp(IPB3/(R*(T + T0)))*(y - 1) + 1)...
    - (IPA1*IPA2*p*y*exp(IPA3/(R*(T + T0)))*(IPA2*p*exp(IPA3/(R*(T + T0))) - IPB2*p*exp(IPB3/(R*(T + T0)))))/...
    (IPA2*p*y*exp(IPA3/(R*(T + T0))) - IPB2*p*exp(IPB3/(R*(T + T0)))*(y - 1) + 1)^2;
 
DN2Y =@(y,T)      (IPB1*IPB2*p*exp(IPB3/(R*(T + T0)))*(y - 1)*(IPA2*p*exp(IPA3/(R*(T + T0))) - IPB2*p*exp(IPB3/(R*(T + T0)))))/...
    (IPA2*p*y*exp(IPA3/(R*(T + T0))) - IPB2*p*exp(IPB3/(R*(T + T0)))*(y - 1) + 1)^2 - (IPB1*IPB2*p*exp(IPB3/(R*(T + T0))))...
    /(IPA2*p*y*exp(IPA3/(R*(T + T0))) - IPB2*p*exp(IPB3/(R*(T + T0)))*(y - 1) + 1);
  
A=@(y,T)   Cg*(DN1T(y,T)-y*(DN1T(y,T)+DN2T(y,T)));
B=@(y,T)   (Cg*(T)*(DN1T(y,T)+DN2T(y,T))-Cs)+Cg*(DN1Y(y,T)-y*(DN1Y(y,T)+DN2Y(y,T)))+cA(i,3)*DN1T(y,T)+cB(i,3)*DN2T(y,T);
C=@(y,T)    Cg*(T)*(DN1Y(y,T)+DN2Y(y,T))+cA(i,3)*DN1Y(y,T)+cB(i,3)*DN2Y(y,T);
 
 S4=@(y,T)   (-B(y,T) +(B(y,T) ^2-4*A(y,T) *C(y,T) )^0.5)/2/A(y,T) ;
 S2=@(y,T)   (-B(y,T) -(B(y,T) ^2-4*A(y,T) *C(y,T) )^0.5)/2/A(y,T) ;
S4sol=dsolve([diff(T,y)==S4(y,T) , T(y0)==0])
S2sol=dsolve([diff(T,y)==S2(y,T) , T(ys)==0])