0
$\begingroup$

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); $$

  • 3
    @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](http://reference.wolfram.com/mathematica/ref/NDSolve.html) online, and look for possible settings of `Method` options2011-11-15
  • 1
    @Sasha, why not post this as an answer so we can vote for it?2011-11-15
  • 1
    @Georde, also the `NDSolve` [plug-ins tutorial](http://reference.wolfram.com/mathematica/tutorial/NDSolvePlugIns.html) 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.2011-11-15
  • 1
    @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.2011-11-16
  • 1
    @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.2011-11-18
  • 1
    @George: Finally, I keep telling you to try to reduce your problems down to the bare minimum. This is both good for you and for the forum. The act of isolating your sticking point will normally help you solve it. And if you're still stuck, you can post a nice clear question that other people can work on without having to dig through the muck for you.2011-11-18

1 Answers 1

6

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"}
  • 0
    Yes, but I have matrices in basic equations. Code is different and I want to apply http://reference.wolfram.com/mathematica/tutorial/NDSolvePlugIns.html but the problem is how for my case?2011-11-15
  • 0
    @George I have added an example that solve the vector equation using the built-in method.2011-11-15
  • 0
    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 programming2011-11-15
  • 0
    You can not help to finish this problem?2011-11-15
  • 0
    @George I do not have time to retype your matrices. Besides, your post does not state differentiation equations for $X2$ and $X3$. Could you post your matrices copy-paste-ready somewhere ?2011-11-15
  • 0
    Ok Sasha. Here you can find all of them and system http://www.sendspace.com/file/6kbeaa2011-11-15
  • 0
    @George OK, I got the input. However, the problem is not fully determined. What is lacking is expression for $\frac{\mathrm{d} X2}{\mathrm{d}t}$ and $\frac{\mathrm{d} X3}{\mathrm{d}t}$. Unless these are known, no further progress is possible.2011-11-15
  • 0
    Yes, I forgot, dX2/dt and dX3/dt are zero.2011-11-15
  • 0
    @George If $X_2(0) = 0$ and $\frac{\mathrm{d} X2}{\mathrm{d}t}(t) = 0$, it follows $X_2(t) = 0$, and similarly $X_3(t) = 0$. At this point I am signing off. I have already spent a lot more time on this post of yours than it merits. Try to do some work for yourself. Try asking your peers for help.2011-11-15
  • 0
    I can not understand what is a problem. I transform my system, take a look on beginning. Everything is known. Just to start some code, but which and how? I don't understand commands in mathematica. Take a look on a beginning, system.2011-11-15
  • 2
    @George, I just downloaded your notebook, and in it you're using the form `\[DifferentialD]f/\[DifferentialD]t` which is incorrect as [`DifferentialD`](http://reference.wolfram.com/mathematica/ref/DifferentialD.html) has no meaning by itself. Instead, you're looking for `\[PartialD]`.2011-11-16
  • 3
    @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.2011-11-17