Solve a linear system of ODEs?

Questions
, ,
1

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];

debug.nb (129 KB)

2

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}}
3

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

4

DSolve may be working very hard on finding more complicated solutions and given enough time it might find those for you. There are other simple solutions, like each function being a scalar*Sin[t] except for IC which needs to be scalar1*Sin[t]+scalar2*Cos[t]. I expect there are complex exponential solutions. I expect there are more complicated solutions than the simple ones that I can guess.

5

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.

6

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.

7

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
8

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
9

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