# 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]-Cx[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 suppose to fluctuate around the equilibrium point but in the code I have so far, this doesn't happen.

-
Welcome. Why don't you show what you have done so far? – user21 Apr 18 '13 at 4:48
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 – andre Apr 18 '13 at 7:53

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 this 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! – Isaias Cordova 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

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[0] == 1, y[0] == 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}}


ListLinePlot[{xl, yl}, Mesh -> All, MeshStyle -> PointSize[Medium]]