8
$\begingroup$

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}]; (* 
    Equivalent to step* f/.Thread[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!

$\endgroup$
  • 3
    $\begingroup$ 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. $\endgroup$ – b3m2a1 Mar 26 at 21:31
  • $\begingroup$ 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! $\endgroup$ – Shinaolord Mar 26 at 21:33
  • $\begingroup$ 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. $\endgroup$ – Henrik Schumacher Mar 26 at 21:40
  • 3
    $\begingroup$ Why not just get NDSolve[] to use fourth-order Runge-Kutta to begin with? $\endgroup$ – J. M. is away Mar 27 at 2:46
  • 1
    $\begingroup$ @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!! $\endgroup$ – Shinaolord Mar 27 at 13:27
17
$\begingroup$

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[Compile`GetElement[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[Compile`GetElement[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)[[1]],
     code2 = (YY + (k1 + 2. (k2 + k3) + k4)/6)[[2]]
     },
    Compile[{{Y0, _Real, 1}, {τ, _Real}, {n, _Integer}},
     Block[{Ylist},
      Ylist = Table[0., {n + 1}, {Length[Y0]}];
      Ylist[[1]] = 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.

$\endgroup$
  • 1
    $\begingroup$ Damn, you definitely know how to use Mathematica A LOT more efficiently than I do. Thanks! $\endgroup$ – Shinaolord Mar 26 at 22:07
  • $\begingroup$ You're welcome. $\endgroup$ – Henrik Schumacher Mar 26 at 22:08
  • $\begingroup$ I'll have to play around with Compile, it definitely seems like a massive speed up if used correctly. $\endgroup$ – Shinaolord Mar 26 at 22:08
  • 1
    $\begingroup$ 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. $\endgroup$ – b3m2a1 Mar 27 at 7:45
  • 1
    $\begingroup$ Might be of interest: JM's blog article tpfto.wordpress.com/2019/04/03/the-runge-kutta-gill-method $\endgroup$ – anderstood Apr 3 at 14:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.