# another function for solving differential equation other than dsolve in MATLAB

3 views (last 30 days)
Leila Dehdari on 5 Oct 2020
Commented: Leila Dehdari on 12 Oct 2020
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])

Leila Dehdari on 5 Oct 2020
Hi Bjorn, thanks for the response. Can you give me an example of your suggestion? if you mean diff(T,y) after dsolve function, this is the way that MATLAB help document shows.
Bjorn Gustavsson on 5 Oct 2020
No, I mean write the equations like you would see them in a book or article:
That makes it far easier for others to read than your long lines of code.
Leila Dehdari on 5 Oct 2020
Thanks Bjorn, I added it to the question

Kiran Felix Robert on 8 Oct 2020
Hi Leila,
One way to speed-up the execution is to add a limit on the maximum degree of radicals. The following shows you an example,
S = dsolve(..,'MaxDegree',2);
This will force dsolve to assume implicit formulas for polynomials of degree greater than the specified value.
Kiran Felix Robert

#### 1 Comment

Leila Dehdari on 12 Oct 2020
Thanks Kiran