# Parallelizing NDSolve

I have a system of 300 coupled linear ODE's, which takes days to solve with NDSolve (i.e. but it's basically not possible to solve it since it requires >100 gigabytes of RAM). In other posts, I read that NDSolve isn't really the best option. Since these posts were a while ago, I was wondering if in the latest versions of Mathematica there are some more capabilities/improvements to solve large systems of ODE's more efficiently (e.g. by running it in parallel)?

If yes, can I get an example code?

Here is the system reduced to 10 equations:

eqn={Derivative[1][l1][t] ==
1. l1[t] - 0.0176237 l10[t] - 0.166071 l2[t] - 0.0331033 l3[t] -
0.0851187 l4[t] - 0.0293385 l5[t] - 0.142254 l6[t] -
0.133462 l7[t] - 0.0721599 l8[t] - 0.185757 l9[t],
Derivative[1][l2][t] == -0.166071 l1[t] - 0.167883 l10[t] +
1. l2[t] - 0.117562 l3[t] - 0.028769 l4[t] - 0.00579195 l5[t] -
0.125756 l6[t] - 0.0949809 l7[t] - 0.000568083 l8[t] -
0.0762299 l9[t],
Derivative[1][l3][t] == -0.0331033 l1[t] - 0.103511 l10[t] -
0.117562 l2[t] + 1. l3[t] - 0.0923666 l4[t] - 0.106864 l5[t] -
0.0459526 l6[t] - 0.194075 l7[t] - 0.017699 l8[t] - 0.126495 l9[t],
Derivative[1][l4][t] == -0.0851187 l1[t] - 0.0060876 l10[t] -
0.028769 l2[t] - 0.0923666 l3[t] + 1. l4[t] - 0.0913889 l5[t] -
0.164316 l6[t] - 0.15856 l7[t] - 0.189563 l8[t] - 0.199029 l9[t],
Derivative[1][l5][t] == -0.0293385 l1[t] - 0.095501 l10[t] -
0.00579195 l2[t] - 0.106864 l3[t] - 0.0913889 l4[t] + 1. l5[t] -
0.113142 l6[t] - 0.198038 l7[t] - 0.142897 l8[t] - 0.201798 l9[t],
Derivative[1][l6][t] == -0.142254 l1[t] - 0.211227 l10[t] -
0.125756 l2[t] - 0.0459526 l3[t] - 0.164316 l4[t] -
0.113142 l5[t] + 1. l6[t] - 0.0224247 l7[t] - 0.0135337 l8[t] -
0.163697 l9[t],
Derivative[1][l7][t] == -0.133462 l1[t] - 0.208206 l10[t] -
0.0949809 l2[t] - 0.194075 l3[t] - 0.15856 l4[t] -
0.198038 l5[t] - 0.0224247 l6[t] + 1. l7[t] - 0.0811846 l8[t] -
0.0446819 l9[t],
Derivative[1][l8][t] == -0.0721599 l1[t] - 0.0941197 l10[t] -
0.000568083 l2[t] - 0.017699 l3[t] - 0.189563 l4[t] -
0.142897 l5[t] - 0.0135337 l6[t] - 0.0811846 l7[t] + 1. l8[t] -
0.100951 l9[t],
Derivative[1][l9][t] == -0.185757 l1[t] - 0.202612 l10[t] -
0.0762299 l2[t] - 0.126495 l3[t] - 0.199029 l4[t] -
0.201798 l5[t] - 0.163697 l6[t] - 0.0446819 l7[t] -
0.100951 l8[t] + 1. l9[t],
Derivative[1][l10][t] == -0.0176237 l1[t] + 1. l10[t] -
0.167883 l2[t] - 0.103511 l3[t] - 0.0060876 l4[t] -
0.095501 l5[t] - 0.211227 l6[t] - 0.208206 l7[t] -
0.0941197 l8[t] - 0.202612 l9[t],
Derivative[1][x1][t] == -l1[t] - l10[t] - l2[t] - l3[t] - l4[t] -
l5[t] - l6[t] - l7[t] - l8[t] - l9[t] - 1. x1[t] +
0.0176237 x10[t] + 0.166071 x2[t] + 0.0331033 x3[t] +
0.0851187 x4[t] + 0.0293385 x5[t] + 0.142254 x6[t] +
0.133462 x7[t] + 0.0721599 x8[t] + 0.185757 x9[t],
Derivative[1][x2][t] == -l1[t] - l10[t] - l2[t] - l3[t] - l4[t] -
l5[t] - l6[t] - l7[t] - l8[t] - l9[t] + 0.166071 x1[t] +
0.167883 x10[t] - 1. x2[t] + 0.117562 x3[t] + 0.028769 x4[t] +
0.00579195 x5[t] + 0.125756 x6[t] + 0.0949809 x7[t] +
0.000568083 x8[t] + 0.0762299 x9[t],
Derivative[1][x3][t] == -l1[t] - l10[t] - l2[t] - l3[t] - l4[t] -
l5[t] - l6[t] - l7[t] - l8[t] - l9[t] + 0.0331033 x1[t] +
0.103511 x10[t] + 0.117562 x2[t] - 1. x3[t] + 0.0923666 x4[t] +
0.106864 x5[t] + 0.0459526 x6[t] + 0.194075 x7[t] +
0.017699 x8[t] + 0.126495 x9[t],
Derivative[1][x4][t] == -l1[t] - l10[t] - l2[t] - l3[t] - l4[t] -
l5[t] - l6[t] - l7[t] - l8[t] - l9[t] + 0.0851187 x1[t] +
0.0060876 x10[t] + 0.028769 x2[t] + 0.0923666 x3[t] - 1. x4[t] +
0.0913889 x5[t] + 0.164316 x6[t] + 0.15856 x7[t] +
0.189563 x8[t] + 0.199029 x9[t],
Derivative[1][x5][t] == -l1[t] - l10[t] - l2[t] - l3[t] - l4[t] -
l5[t] - l6[t] - l7[t] - l8[t] - l9[t] + 0.0293385 x1[t] +
0.095501 x10[t] + 0.00579195 x2[t] + 0.106864 x3[t] +
0.0913889 x4[t] - 1. x5[t] + 0.113142 x6[t] + 0.198038 x7[t] +
0.142897 x8[t] + 0.201798 x9[t],
Derivative[1][x6][t] == -l1[t] - l10[t] - l2[t] - l3[t] - l4[t] -
l5[t] - l6[t] - l7[t] - l8[t] - l9[t] + 0.142254 x1[t] +
0.211227 x10[t] + 0.125756 x2[t] + 0.0459526 x3[t] +
0.164316 x4[t] + 0.113142 x5[t] - 1. x6[t] + 0.0224247 x7[t] +
0.0135337 x8[t] + 0.163697 x9[t],
Derivative[1][x7][t] == -l1[t] - l10[t] - l2[t] - l3[t] - l4[t] -
l5[t] - l6[t] - l7[t] - l8[t] - l9[t] + 0.133462 x1[t] +
0.208206 x10[t] + 0.0949809 x2[t] + 0.194075 x3[t] +
0.15856 x4[t] + 0.198038 x5[t] + 0.0224247 x6[t] - 1. x7[t] +
0.0811846 x8[t] + 0.0446819 x9[t],
Derivative[1][x8][t] == -l1[t] - l10[t] - l2[t] - l3[t] - l4[t] -
l5[t] - l6[t] - l7[t] - l8[t] - l9[t] + 0.0721599 x1[t] +
0.0941197 x10[t] + 0.000568083 x2[t] + 0.017699 x3[t] +
0.189563 x4[t] + 0.142897 x5[t] + 0.0135337 x6[t] +
0.0811846 x7[t] - 1. x8[t] + 0.100951 x9[t],
Derivative[1][x9][t] == -l1[t] - l10[t] - l2[t] - l3[t] - l4[t] -
l5[t] - l6[t] - l7[t] - l8[t] - l9[t] + 0.185757 x1[t] +
0.202612 x10[t] + 0.0762299 x2[t] + 0.126495 x3[t] +
0.199029 x4[t] + 0.201798 x5[t] + 0.163697 x6[t] +
0.0446819 x7[t] + 0.100951 x8[t] - 1. x9[t],
Derivative[1][x10][t] == -l1[t] - l10[t] - l2[t] - l3[t] - l4[t] -
l5[t] - l6[t] - l7[t] - l8[t] - l9[t] + 0.0176237 x1[t] -
1. x10[t] + 0.167883 x2[t] + 0.103511 x3[t] + 0.0060876 x4[t] +
0.095501 x5[t] + 0.211227 x6[t] + 0.208206 x7[t] +
0.0941197 x8[t] + 0.202612 x9[t], x1[0] == 1, x2[0] == 1,
x3[0] == 0, x4[0] == 0, x5[0] == 1, x6[0] == 0, x7[0] == 0,
x8[0] == 1, x9[0] == 1, x10[0] == 0, x1[10] == 0, x2[10] == 1,
x3[10] == 1, x4[10] == 0, x5[10] == 1, x6[10] == 0, x7[10] == 1,
x8[10] == 1, x9[10] == 1, x10[10] == 0}


This is the code for NDSolve:

   NDSolve[eqn, {x1[t], x2[t], x3[t], x4[t], x5[t], x6[t], x7[t], x8[t],
x9[t], x10[t], l1[t], l2[t], l3[t], l4[t], l5[t], l6[t], l7[t],
l8[t], l9[t], l10[t]}, {t, 0, 10}]

• What type are these ODE's? Maybe just linear? Oct 21, 2017 at 16:19
• There is really not enough information here to attempt an informed response. I suggest posting a smaller example, say 10 variables, that is representative in terms of the type of equations, time axis, method(s) used, etc. Oct 21, 2017 at 16:23
• The example appears to be a BVP and is not so nicely behaved, which may account for the trouble in ramping to a much larger system. I will suggest having a look at the [NDSolve advanced docs](tutorial/NDSolveBVP) if you've not already done so. Also I did notice that setting Method->"Chasing" sped things noticeably (possibly at expense of quality though). Oct 21, 2017 at 19:36
• @DanielLichtblau: Ah yes, forgot to mention it is in deed a BVP. Thank very much, using the Chasing-method sped up things considerably and I can integrate the 300 equations now! Oct 21, 2017 at 21:18
• There was a talk on this at WTC 2017 this week. You should try to get the notebook. Oct 22, 2017 at 13:54

The example appears to be a BVP and is not so nicely behaved, which may account for the trouble in ramping to a much larger system. I will suggest having a look at the NDSolve advanced docs if you've not already done so. Also I did notice that setting Method -> "Chasing" sped things noticeably (possibly at expense of quality though). – Daniel Lichtblau Oct 21 '17 at 19:36