MATLAB Answers

another function for solving differential equation other than dsolve in MATLAB

2 views (last 30 days)
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])

  5 Comments

Show 2 older comments
Leila Dehdari
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
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.

Sign in to comment.

Accepted Answer

Kiran Felix Robert
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.
Refer to the dsolve documentation for more information.
Kiran Felix Robert

More Answers (0)

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!