# 4th-order Runge-Kutta method to solve a system of coupled ODEs [duplicate]

I am a beginner at Mathematica programming and with the Runge-Kutta method as well. I'm trying to solve a system of coupled ODEs using a 4th-order Runge-Kutta method for my project work.

I have solved it by NDSolve, but I want to solve this by 4th-order Runge-Kutta method. Here is my problem:

Γ = 1.4
k = 0
z = 0
β = 0.166667

k1 = (d[η] v[η] η (1 - z d[η]) (v[η] - η) - 2 p[η] η (1 - z d[η]) - ϕ[η]^2 d[η]
(1 - z d[η]) - Γ p[η] v[η])/((Γ p[η] - (v[η] - η)^2 d[η] (1 - z d[η])) η)

k2 = (d[η] (1 - z d[η]) (v[η] d[η] (v[η] - 2 η) (v[η] - η) + 2 p[η] η + ϕ[η]^2
d[η]))/((Γ p[η] - (v[η] - η)^2 d[η] (1 - z d[η])) (v[η] - η) η)

k3 = (p[η] d[η] (2 η (v[η] - η)^2 (1 - z d[η]) + Γ v[η] (v[η] - 2 η) (v[η] - η) +
ϕ[η]^2 Γ))/((Γ p[η] - (v[η] - η)^2 d[η] (1 - z d[η])) (v[η] - η) η)

k4 = -((ϕ[η] (v[η] + η))/(η (v[η] - η)))

k5 = -(w[η]/(η (v[η] - η)))

sol = NDSolve[{v'[η] == k1, d'[η] == k2, p'[η] == k3, ϕ'[η] == k4, w'[η] == k5,
v == (1 - β), d == 1/β, p  == (1 - β), ϕ == 0.01, w == 0.02},
{v, d, p, ϕ, w}, {η, 0, 1}, MaxSteps -> 30000]


Please guide me how can I solve the above problem with 4th-order Runge-Kutta method, thanks.

code for RK4 method are given in Solving a system of ODEs with the Runge-Kutta method

but how can I apply those codes to my problem...please guide me...

• Can you show the equations themselves you are trying to solve with the initial/boundary conditions? If you do not know latex, you can scan them from the book and paste the image. – Nasser Mar 28 '15 at 5:53
• @ Nasser... my problem is $${v'[\eta] = k1, d'[\eta] =k2, p'[\eta] = k3, \phi'[\eta] = k4, w'[\eta] = k5, v = (1 - \beta), d = 1/\beta, p  = (1 - \beta), \phi = 0.01, w = 0.02}$$ – user45799 Mar 28 '15 at 5:55
• This post has code for Runge-Kutta solving-a-system-of-odes-with-the-runge-kutta-method – Nasser Mar 28 '15 at 6:23
• After pasting your code here, select them and press Ctrl+K, or add four whitespaces before each line of code. Things like \[Eta] isn't unreadable for Mathematica user, when pasting back to the notebook, it'll get back to the normal greek letters, and if you want to make your code look better in this site (of course it's encouraged), use the link I gave above to convert those "unreadable" symbol to greek letters! – xzczd Mar 28 '15 at 8:35
• If RunnyKine doesn't appear in the comment under this question, then the "@" won't work. And why are you still adding whitespace between "@" and the name? – xzczd Mar 28 '15 at 9:59

## 1 Answer

According to your statement, I think what you need is just 4th-order Runge-Kutta method, and a completely self-made implementation of 4th-order Runge-Kutta method isn't necessary, then the answer from J.M. has shown you the optimal direction:

(* Unchanged part omitted. *)

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

sol = NDSolve[{v'[η] == k1, d'[η] == k2, p'[η] == k3, ϕ'[η] == k4, w'[η] == k5,
v == (1 - β), d == 1/β, p == (1 - β), ϕ == 0.01, w == 0.02},
{v, d, p, ϕ, w}, {η, 0, 1},
Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4,
"Coefficients" -> ClassicalRungeKuttaCoefficients}, StartingStepSize -> 1/10000]


However, what I really want to point out is, despite the above code seems to solve your ODE set up to η = 0.0001, I'm afraid it's not reliable at all:

 {{nl, nr}} = (v /. sol)[]["Domain"];
Plot[{v@n, d@n, p@n, ϕ@n, w@n} /. sol // Evaluate, {n, nl, nr}] NDSolve by default setting doesn't manage to solve this set of equation, too. It stopped at about η = 0.9576. (I'm not sure what do you mean by saying you have solved it by NDSolve.) I'm not surprised though, your ODEs are non-linear. As for how to solve the ODEs, it's another question. I vote to close this question as a duplicate.

• @xzczd...if answer is not reliable ...please tell me how should I solve my problem....I want to solve this problem numerically for my project... – user45799 Mar 28 '15 at 13:20
• @xzczd...when I try to solve the above problem by 'Ndsolve'...is there any command to find which method mathematica has used to solve it – user45799 Mar 28 '15 at 13:21
• @Prasanta Since I'm the author of this post, you don't need to add "@xzczd" to remind me. In my view the most vital part of your problem is, you haven't even find a way to solve your nonlinear ODE. It's easy to create nonlinear ODE, but solving it can be really hard. I suggest you to first make sure if your equations are correct (incorrect translation for the real problem isn't rare!), and consider carefully if you really need your equations to be so complicated. Simplify your model if possible. Diving into the hell of looking for proper combination of options of NDSolve is the last choice. – xzczd Mar 29 '15 at 9:10
• range of η can be between o to 1....as NDSolve stopped at about η = 0.9576. while solving the system.....is there any command to find when we solve the above system by RK4 where it stopped solving .i.e. about η = ? – user45799 Mar 29 '15 at 9:39
• @Prasanta If I correctly understood those material I found, RK4 i.e. classical RK doesn't have an error estimator. Notice the document also gives a example for a fourth order RK called Fehlberg method, which owns an error estimator and can find the singularity at 0.9576. – xzczd Apr 1 '15 at 5:43