# Ways to speed up user implemented RK4

So, I've implemented RK4, and I'm wondering what I can do to make it more efficient? What I've got so far is below. I wish to still record all steps. I think AppendTo is doing the most damage to the time, is there a faster alternative?

rk4[f_, variables_, valtinit_, tinit_, tfinal_, nsteps_] :=
Module[{table, xlist, ylist, step, k1, k2, k3, k4},
xlist = tinit;
step = N[(tfinal - tinit)/(nsteps)];
ylist = valtinit;
table = {{xlist, ylist}};
Table[
k1 = step* f /. MapThread[Rule, {variables, ylist}]; (*
k2 = step*f /. MapThread[Rule, {variables, k1/2 + ylist}];
k3 = step*f /. MapThread[Rule, {variables, k2/2 + ylist}];
k4 = step*f /. MapThread[Rule, {variables, k3 + ylist}];
ylist += 1/6 (k1 + 2 (k2 + k3) + k4);
xlist += step;
AppendTo[table, {xlist, ylist}];
{xlist, ylist}, nsteps];
table
];


Example Input:

funclist = {-x + y, x - y};
initials = {1, 2};
variables = {x, y};
init = 0;
final = 200;
nstep = 20000;
approx = rk4[funclist, variables, initials, init, final, nstep]//AbsoluteTiming;


{3.59932,{...}}

I'd love some suggestions!

• AppendTo is quadratic time complexity. Might be better to preallocate and set by index. Also it'll be much faster to not use Rule and instead code stuff up a little bit more explicitly. As a general rule, too, use vectorized operators. Those can be very fast. And if everything can be totally functional over "packed arrays" (look them up here) it'll be very quick too. – b3m2a1 Mar 26 '19 at 21:31
• I'll work on implementing it more explicity, this is what came to find first though. It'll require some changes to the inputs, I'll have to ponder this. And preallocating the list is a quick change that won't be an issue to do, I can't believe I forgot that's faster :(. Thanks though! – Shinaolord Mar 26 '19 at 21:33
• Shinaoloard, using Join[ {{xlist, ylist}}, Table[ k1 = step*f /. MapThread[Rule, {variables, ylist}]; k2 = step*f /. MapThread[Rule, {variables, k1/2 + ylist}]; k3 = step*f /. MapThread[Rule, {variables, k2/2 + ylist}]; k4 = step*f /. MapThread[Rule, {variables, k3 + ylist}]; ylist += 1/6 (k1 + 2 (k2 + k3) + k4); xlist += step; {xlist, ylist}, nsteps ] ] as return value is already a first step. There is no point in appending if you use a Table anyways. – Henrik Schumacher Mar 26 '19 at 21:40
• Why not just get NDSolve[] to use fourth-order Runge-Kutta to begin with? – J. M.'s discontentment Mar 27 '19 at 2:46
• @J.M.isslightlypensive I know it can, I just wanted to make sure I could actually code it myself, instead of just using options to get Mathematica to do it for me (:. Thanks for trying to help though!! – Shinaolord Mar 27 '19 at 13:27

Just to give you an impression how fast things may get when you use the right tools.

For given stepsize τ and given vector field F, this creates a CompiledFunction cStep  that computes a single Runge-Kutta step

F = X \[Function] {-Indexed[X, 2], Indexed[X, 1]};

τ = 0.01;
Block[{YY, Y, k1, k2, k3, k4},

YY = Table[CompileGetElement[Y, i], {i, 1, 2}];
k1 = τ F[YY];
k2 = τ F[0.5 k1 + YY];
k3 = τ F[0.5 k2 + YY];
k4 = τ F[k3 + YY];

cStep = With[{code = YY + (k1 + 2. (k2 + k3) + k4)/6. },
Compile[{{Y, _Real, 1}},
code,
CompilationTarget -> "C",
RuntimeOptions -> "Speed"
]
]
];


Now we can apply it 20 million times with NestList and it stills takes only about 2 seconds.

nsteps = 20000000;
xlist = Range[0., τ nsteps, τ];
Ylist = NestList[cStep, {1., 0.}, nsteps]; // AbsoluteTiming // First


2.08678

# Edit

This can be sped up even more my avoiding NestList (the loop behind it can also be compiled which saves several calls to libraries) and by utilizing that the dimension of the ODE is known at compile time. For low dimensional systems, it may be also beneficial to avoid tensor operations altogether and to perform computations in scalar registers as done below.

τ = 0.01;
cFlow = Block[{YY, Y, k1, k2, k3, k4, τ, Ylist, j},
YY = Table[CompileGetElement[Ylist, j, i], {i, 1, 2}];
k1 = τ F[YY];
k2 = τ F[0.5 k1 + YY];
k3 = τ F[0.5 k2 + YY];
k4 = τ F[k3 + YY];
With[{
code1 = (YY + (k1 + 2. (k2 + k3) + k4)/6)[],
code2 = (YY + (k1 + 2. (k2 + k3) + k4)/6)[]
},
Compile[{{Y0, _Real, 1}, {τ, _Real}, {n, _Integer}},
Block[{Ylist},
Ylist = Table[0., {n + 1}, {Length[Y0]}];
Ylist[] = Y0;
Do[
Ylist[[j + 1, 1]] = code1;
Ylist[[j + 1, 2]] = code2;
,
{j, 1, n}];
Ylist
],
CompilationTarget -> "C", RuntimeOptions -> "Speed"
]
]
];
Ylist2 = cFlow[{1., 0.}, τ, nsteps]; // AbsoluteTiming // First


1.06549

Don't be too upset by parts of the code being highlighted in red; this is on purpose.

• Damn, you definitely know how to use Mathematica A LOT more efficiently than I do. Thanks! – Shinaolord Mar 26 '19 at 22:07
• You're welcome. – Henrik Schumacher Mar 26 '19 at 22:08
• I'll have to play around with Compile, it definitely seems like a massive speed up if used correctly. – Shinaolord Mar 26 '19 at 22:08
• This is exactly the kind of thing I like to show people when they complain about the slowness of Mathematica. Of course with some cleverness in vectorized operations Compile could probably be avoided altogether if one so desired. – b3m2a1 Mar 27 '19 at 7:45
• Might be of interest: JM's blog article tpfto.wordpress.com/2019/04/03/the-runge-kutta-gill-method – anderstood Apr 3 '19 at 14:58