# Runge-Kutta implemented on Mathematica

I am trying to solve differential equations numerically, so I am trying to write a 4th -order Runge-Kutta program for Mathematica (I know NDSolve does this, but I want to do my own). I ran into some trouble though, as my program just loops infinitely.

RK[a_,b_,y0_,n_,f_]:= Module[{},
h=(b-a)/n;
X = Table[a+k*h, {k,0,n} ];
Y = Table[y0, {k,0,n} ];
For[j=1, j<n, j++,
k1 = f[X[[j]],Y[[j]]];
k2 = f[X[[j]]+(h/2),Y[[j]]+h*(k1/2)];
k3 = f[X[[j]]+(h/2),Y[[j]]+h*(k2/2)];
k4 = f[X[[j+1]],Y[[j]]+h*k3];
Y[[j+1]]= Y[[j]]+(h/6)(k1+2*k2+2*k3+k4);
];
Return[Transpose[{X,Y}]];
];


I don't think my issue is with the algortithm though... I think its with my definition of the differential equation. I was honestly pretty lost on how I do this, but this is what I came up with:

f[x_,y_] = y - (x^2)(y)^2;
RK[0,10,2,50,f[x,Function[x,y[x]]]]


I tried defining it as a function of two variables... but I think I might have done some thing wrong.

If this is wrong...how do I define a differential equation as a function of two variables?

• there are many problems in your code. First you do not find k4. Second, you need ; between end of the For loops and the Return[] statement. Also you do not need an explicit Return. Also the call is wrong. Why not just RK[0, 2, 2, 5, f[x, x]] ? Try to correct these first and see. btw, there is lots of RK4 code in this forum, many questions were asked about it before. If you google, you'll find examples. Apr 22, 2020 at 17:26
• The k4 I just messed up in post...its correct in my code. Ill fix the others, thank you. Apr 22, 2020 at 17:28
• I am not sure how f[x,x] is equivalent...or correct. Wont that just solve x - x^4? Apr 22, 2020 at 17:30
• You also failed to localize any variables within Module. Apr 22, 2020 at 17:34

This works for me

RK[a_, b_, y0_, n_, f_] := Module[{X, Y, j, k1, k2, k3, k4, h},
h = (b - a)/n;
X = Table[a + k*h, {k, 0, n}];
Y = Table[y0, {k, 0, n}];
For[j = 1, j < n, j++, k1 = f[X[[j]], Y[[j]]];
k2 = f[X[[j]] + (h/2), Y[[j]] + h*(k1/2)];
k3 = f[X[[j]] + (h/2), Y[[j]] + h*(k2/2)];
k4 = f[X[[j + 1]], Y[[j]] + h*k3];
Y[[j + 1]] = Y[[j]] + (h/6) (k1 + 2*k2 + 2*k3 + k4);
];

Transpose[{X, Y}]
];

f[x_, y_] := y - (x^2) (y)^2;
RK[0, 2, 2, 5, f] // N


• Nice fix on the output... the first one destroyed my browser! Thanks for being so helpful! Apr 22, 2020 at 17:47

Nasser already pointed out many mistakes, so I won't go into those.

NestList would allow for a much cleaner implementation.

Below, RK4step[f,h] denotes a function which takes a pair of $$\{t,y(t)\}$$ values, and produces the next one at $$t+h$$, assuming that $$y'(t) = f(t, y(t))$$.

ClearAll[RK4step]
RK4step[f_, h_][{t_, y_}] :=
Module[{k1, k2, k3, k4},
k1 = f[t,       y];
k2 = f[t + h/2, y + h k1/2];
k3 = f[t + h/2, y + h k2/2];
k4 = f[t + h,   y + h k3];
{t + h, y + h/6*(k1 + 2 k2 + 2 k3 + k4)}
]


We can use NestList to take a starting pair $$\{t_0, y(t_0)\}$$, and repeatedly propagate the time using RK4step.

res =
NestList[
RK4step[-#2 &, 0.1], (* #2 & is short for f where f[t_, y_] := -y, look up Function *)
{0.0, 1.0}, (* this is {t0, y(t0)} *)
100 (* compute this many steps *)
]

ListPlot[res, PlotRange -> All]


More complex example, a harmonic oscillator:

f[t_, {x_, v_}] := {v, -x}

res = NestList[
RK4step[f, 0.1],
{0.0, {1.0, 0.0}},
100
];

ListPlot[
Transpose[{res[[All, 1]], res[[All, 2, 1]]}]
]

• Thanks for the spunked up version. I am a Mathematica Ignoramus, so I am gonna stick with Nasser's fix on my code. Yours looks cool, but its beyond me! Thanks for being so helpful! Apr 22, 2020 at 17:43
• @ragnvaldr.js Recommended reading: mathematica.stackexchange.com/q/134609/12 Apr 22, 2020 at 17:44
• Thanks... I had no clue that the For loop should be avoided in Mathematica. For what I am doing, this will work, but Ill try and work through your solution for future use... Thanks! Apr 22, 2020 at 17:49
• I'm sure it is... but I'm an amatuer Java Developer so the For loop is just second nature to me. I guess I am confused as the res = ... stuff. The added info helps... The syntax is still just unfamiliar. Apr 22, 2020 at 17:56
• @Szabolcs Your code is perfect (+1), see my answer on mathematica.stackexchange.com/questions/221848/… May 14, 2020 at 22:32