2
$\begingroup$

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] == (1 - β), d[1] == 1/β, p [1] == (1 - β), ϕ[1] == 0.01, w[1] == 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...

|improve this question
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Nasser Mar 28 '15 at 5:53
  • $\begingroup$ @ Nasser... my problem is $${v'[\eta] = k1, d'[\eta] =k2, p'[\eta] = k3, \phi'[\eta] = k4, w'[\eta] = k5, v[1] = (1 - \beta), d[1] = 1/\beta, p [1] = (1 - \beta), \phi[1] = 0.01, w[1] = 0.02}$$ $\endgroup$ – user45799 Mar 28 '15 at 5:55
  • $\begingroup$ @ Nasser....where $k_i$ are function of $\eta$...given in question... $\endgroup$ – user45799 Mar 28 '15 at 6:06
  • $\begingroup$ I see a few cases of extra ] brackets thrown in in the "k2=" line. $\endgroup$ – Sjoerd C. de Vries Mar 28 '15 at 6:14
  • 2
    $\begingroup$ This post has code for Runge-Kutta solving-a-system-of-odes-with-the-runge-kutta-method $\endgroup$ – Nasser Mar 28 '15 at 6:23
5
$\begingroup$

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] == (1 - β), d[1] == 1/β, p[1] == (1 - β), ϕ[1] == 0.01, w[1] == 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)[[1]]["Domain"];
 Plot[{v@n, d@n, p@n, ϕ@n, w@n} /. sol // Evaluate, {n, nl, nr}]

enter image description here

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.

|improve this answer
$\endgroup$
  • $\begingroup$ @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... $\endgroup$ – user45799 Mar 28 '15 at 13:20
  • $\begingroup$ @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 $\endgroup$ – user45799 Mar 28 '15 at 13:21
  • $\begingroup$ @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. $\endgroup$ – xzczd Mar 29 '15 at 9:10
  • $\begingroup$ 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 η = ? $\endgroup$ – user45799 Mar 29 '15 at 9:39
  • $\begingroup$ @Prasanta NDSolve does stop at 0.9576, but have you noticed that it begins from 1 because your initial condition is given at η = 1? BTW, what's your first language? $\endgroup$ – xzczd Mar 29 '15 at 9:49

Not the answer you're looking for? Browse other questions tagged or ask your own question.