10
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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 $\endgroup$ – István Zachar Apr 18 '13 at 7:50
  • 2
    $\begingroup$ Runge-Kutta 4 is described in the documentation as a example here $\endgroup$ – andre314 Apr 18 '13 at 7:53
18
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ I appreciate your help, I have achieved my simulation. With the huge help of your answer! $\endgroup$ – user6983 Apr 19 '13 at 0:34
  • $\begingroup$ @IsaiasCordova, glad I could help. $\endgroup$ – RunnyKine Apr 19 '13 at 0:51
  • 2
    $\begingroup$ 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.) $\endgroup$ – xzczd Jun 26 '13 at 10:51
  • $\begingroup$ @IsaiasCordova Can you please post the complete answer here? $\endgroup$ – TMH Apr 9 '14 at 14:11
  • 2
    $\begingroup$ @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. $\endgroup$ – J. M. will be back soon Mar 28 '18 at 0:49
14
$\begingroup$

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}}

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]]

plot of approximate solution from RK

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy