How to speed up the simulation time of a big system of ODEs?

Hi Everyone,

I have a big system of equations (378 diff. equations) and I have to simulate it using an ODE solver (I am using ODE23s) but the simulation time is too long.

Is there a way to improve my code in order to have a small simulation time?

Attached you can find the code.

Thank you!

Gabri

Answer 1

[J. Alex Lee]

don't use syms

[Gabriele Galli]

Thank you Alex,

My ODE system is obtained using symbolic steps hence the final matrix, the one that is the output of the function that I pass to the ODE solver is also symbolic. How can I convert change this matrix (Sxd) into a non symbol one? Because in my code I tried to convert it using double() (Sxd=double(Sxd)) but the simulation is still endless. Thank you!

[J. Alex Lee]

i would start from the paper-and-pencil and see if you can define equations numerically, and also vectorially (your variables suggest that). What do your equations look like if you write them down?

[Gabriele Galli]

Hi Alex,

Unfortunately using a paper-and-pencil approach for this problem is very difficult due to its dimension.

These 378 equations are obtained starting from a system of just 18 equations (and 18 variables) and 21 parameters and performing partial differentiation of each equation w.r.t. each of the 21 parameters that is why the final system is composed by 21x18=378 equations.

So far I solved this problem manually but that system was much smaller and easy to handle than this one. That is why, for this one I need to use symbolic calculations (or, maybe, there are other approaches that I don't know).

By the way, I though to convert my symbolic matrix into a double one (see code attached here) but, unfortunately, it contains the variables x1 through x378 that are not numerical values hence I cannot convert it.

Due you have any idea either on how I can perform this conversion or other approaches in order to obtain a faster simulation?

Thank you so much for your help!

[J. Alex Lee]

well, if you insist on the symbolics approach, i'm afraid i can't help because i don't know how to use the symbolic toolbox (and don't have it).

if you want to step back, the only anyone can help you is if you show the original 18 equations, and why/how you are taking partial derivatives w.r.t. parameters.

[Gabriele Galli]

Yes, sure. I mean maybe we can find a better solution. I am not familiar either with the symbolic toolbox. The equations are contained in the f matrix below:

r1 = k1*TGFb_TGFbR;
r2 = k2*TGFbR*TGFb;
r3 = k3*TGFb_TGFbR*(1-exp(-((t-k20)/(k21))^10));  
r4 = k4*TGFb_TGFbR_P;
r5 = k5*TGFb_TGFbR_P*I_Smad;
r6 = k6*I_Smad_TGFb_TGFbR_P;
r7 = k7*Smad*TGFb_TGFbR_P;
r8 = k8*Smad;
r9 = k9*Smad_N;
r10 = k10*2*Smad_P*Smad_P;
r11 = k11*Smad_P_Smad_P;
r12 = k10*Smad_P*CoSmad;
r13 = k11*Smad_P_CoSmad;
r14 = k8*CoSmad;
r15 = k9*CoSmad_N;
r16 = k12*k8*Smad_P_Smad_P;
r17 = k8*Smad_P;
r18 = k9*Smad_P_N;
r19 = k12*k8*Smad_P_CoSmad;
r20 = k13*Smad_P_N;
r21 = k10*2*Smad_P_N*Smad_P_N;
r22 = k11*Smad_P_Smad_P_N;
r23 = k10*Smad_P_N*CoSmad_N;
r24 = k11*Smad_P_CoSmad_N;
r25 = k14*Smad_P_CoSmad_N^2/(Smad_P_CoSmad_N^2 + k15^2);
r26 = k16*I_Smad_mRNA1;
r27 = k17*I_Smad_mRNA2;
r28 = k18*I_Smad_mRNA2;
r29 = k19*I_Smad;

f=[r1 - r2;
r1 - r2;
- r1 + r2 - r3 + r4 + r6;
r3 - r4 - r5;
r5 - r6;
- r7 - r8 + r9;
r7 - r10 + r11 - r12 + r13 - r17 + r18;
-r12 + r13 - r14 + r15;
r10 - r11 - r16;
r12 - r13 - r19;
r8 - r9 + r20;
r16 + r21 - r22;
r17 - r18 - r20 - r21 + r22 - r23 + r24;
r19 + r23 - r24;
r14 - r15 - r23 + r24;
r25 - r26;
r26 - r27;
r28 - r29 - r5 + r6];

In particular, k1 through k21 are the parameters with nominal values known and TGFb_TGFbR through I_Smad are the variables, for semplicity I called them x1 thorough x18.

What I want to find is the sensitivity matrix Sxd:

Sxd=dfdx*Sx+dfdk;
%dfdx is the Jacobian of f wrt x(=vector contained the 18 parameters), dim(dfdx=18x18)
%dfdk= is the Jacobian of f wrt k(=vector contained the 21 parameters), dim(dfdk=18x21)
%Sx=sensitivity matrix=dx(t,k)/dk. This is what I want to calculate using the ODE...
... I for simplicity I created I called each of these x1, x2, ..., x378 and reshape...
... them in a matrix 18x21. So basically, each line of the matrix contains the sensitivity...
... of a variable x wrt to each of the 21 parameters. 
%Sxd is the first derivative wrt t of Sx. dim(Sxd)=18x21. For this reason, I reshape...
... this matrix in a column vector (dim=378x1) in order to be a suitable input for the ODE solver.

The time span I want to simulate is [0, 180] and we can assume all the 378 initial conditions = 0 except for:

 IC(1)=1; IC(2)=1; IC(6)=40.98; IC(8)=34.15; IC(11)=19.02; IC(15)=15.85;

Please, let me know if you need more information.

Thank you!