Solving a system of ODEs with the Runge-Kutta method

I´m trying to solve a system of ODEs using a fourth-order Runge-Kutta method. I have to recreate certain results to obtain my degree. But I'm a beginner at Mathematica programming and with the Runge-Kutta method as well.

{A = 0.30, B = 1, C = 40, D = 1, E = 0.75, F = 0.11,
r = 2.5, a = 2, e = 0.475, g = 2, d = 0.03, n = 0.01, p = -0.00005}

x'[t]/x[t] = (A + B x[t] - C x[t]^2 - F)/(D + E^(r y[t]));
y'[t]/y[t] = g (((a s[t] x[t] k[t])/m[t]) - e) (1 - y[t]);
m'[t]/m[t] = n;
k'[t]/k[t] = x[t] - d;
s'[t]/s[t] = -p;

I'd appreciate any kind of help. For over a month now, I've tried to solve this system myself but have only gotten bad results. This model is supposed to fluctuate around the equilibrium point, but in the code I have so far, this doesn't happen.

• Please don't use single capital letters as variables: both C and E are protected system symbols, and no value should be assigned to them. Also note that multiplication requires at least a space, thus Cx is a symbol while C x is C times x – István Zachar Apr 18 '13 at 7:50
• Runge-Kutta 4 is described in the documentation as a example here – andre314 Apr 18 '13 at 7:53

Here is a functional approach. The following will give you one step of the Runge-Kutta formula:

RungeKutta[func_List, yinit_List, y_List, step_] :=
Module[{k1, k2, k3, k4},
k1 = step N[func /. MapThread[Rule, {y, yinit}]];
k2 = step N[func /. MapThread[Rule, {y, k1/2 + yinit}]];
k3 = step N[func /. MapThread[Rule, {y, k2/2 + yinit}]];
k4 = step N[func /. MapThread[Rule, {y, k3 + yinit}]];
yinit + Total[{k1, 2 k2, 2 k3, k4}]/6]

Here, func is a list of functions, yinit a list of initial values, y a list of function variables (in your case that will be {x, y, m, k, s}, and step is the step size of the numerical simulation. You can then use something like NestList to iterate it many times like this:

NestList[RungeKutta[func, #, y, step] &, N[yinit], Round[t/step]]

Here, t is the maximum value of your independent variable. You can also include conditionals to check that the lists are of equal length.

• I appreciate your help, I have achieved my simulation. With the huge help of your answer! – user6983 Apr 19 '13 at 0:34
• @IsaiasCordova, glad I could help. – RunnyKine Apr 19 '13 at 0:51
• Er… I think it's better to mention that the y for RungeKutta[] usually include t for more general cases. (And the counterpart of t is 1 in yinit.) – xzczd Jun 26 '13 at 10:51
• @IsaiasCordova Can you please post the complete answer here? – TMH Apr 9 '14 at 14:11
• @xzczd, yes, that version as written in this answer is intended for autonomous equations; Maeder's book implemented it this way initially, but then also gave a version for non-autonomous DEs. – J. M. will be back soon Mar 28 '18 at 0:49

andre has linked you to the method plug-in framework in the comments, but there is a more direct way to implement classical Runge-Kutta, just by supplying its Butcher tableau to "ExplicitRungeKutta". Here's an adaptation from the docs:

ClassicalRungeKuttaCoefficients[4, prec_] := With[{amat = {{1/2}, {0, 1/2}, {0, 0, 1}},
bvec = {1/6, 1/3, 1/3, 1/6}, cvec = {1/2, 1/2, 1}}, N[{amat, bvec, cvec}, prec]]

{xf, yf} = {x, y} /. First @
NDSolve[{x'[t] == -y[t], y'[t] == x[t], x == 1, y == 0},
{x, y}, {t, 0, 6},
Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4,
"Coefficients" -> ClassicalRungeKuttaCoefficients},
StartingStepSize -> 1/2];

To obtain the values computed by RK, you can then do this:

xl = MapThread[Append, {xf["Grid"], xf["ValuesOnGrid"]}]
{{0., 1.}, {0.5, 0.877604}, {1., 0.540588}, {1.5, 0.0714256}, {2., -0.415108},
{2.5, -0.800012}, {3., -0.989166}, {3.5, -0.936349}, {4., -0.65453}, {4.5, -0.212684},
{5., 0.281088}, {5.5, 0.706007}, {6., 0.95816}}

yl = MapThread[Append, {yf["Grid"], yf["ValuesOnGrid"]}]
{{0., 0.}, {0.5, 0.479167}, {1., 0.841037}, {1.5, 0.99713}, {2., 0.90931},
{2.5, 0.599108}, {3., 0.142441}, {3.5, -0.348969}, {4., -0.754924}, {4.5, -0.976153},
{5., -0.958587}, {5.5, -0.706572}, {6., -0.281796}}

Plot them:

ListLinePlot[{xl, yl}, Mesh -> All, MeshStyle -> PointSize[Medium]] • Hey, I have a question. why do they call the vector on the inner part of the butchear tableau amac and not avec? I don't quite see why they call it that. Not a big issue, but just conufsed by it. It is because it's a matrix and not a vector? and if so, why doesn't it have for example, as the entries, {{1/2,0,0},{0,1/2,0},{0,0,1}}? Is it because we add the factors in as necessary, so exclude and 0s after the last entry for that step of the algorithm? – Shinaolord Oct 13 at 18:00