# 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. – Nasser Apr 22 '20 at 17:26
• The k4 I just messed up in post...its correct in my code. Ill fix the others, thank you. – ragnvaldr.js Apr 22 '20 at 17:28
• I am not sure how f[x,x] is equivalent...or correct. Wont that just solve x - x^4? – ragnvaldr.js Apr 22 '20 at 17:30
• You also failed to localize any variables within Module. – Szabolcs Apr 22 '20 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! – ragnvaldr.js Apr 22 '20 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! – ragnvaldr.js Apr 22 '20 at 17:43
• @ragnvaldr.js Recommended reading: mathematica.stackexchange.com/q/134609/12 – Szabolcs Apr 22 '20 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! – ragnvaldr.js Apr 22 '20 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. – ragnvaldr.js Apr 22 '20 at 17:56
• @Szabolcs Your code is perfect (+1), see my answer on mathematica.stackexchange.com/questions/221848/… – Alex Trounev May 14 '20 at 22:32