Group Abstract

Message Boards

WOLFRAM COMMUNITY

8.4K Views

8 Replies

0 Total Likes

View groups...

Follow this post

Share this post:

GROUPS:

Solve a linear system of ODEs?

Abed Alnaif

Posted 6 years ago

Hello, I'm trying to solve the following system of ODEs, but it is taking very long to solve. Does anyone know why it's taking so long to solve, and whether there's a way to speed up the computation? Please refer to the code below or the attached file. Thank you, Abed de = {Derivative[1][DC][t] == -kE DC[t] - kPT DC[t] + kTP DP[t] + kD MC[t], Derivative[1][DP][t] == kPT DC[t] - kTP DP[t], Derivative[1][IC][t] == -2 kD IC[t] - kE IC[t] - kPT IC[t] + kTP IP[t], Derivative[1][IP][t] == kPT IC[t] - kTP IP[t], Derivative[1][MC][t] == 2 kD IC[t] - kD MC[t] - kE MC[t] - kPT MC[t] + kTP MP[t], Derivative[1][MP][t] == kPT MC[t] - kTP MP[t]} ic = {IC[0] == IC0, IP[0] == 0, MC[0] == 0, MP[0] == 0, DC[0] == 0, DP[0] == 0}; sol = DSolve[Join[de, ic], {IC[t], IP[t], MC[t], MP[t], DC[t], DP[t]}, t]; Attachments: debug.nb

Hello,

I'm trying to solve the following system of ODEs, but it is taking very long to solve. Does anyone know why it's taking so long to solve, and whether there's a way to speed up the computation? Please refer to the code below or the attached file.

Thank you,

Abed

de = {Derivative[1][DC][t] == -kE DC[t] - kPT DC[t] + kTP DP[t] + 
    kD MC[t], Derivative[1][DP][t] == kPT DC[t] - kTP DP[t], 
  Derivative[1][IC][t] == -2 kD IC[t] - kE IC[t] - kPT IC[t] + 
    kTP IP[t], Derivative[1][IP][t] == kPT IC[t] - kTP IP[t], 
  Derivative[1][MC][t] == 
   2 kD IC[t] - kD MC[t] - kE MC[t] - kPT MC[t] + kTP MP[t], 
  Derivative[1][MP][t] == kPT MC[t] - kTP MP[t]}

ic = {IC[0] == IC0, IP[0] == 0, MC[0] == 0, MP[0] == 0, DC[0] == 0, 
   DP[0] == 0};

sol = DSolve[Join[de, ic], {IC[t], IP[t], MC[t], MP[t], DC[t], DP[t]},
    t];

POSTED BY: Abed Alnaif

8 Replies

Sort By:

Abed Alnaif

Posted 6 years ago

Hans, Thank you for sharing this method! I tried it out and it works. The reason for the zero solutions is because the initial condition was placed on the wrong state. Placing it on the correct state (the third state) results in non-zero solutions. Thank you! Abed

POSTED BY: Abed Alnaif

Hans Dolhaine

Hans Dolhaine, retired

Posted 6 years ago

You get the same solution by diagonalising the matrix of coefficients: eigv = Transpose[Eigenvectors[mat /. values]] // Chop; diag = Inverse[eigv].(mat /. values).eigv // Chop; % // MatrixForm Then matnew = Map[If[# != 0, Exp[# t], 0] &, diag, {2}]; % // MatrixForm and xxv2 = eigv.matnew.Inverse[eigv].{IC0, 0, 0, 0, 0, 0} // Chop; % // MatrixForm compared to xxv // Chop // MatrixForm And this works even with your general (non-numerical) system: eigv = Transpose[Eigenvectors[mat]] // FullSimplify; diag = Inverse[eigv].(mat).eigv // FullSimplify; % // MatrixForm matnew = Map[If[# =!= 0, Exp[# t], 0] &, diag, {2}]; % // MatrixForm xxv2 = eigv.matnew.Inverse[eigv].{IC0, 0, 0, 0, 0, 0} // FullSimplify; % // MatrixForm Test if xxv2 is a solution: D[ xxv2 , t] - mat. xxv2 // FullSimplify

You get the same solution by diagonalising the matrix of coefficients:

eigv = Transpose[Eigenvectors[mat /. values]] // Chop;
diag = Inverse[eigv].(mat /. values).eigv // Chop;
% // MatrixForm

Then

matnew = Map[If[# != 0, Exp[# t], 0] &, diag, {2}];
% // MatrixForm

and

xxv2 = eigv.matnew.Inverse[eigv].{IC0, 0, 0, 0, 0, 0} // Chop;
% // MatrixForm

compared to

xxv // Chop // MatrixForm

And this works even with your general (non-numerical) system:

eigv = Transpose[Eigenvectors[mat]] // FullSimplify;
diag = Inverse[eigv].(mat).eigv // FullSimplify;
% // MatrixForm
matnew = Map[If[# =!= 0, Exp[# t], 0] &, diag, {2}];
% // MatrixForm

xxv2 = eigv.matnew.Inverse[eigv].{IC0, 0, 0, 0, 0, 0} // FullSimplify;
% // MatrixForm

Test if xxv2 is a solution:

D[ xxv2 , t] - mat. xxv2 // FullSimplify

POSTED BY: Hans Dolhaine

Hans Dolhaine

Hans Dolhaine, retired

Posted 6 years ago

As you said you have a linear system of ODE's. You may write it in MatrixForm: x' == A x like this: your equations de = {Derivative[1][DC][t] == -kE DC[t] - kPT DC[t] + kTP DP[t] + kD MC[t], Derivative[1][DP][t] == kPT DC[t] - kTP DP[t], Derivative[1][IC][t] == -2 kD IC[t] - kE IC[t] - kPT IC[t] + kTP IP[t], Derivative[1][IP][t] == kPT IC[t] - kTP IP[t], Derivative[1][MC][t] == 2 kD IC[t] - kD MC[t] - kE MC[t] - kPT MC[t] + kTP MP[t], Derivative[1][MP][t] == kPT MC[t] - kTP MP[t]} Get rid of the left-hand sides de1 = (de /. a_ == b_ -> b); and make the matrix A (here called mat) mat = Table[ Coefficient[de1[[j]], #] & /@ {DC[t], DP[t], IC[t], IP[t], MC[t], MP[t]}, {j, 1, Length[de1]}]; % // MatrixForm And really we have a system equivalent to the first one sys = mat.{DC[t], DP[t], IC[t], IP[t], MC[t], MP[t]}; (de /. a_ == b_ -> b) - sys // Simplify Then your solution should be x = Exp[ A t ] . x0 or xx = MatrixExp[mat t].{IC0, 0, 0, 0, 0, 0} which gives a very large output which can't be handled conveniently. So I define some numerical values values = {kD -> .03, kE -> 1.4, kPT -> .75, kTP -> 3.1} and the numerical solution xxv = xx /. values which interestingly has only two components. But it fulfills your system of equations: (mat /. values) .xxv - D[ xxv, t] // FullSimplify // Chop

As you said you have a linear system of ODE's.

You may write it in MatrixForm:

x' == A  x

like this:

your equations

de = {Derivative[1][DC][t] == -kE DC[t] - kPT DC[t] + kTP DP[t] + 
    kD MC[t], Derivative[1][DP][t] == kPT DC[t] - kTP DP[t], 
  Derivative[1][IC][t] == -2 kD IC[t] - kE IC[t] - kPT IC[t] + 
    kTP IP[t], Derivative[1][IP][t] == kPT IC[t] - kTP IP[t], 
  Derivative[1][MC][t] == 
   2 kD IC[t] - kD MC[t] - kE MC[t] - kPT MC[t] + kTP MP[t], 
  Derivative[1][MP][t] == kPT MC[t] - kTP MP[t]}

Get rid of the left-hand sides

de1 = (de /. a_ == b_ -> b);

and make the matrix A (here called mat)

mat = Table[
   Coefficient[de1[[j]], #] & /@ {DC[t], DP[t], IC[t], IP[t], MC[t],  MP[t]}, {j, 1, Length[de1]}];
% // MatrixForm

And really we have a system equivalent to the first one

sys = mat.{DC[t], DP[t], IC[t], IP[t], MC[t], MP[t]};
(de /. a_ == b_ -> b) - sys // Simplify

Then your solution should be

x = Exp[ A t ] . x0

xx =   MatrixExp[mat t].{IC0, 0, 0, 0, 0, 0}

which gives a very large output which can't be handled conveniently.

So I define some numerical values

values = {kD -> .03, kE -> 1.4, kPT -> .75, kTP -> 3.1}

and the numerical solution

xxv = xx /. values

which interestingly has only two components.

But it fulfills your system of equations:

(mat /. values) .xxv -  D[ xxv, t] // FullSimplify // Chop

POSTED BY: Hans Dolhaine

Bill Nelson

Posted 6 years ago

I experimented a bit with forms similar to what you expect and didn't find anything where Solve could find a solution. Perhaps you can experiment by substituting various forms into the example I showed and see if you can find a solution. Or perhaps someone else can see something in this that we have missed.

POSTED BY: Bill Nelson

Abed Alnaif

Posted 6 years ago

Thanks again, Bill. Also, one more question that I'm confused about: I noticed that the solution contains exponentials of the form e^t. How can this be possible, physically, given that the unit of t in this problem is time? Then what are the units of e^t? I was expecting a solution of the form e^(k*t), where k represents the various parameters in the system with units 1/time.

POSTED BY: Abed Alnaif

Bill Nelson

Posted 6 years ago

POSTED BY: Bill Nelson

Abed Alnaif

Posted 6 years ago

Thank you, Bill. Out of curiosity, and for future reference, do we know why my original code didn't work?

POSTED BY: Abed Alnaif

Bill Nelson

Posted 6 years ago

Because the structure of the system is so simple, by inspection guess the form of one solution DC[t_]:=mDC E^-t+bDC; DP[t_]:=mDP E^-t+bDP; IC[t_]:=mIC E^-t+bIC; IP[t_]:=mIP E^-t+bIP; MC[t_]:=mMC E^-t+bMC; MP[t_]:=mMP E^-t+bMP; de = {Derivative[1][DC][t] == -kE DC[t] - kPT DC[t] + kTP DP[t] + kD MC[t], Derivative[1][DP][t] == kPT DC[t] - kTP DP[t], Derivative[1][IC][t] == -2 kD IC[t] - kE IC[t] - kPT IC[t] + kTP IP[t], Derivative[1][IP][t] == kPT IC[t] - kTP IP[t], Derivative[1][MC][t] == 2 kD IC[t] - kD MC[t] - kE MC[t] - kPT MC[t] + kTP MP[t], Derivative[1][MP][t] == kPT MC[t] - kTP MP[t]}; ic = {DC[0] == 0, DP[0] == 0, IC[0] == IC0, IP[0] == 0, MC[0] == 0, MP[0] == 0}; sol = Simplify[Solve[Join[de, ic],{mDC,bDC,mDP,bDP,mIC,bIC,mIP,bIP,mMC,bMC,mMP,bMP}]] finds a solution in a second and Join[de,ic]//.sol//Simplify returns {{True,True,True,True,True,True,True,True,True,True,True,True}}

Because the structure of the system is so simple, by inspection guess the form of one solution

DC[t_]:=mDC E^-t+bDC;
DP[t_]:=mDP E^-t+bDP;
IC[t_]:=mIC E^-t+bIC;
IP[t_]:=mIP E^-t+bIP;
MC[t_]:=mMC E^-t+bMC;
MP[t_]:=mMP E^-t+bMP;
de = {Derivative[1][DC][t] == -kE DC[t] - kPT DC[t] + kTP DP[t] + kD MC[t],
  Derivative[1][DP][t] == kPT DC[t] - kTP DP[t],
  Derivative[1][IC][t] == -2 kD IC[t] - kE IC[t] - kPT IC[t] + kTP IP[t],
  Derivative[1][IP][t] == kPT IC[t] - kTP IP[t],
  Derivative[1][MC][t] == 2 kD IC[t] - kD MC[t] - kE MC[t] - kPT MC[t] + kTP MP[t],
  Derivative[1][MP][t] == kPT MC[t] - kTP MP[t]};
ic = {DC[0] == 0, DP[0] == 0, IC[0] == IC0, IP[0] == 0, MC[0] == 0, MP[0] == 0};
sol = Simplify[Solve[Join[de, ic],{mDC,bDC,mDP,bDP,mIC,bIC,mIP,bIP,mMC,bMC,mMP,bMP}]]

finds a solution in a second and

Join[de,ic]//.sol//Simplify

returns

{{True,True,True,True,True,True,True,True,True,True,True,True}}

POSTED BY: Bill Nelson

Reply to this discussion

Reply Preview

Attachments

Remove Add a file to this post

Follow this discussion

or Discard

Feedback