3

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

Community
  • 1
user45799
  • 133
  • 1
  • 6
  • 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 05:53
  • @ 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}$$ – user45799 Mar 28 '15 at 05:55
  • @ Nasser....where $k_i$ are function of $\eta$...given in question... – user45799 Mar 28 '15 at 06:06
  • I see a few cases of extra ] brackets thrown in in the "k2=" line. – Sjoerd C. de Vries Mar 28 '15 at 06:14
  • Please post Mathematica code, not Latex – Sjoerd C. de Vries Mar 28 '15 at 06:16
  • 2
    This post has code for Runge-Kutta solving-a-system-of-odes-with-the-runge-kutta-method – Nasser Mar 28 '15 at 06:23
  • @ Nasser...how can I apply those codes to my problem...please guide me... – user45799 Mar 28 '15 at 06:30
  • @ Sjoerd C. de Vries...please guide me how can I post my Mathematica code....i'm unable to doing so...If I copy-paste the codes...its not readable here... – user45799 Mar 28 '15 at 06:32
  • Select your code and press Ctrl+Shift+I before copying. And you can also use this to convert your code to make it look better. – xzczd Mar 28 '15 at 07:29
  • Just copy and paste (Ctrl+C and Ctrl+V under Windows), as you copy and paste any thing. BTW, if you add a white space before "@" and "xzczd", the reminder won't work. It should be "@xzczd". – xzczd Mar 28 '15 at 08:05
  • @xzczd ...but codes are not readable here... – user45799 Mar 28 '15 at 08:20
  • 1
    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 08:35
  • (Sigh…) Just found your code in the edit history. This time I've done the edit for you, please check what I've done carefully. – xzczd Mar 28 '15 at 08:51
  • @ Sjoerd C. de Vries...mathematica code is added to my question... – user45799 Mar 28 '15 at 09:02
  • @ RunnyKine...how can I apply code of http://mathematica.stackexchange.com/questions/23516/solving-a-system-of-odes-with-the-runge-kutta-method to my problem...please guide me... – user45799 Mar 28 '15 at 09:41
  • 1
    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 09:59

1 Answers1

7

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.

Community
  • 1
xzczd
  • 65,995
  • 9
  • 163
  • 468
  • @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
  • 1
    @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 09: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 09:39
  • @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? – xzczd Mar 29 '15 at 09:49
  • NDSolve shows that system of given ode has singularity at 0.9576.....how can we find the singularity of the system in RK4 algorithm... – user45799 Apr 01 '15 at 04:50
  • 1
    @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 01 '15 at 05:43