# How to make Runge-Kutta for solving nonlinear ODE system in Mathematica

How to make Runge Kutta method for system of non linear equations in this Matrix form. Matrices A3, B3 and B4 are functions of matrix X3. Initial conditions are X1=X2=X3=X4=0 and system is here

$$\frac{d}{dt}\left( \begin{array}{c} \text{X1} \\ \text{X2} \\ \text{X3} \\ \text{X4} \end{array} \right)=\left( \begin{array}{cccc} -\text{A2} & -\text{A1} & -\text{A3} & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ -\text{B3} & 0 & -\text{B2}-\text{B4} & \text{B1} \end{array} \right).\left( \begin{array}{c} \text{X1} \\ \text{X2} \\ \text{X3} \\ \text{X4} \end{array} \right)+\left( \begin{array}{c} 0 \\ 0 \\ 0 \\ \text{B5} \end{array} \right).\left( \begin{array}{cccc} 0 & 0 & 0 & 1 \end{array} \right)$$

$$\text{X1}=\left( \begin{array}{c} \text{x11} \\ \text{x12} \\ \text{x13} \end{array} \right);$$

$$\text{X2}=\left( \begin{array}{c} \text{x21} \\ \text{x22} \\ \text{x23} \end{array} \right);$$

$$\text{X3}=\left( \begin{array}{c} \text{x31} \\ \text{x32} \\ \text{x33} \end{array} \right);$$

$$\text{X4}=\left( \begin{array}{c} \text{x41} \\ \text{x42} \\ \text{x43} \end{array} \right);$$

Where we have matrices

$$\text{A1}=\left( \begin{array}{ccc} \frac{640576501799}{76356} & 0 & \frac{311934265235}{57839} \\ 0 & \frac{1285802795705}{40871} & 0 \\ \frac{388065734765}{30838} & 0 & \frac{523620702496}{6935} \end{array} \right);$$

$$\text{A2}=\left( \begin{array}{ccc} \frac{980000000000000}{1168149} & 0 & \frac{210000000000000}{389383} \\ 0 & \frac{1225000000000000}{389383} & 0 \\ \frac{490000000000000}{389383} & 0 & \frac{2940000000000000}{389383} \end{array} \right);$$

$$\text{A3}=\left( \begin{array}{ccc} 0 & \frac{152600000000000000 \text{x31}}{14532941709}+\frac{76300000000000000 \text{x33}}{4844313903} & 0 \\ 0 & \frac{953750000000000000 \text{x32}}{14532941709} & 0 \\ 0 & \frac{76300000000000000 \text{x31}}{4844313903}+\frac{1068200000000000000 \text{x33}}{4844313903} & 0 \end{array} \right);$$

$$\text{B1}=\left( \begin{array}{ccc} \frac{5874951206}{12955399} & 0 & \frac{3212738087}{4414317} \\ 0 & \frac{8215070163}{2795011} & 0 \\ \frac{8400845503}{6412667} & 0 & \frac{6246501886}{520165} \end{array} \right);$$

$$\text{B2}=\left( \begin{array}{ccc} \frac{14850000000000}{327471103} & 0 & \frac{71500000000000}{982413309} \\ 0 & \frac{96250000000000}{327471103} & 0 \\ \frac{42900000000000}{327471103} & 0 & \frac{393250000000000}{327471103} \end{array} \right);$$

$$\text{B3}=\left( \begin{array}{ccc} 0 & -\frac{180000000000000000 \text{x31}}{1027581737}-\frac{40000000000000000 \text{x33}}{79044749} & 0 \\ 0 & -\frac{700000000000000000 \text{x32}}{440392173} & 0 \\ 0 & -\frac{40000000000000000 \text{x31}}{79044749}-\frac{660000000000000000 \text{x33}}{79044749} & 0 \end{array} \right);$$

$$\text{B4}=\left( \begin{array}{ccc} \frac{55687500000000000000 \left(\frac{16 \text{x31}^2}{11781}-\frac{640 \text{x31} \text{x33}}{2909907}\right)}{327471103}+\frac{160875000000000000000 \left(\frac{16 \text{x31} \text{x33}}{11781}-\frac{640 \text{x33}^2}{2909907}\right)}{327471103} & \frac{540900000000000000000 \text{x31} \text{x32}}{9625358130479}+\frac{120200000000000000000 \text{x32} \text{x33}}{740412163883} & \frac{55687500000000000000 \left(-\frac{640 \text{x31}^2}{2909907}+\frac{2456 \text{x31} \text{x33}}{14549535}\right)}{327471103}+\frac{160875000000000000000 \left(-\frac{640 \text{x31} \text{x33}}{2909907}+\frac{2456 \text{x33}^2}{14549535}\right)}{327471103} \\ \frac{505312500000000000000 \left(\frac{16 \text{x31} \text{x32}}{11781}-\frac{640 \text{x32} \text{x33}}{2909907}\right)}{327471103} & \frac{2103500000000000000000 \text{x32}^2}{4125153484491} & \frac{505312500000000000000 \left(-\frac{640 \text{x31} \text{x32}}{2909907}+\frac{2456 \text{x32} \text{x33}}{14549535}\right)}{327471103} \\ \frac{160875000000000000000 \left(\frac{16 \text{x31}^2}{11781}-\frac{640 \text{x31} \text{x33}}{2909907}\right)}{327471103}+\frac{2654437500000000000000 \left(\frac{16 \text{x31} \text{x33}}{11781}-\frac{640 \text{x33}^2}{2909907}\right)}{327471103} & \frac{120200000000000000000 \text{x31} \text{x32}}{740412163883}+\frac{1983300000000000000000 \text{x32} \text{x33}}{740412163883} & \frac{160875000000000000000 \left(-\frac{640 \text{x31}^2}{2909907}+\frac{2456 \text{x31} \text{x33}}{14549535}\right)}{327471103}+\frac{2654437500000000000000 \left(-\frac{640 \text{x31} \text{x33}}{2909907}+\frac{2456 \text{x33}^2}{14549535}\right)}{327471103} \end{array} \right);$$

$$\text{B5}=\left( \begin{array}{c} \frac{37125000}{3241} \\ 0 \\ \frac{107250000}{3241} \end{array} \right);$$

-
@Georde NDSolve has the method built-in. NDSolve[{x1'[t] == -x1[t] x2[t]^2, x2'[t] == -x1[t] (1 + x2[t]^2), x1[0] == 1, x2[0] == 1}, {x1, x2}, {t, 0, 2}, Method -> "ExplicitRungeKutta"]. See ref-page online, and look for possible settings of Method options – Sasha Nov 15 '11 at 3:00
@Sasha, why not post this as an answer so we can vote for it? – Mr.Wizard Nov 15 '11 at 4:15
@Georde, also the NDSolve plug-ins tutorial gives explicit detail of setting up an RK4 integration. Admittedly, not all of it is germane to your question, but it does lay out the algorithm for you. – rcollyer Nov 15 '11 at 4:27
@George You should better spend your time thinking about this extremely simple problem, which admits a closed form solution, that seek to solve it using Runge-Kutta. The solution can be worked out by hand, in terms of matrix exponentials and integrals theoreof. You could then use software to compute those. – Sasha Nov 16 '11 at 18:40
@George: You really should learn how to use Mathematica for the simple examples before you try to apply it to your real work. I've been helping you since August and neither your ability to program nor your ability to ask clear questions on a forum has increased. People want to help, but are not going to do ALL of your work for you. Make your questions either interesting or short. Long, messy, stupid and frustrating questions will get downvoted or ignored. – Simon Nov 18 '11 at 3:03

Function NDSolve in Mathematica has the Runge-Kutta methods built-in:

NDSolve[{x1'[t] == -x1[t] x2[t]^2, x2'[t] == -x1[t] (1 + x2[t]^2),
x1[0] == 1, x2[0] == 1}, {x1, x2}, {t, 0, 2},
Method -> "ExplicitRungeKutta"]


According to NDSolve's online reference page, method "ExplicitRungeKutta" gives explicit Runge-Kutta methods with adaptive embedded pairs of 2(1) through 9(8), method "ImplicitRungeKutta" gives families of arbitrary-order implicit Runge-Kutta methods.

Order of the Runge-Kutta method can be set as follows:

NDSolve[{x1'[t] == -x1[t] x2[t]^2, x2'[t] == -x1[t] (1 + x2[t]^2),
x1[0] == 1, x2[0] == 1}, {x1, x2}, {t, 0, 2},
Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 7}]


NDSolve permits writing your own method, and an example of implementing classic Runge-Kutta method of order 4 is given in online tutorial "NDSolve Method Plugin Framework".

Added: This example solve a vector system in NDSolve using Runge-Kutta method:

xvec = {x1[t], x2[t], x3[t]};
x0vec = {1, 2, 3};
amat = {{0, -x3[t], x2[t]}, {-x3[t], 0, x1[t]}, {1, 1, 1}};
NDSolve @@ {Flatten@{Thread[D[xvec, t] == amat.xvec],
Thread[x0vec == (xvec /. t -> 0)]}, {x1, x2, x3}, {t, 0, 2},
Method -> "ExplicitRungeKutta"}

-
Yes, but I have matrices in basic equations. Code is different and I want to apply reference.wolfram.com/mathematica/tutorial/NDSolvePlugIns.html but the problem is how for my case? – George Nov 15 '11 at 12:54
@George I have added an example that solve the vector equation using the built-in method. – Sasha Nov 15 '11 at 13:13
Thank you Sasha, but can you make me code like here but exactly for my case reference.wolfram.com/mathematica/tutorial/NDSolvePlugIns.html because I don't understand programming – George Nov 15 '11 at 14:17
@George, I just downloaded your notebook, and in it you're using the form \[DifferentialD]f/\[DifferentialD]t which is incorrect as DifferentialD has no meaning by itself. Instead, you're looking for \[PartialD]. – rcollyer Nov 16 '11 at 19:22
@George You might want to gain a modicum of familiarity with the tools you're using before attacking a problem head on. Else you'll not know whether to trust the output (because you can't trust your input) and you will spend a lot of time fretting over small issues. – user4423 Nov 17 '11 at 18:57