Ttgjyuk

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!

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

For given stepsize τ and given vectorfield 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 2 million times with NestList and still need only 2 seconds.

nsteps = 20000000;
xlist = Range[0., step nsteps, step];
Ylist = NestList[cStep, initials, nsteps]; // AbsoluteTiming // First

2.08678