38
$\begingroup$

The original ODE I had was $$ \frac{d^2y}{dx^2}+\frac{dy}{dx}-6y=0$$ with $y(0)=3$ and $y'(0)=1$. Now I can solve this by hand and obtain that $y(1) = 14.82789927$. However I wish to use the 4th order Runge-Kutta method, so I have the system:

$$ \left\{\begin{array}{l} \frac{dy}{dx} = z \\ \frac{dz}{dx} = 6y - z \end{array}\right. $$ With $y(0)=3$ and $z(0)=1$.

Now I know that for two general 1st order ODE's $$ \frac{dy}{dx} = f(x,y,z) \\ \frac{dz}{dx}=g(x,y,z)$$ The 4th order Runge-Kutta formula's for a system of 2 ODE's are: $$ y_{i+1}=y_i + \frac{1}{6}(k_0+2k_1+2k_2+k_3) \\ z_{i+1}=z_i + \frac{1}{6}(l_0+2l_1+2l_2+l_3) $$ Where $$k_0 = hf(x_i,y_i,z_i) \\ k_1 = hf(x_i+\frac{1}{2}h,y_i+\frac{1}{2}k_0,z_i+\frac{1}{2}l_0) \\ k_2 = hf(x_i+\frac{1}{2}h,y_i+\frac{1}{2}k_1,z_i+\frac{1}{2}l_1) \\ k_3 = hf(x_i+h,y_i+k_2,z_i+l_2) $$ and $$l_0 = hg(x_i,y_i,z_i) \\ l_1 = hg(x_i+\frac{1}{2}h,y_i+\frac{1}{2}k_0,z_i+\frac{1}{2}l_0) \\ l_2 = hg(x_i+\frac{1}{2}h,y_i+\frac{1}{2}k_1,z_i+\frac{1}{2}l_1) \\ l_3 = hg(x_i+h,y_i+k_2,z_i+l_2)$$

My problem is I am struggling to apply this method to my system of ODE's so that I can program a method that can solve any system of 2 first order ODE's using the formulas above, I would like for someone to please run through one step of the method, so I can understand it better.

$\endgroup$
1
  • $\begingroup$ For reference, see this answer on SO. $\endgroup$ Commented Mar 24, 2018 at 15:22

5 Answers 5

Reset to default
49
$\begingroup$

I will outline the process and you can fill in the calculations.

We have our system as:

$$ \left\{\begin{array}{l} \frac{dy}{dx} = z \\ \frac{dz}{dx} = 6y - z \end{array}\right. $$ With $y(0)=3$ and $z(0)=1$.

We must do the calculations in a certain order as there are dependencies between the numerical calculations. This order is:

  • $k_0 = hf(x_i,y_i,z_i)$

  • $l_0 = hg(x_i,y_i,z_i)$

  • $k_1 = hf(x_i+\frac{1}{2}h,y_i+\frac{1}{2}k_0,z_i+\frac{1}{2}l_0)$

  • $l_1 = hg(x_i+\frac{1}{2}h,y_i+\frac{1}{2}k_0,z_i+\frac{1}{2}l_0)$

  • $k_2 = hf(x_i+\frac{1}{2}h,y_i+\frac{1}{2}k_1,z_i+\frac{1}{2}l_1)$

  • $l_2 = hg(x_i+\frac{1}{2}h,y_i+\frac{1}{2}k_1,z_i+\frac{1}{2}l_1)$

  • $k_3 = hf(x_i+h,y_i+k_2,z_i+l_2)$

  • $l_3 = hg(x_i+h,y_i+k_2,z_i+l_2)$

  • $y_{i+1}=y_i + \frac{1}{6}(k_0+2k_1+2k_2+k_3)$

  • $z_{i+1}=z_i + \frac{1}{6}(l_0+2l_1+2l_2+l_3)$

We typically need some inputs for the algorithm:

  • A range that we want to do the calculations over: $a \le t \le b$, lets use $a = 0, b = 1$.
  • The number of steps $N$, say $N = 10$.
  • The steps size $h = \dfrac{b-a}{N} = \dfrac{1}{10}$

The system we are solving is:

$$ \frac{dy}{dx} = f(x,y,z) = z \\ \frac{dz}{dx}=g(x,y,z) = 6y - z$$

Doing the calculations using the above order for the first time step $i= 0, t_0 = 0 = x_0$, yields:

  • $k_0 = hf(x_0,y_0,z_0) = \dfrac{1}{10}(z_0) = \dfrac{1}{10}(1) = \dfrac{1}{10}$

  • $l_0 = hg(x_0,y_0,z_0) = \dfrac{1}{10}(6y_0 - z_0) = \dfrac{1}{10}(6 \times 3 - 1) = \dfrac{1}{10}(17)$

  • $k_1 = hf(x_0+\frac{1}{2}h,y_0+\frac{1}{2}k_0,z_0+\frac{1}{2}l_0) = \dfrac{1}{10}(1 + \dfrac{1}{2}\dfrac{1}{10}(17)) ~~$(You please continue the calcs.)

  • $l_1 = hg(x_0+\frac{1}{2}h,y_i+\frac{1}{2}k_0,z_0+\frac{1}{2}l_0)$

  • $k_2 = hf(x_0+\frac{1}{2}h,y_0+\frac{1}{2}k_1,z_0+\frac{1}{2}l_1)$

  • $l_2 = hg(x_0+\frac{1}{2}h,y_0+\frac{1}{2}k_1,z_0+\frac{1}{2}l_1)$

  • $k_3 = hf(x_0+h,y_0+k_2,z_0+l_2)$

  • $l_3 = hg(x_0+h,y_0+k_2,z_0+l_2)$

  • $y_{1}=y_0 + \frac{1}{6}(k_0+2k_1+2k_2+k_3)$

  • $z_{1}=z_0 + \frac{1}{6}(l_0+2l_1+2l_2+l_3)$

You now have $x_1$ and $z_1$ which you need for the next time step after all of the intermediate (in order again). Now, we move on to the next time step $i = 1, t_1 = t_0 + h = \dfrac{1}{10} = x_1$, so we have:

  • $k_0 = hf(x_1,y_1,z_1) = \dfrac{1}{10}(z_1)$

  • $l_0 = hg(x_1,y_1,z_1) = \dfrac{1}{10}(6y_1 - z_1)$

  • $k_1 = hf(x_1+\frac{1}{2}h,y_1+\frac{1}{2}k_0,z_1+\frac{1}{2}l_0)$

  • $l_1 = hg(x_1+\frac{1}{2}h,y_1+\frac{1}{2}k_0,z_1+\frac{1}{2}l_0)$

  • $k_2 = hf(x_1+\frac{1}{2}h,y_1+\frac{1}{2}k_1,z_1+\frac{1}{2}l_1)$

  • $l_2 = hg(x_1+\frac{1}{2}h,y_1+\frac{1}{2}k_1,z_1+\frac{1}{2}l_1)$

  • $k_3 = hf(x_1+h,y_1+k_2,z_1+l_2)$

  • $l_3 = hg(x_1+h,y_1+k_2,z_1+l_2)$

  • $y_{2}=y_1 + \frac{1}{6}(k_0+2k_1+2k_2+k_3)$

  • $z_{2}=z_1 + \frac{1}{6}(l_0+2l_1+2l_2+l_3)$

Continue this for $10$ time steps. Your final result should match closely (assuming the numerical algorithm is stable for this problem) to the exact solution. You will compare $z_{10}$ to the exact result. The exact solution is:

$$y(x) = e^{-3 x}+2 e^{2 x}$$

If we find $y(1) = \dfrac{1}{e^3} + 2 e^2 = 14.8278992662291643974401973...$.

$\endgroup$
7
  • $\begingroup$ Thanks for your thorough response, seeing the start of it I now understand it better. Thanks! $\endgroup$ Commented Mar 21, 2014 at 15:49
  • 1
    $\begingroup$ Recall, in your new system, the first equation $y' = z$ is just a dummy variable in order to use RK4 methods. Regards $\endgroup$ Commented Mar 21, 2014 at 15:53
  • 1
    $\begingroup$ @Michael: Also, you will clearly see when you calculate $y_i, z_i$, which is the correct final result. $\endgroup$ Commented Mar 21, 2014 at 16:02
  • 1
    $\begingroup$ Late reply but, is $y_i$ or $z_i$ the solution to the original ODE? Comparing values it seems like the solution is given by $y_i$, but I'm not sure. $\endgroup$ Commented May 4, 2016 at 15:47
  • 1
    $\begingroup$ @ErikVesterlund if I got this right, then z would be the solution for the derivative of y and y is the solution to the original ODE $\endgroup$ Commented Jan 22, 2017 at 16:11
23
$\begingroup$

Although this answer contains the same content as Amzoti's answer, I think it's worthwhile to see it another way.

In general consider if you had $m$ first-order ODE's (after appropriate decomposition). The system looks like

\begin{align*} \frac{d y_1}{d x} &= f_1(x, y_1, \ldots, y_m) \\ \frac{d y_2}{d x} &= f_2(x, y_1, \ldots, y_m) \\ &\,\,\,\vdots\\ \frac{d y_m}{d x} &= f_m(x, y_1, \ldots, y_m) \\ \end{align*}

Define the vectors $\vec{Y} = (y_1, \ldots, y_m)$ and $\vec{f} = (f_1, \ldots, f_m)$, then we can write the system as

$$\frac{d}{dx} \vec{Y} = \vec{f}(x,\vec{Y})$$

Now we can generalize the RK method by defining \begin{align*} \vec{k}_1 &= h\vec{f}\left(x_n,\vec{Y}(x_n)\right)\\ \vec{k}_2 &= h\vec{f}\left(x_n + \tfrac{1}{2}h,\vec{Y}(x_n) + \tfrac{1}{2}\vec{k}_1\right)\\ \vec{k}_3 &= h\vec{f}\left(x_n + \tfrac{1}{2}h,\vec{Y}(x_n) + \tfrac{1}{2}\vec{k}_2\right)\\ \vec{k}_4 &= h\vec{f}\left(x_n + h, \vec{Y}(x_n) + \vec{k}_3\right) \end{align*}

and the solutions are then given by $$\vec{Y}(x_{n+1}) = \vec{Y}(x_n) + \tfrac{1}{6}\left(\vec{k}_1 + 2\vec{k}_2 + 2\vec{k}_3 + \vec{k}_4\right)$$

with $m$ initial conditions specified by $\vec{Y}(x_0)$. When writing a code to implement this one can simply use arrays, and write a function to compute $\vec{f}(x,\vec{Y})$

For the example provided, we have $\vec{Y} = (y,z)$ and $\vec{f} = (z, 6y-z)$. Here's an example in Fortran90:

program RK4 implicit none integer , parameter :: dp = kind(0.d0) integer , parameter :: m = 2 ! order of ODE real(dp) :: Y(m) real(dp) :: a, b, x, h integer :: N, i ! Number of steps N = 10 ! initial x a = 0 x = a ! final x b = 1 ! step size h = (b-a)/N ! initial conditions Y(1) = 3 ! y(0) Y(2) = 1 ! y'(0) ! iterate N times do i = 1,N Y = iterate(x, Y) x = x + h end do print*, Y contains ! function f computes the vector f function f(x, Yvec) result (fvec) real(dp) :: x real(dp) :: Yvec(m), fvec(m) fvec(1) = Yvec(2) !z fvec(2) = 6*Yvec(1) - Yvec(2) !6y-z end function ! function iterate computes Y(t_n+1) function iterate(x, Y_n) result (Y_nplus1) real(dp) :: x real(dp) :: Y_n(m), Y_nplus1(m) real(dp) :: k1(m), k2(m), k3(m), k4(m) k1 = h*f(x, Y_n) k2 = h*f(x + h/2, Y_n + k1/2) k3 = h*f(x + h/2, Y_n + k2/2) k4 = h*f(x + h, Y_n + k3) Y_nplus1 = Y_n + (k1 + 2*k2 + 2*k3 + k4)/6 end function end program

This can be applied to any set of $m$ first order ODE's, just change m in the code and change the function f to whatever is appropriate for the system of interest. Running this code as-is yields

$ 14.827578509968953 \qquad 29.406156886687729$

The first value is $y(1)$, the second $z(1)$, correct to the third decimal point with only ten steps.

$\endgroup$
5
  • 1
    $\begingroup$ you should use $x_{n}$ instead of $t_{n}$ $\endgroup$ Commented Nov 3, 2018 at 20:39
  • $\begingroup$ Good catch, fixed it $\endgroup$ Commented Nov 4, 2018 at 0:20
  • 1
    $\begingroup$ Fantastic answer @Kai +1. Would give +50 if possible! Many people struggle with systems of ODE's and RK methods. I have a question though regarding your Fortran implementation. If you wanted to be fancy you could write your $k_i$'s using a for loop correct? Essentially placing them in an array? So you would have an array $k(i,n)$ where i was the number of stages and n was the dimension of your state vector? Are you aware of any documentation that does this in Fortran? I am writing something similar at the minute and am a bit stumped!! $\endgroup$ Commented Feb 10, 2019 at 0:03
  • $\begingroup$ I have noticed that there is no explicit dependence of function f on variable x. Is that because of the substitution we used to reduce the order of derivative? Will this dependence on x be automatically picked up by the algorithm through other variables? $\endgroup$ Commented Jun 10, 2022 at 10:31
  • 1
    $\begingroup$ @MadPhysicist the ODE in the original question didn't have explicit $x$ dependence, it's only implicit in $y(x)$ and $z(x)$. You could have explicit $x$ dependence of course $\endgroup$ Commented Jun 10, 2022 at 14:49
8
$\begingroup$

A Matlab implementation is given below:

% It calculates ODE using Runge-Kutta 4th order method % Author Ido Schwartz % Originally available form: http://www.mathworks.com/matlabcentral/fileexchange/29851-runge-kutta-4th-order-ode/content/Runge_Kutta_4.m % Edited by Amin A. Mohammed, for 2 ODEs(April 2016) clc; % Clears the screen clear all; h=0.1; % step size x = 0:h:1; % Calculates upto y(1) y = zeros(1,length(x)); z = zeros(1,length(x)); y(1) = 3; % initial condition z(1) = 1; % initial condition % F_xy = @(t,r) 3.*exp(-t)-0.4*r; % change the function as you desire F_xyz = @(x,y,z) z; % change the function as you desire G_xyz = @(x,y,z) 6*y-z; for i=1:(length(x)-1) % calculation loop k_1 = F_xyz(x(i),y(i),z(i)); L_1 = G_xyz(x(i),y(i),z(i)); k_2 = F_xyz(x(i)+0.5*h,y(i)+0.5*h*k_1,z(i)+0.5*h*L_1); L_2 = G_xyz(x(i)+0.5*h,y(i)+0.5*h*k_1,z(i)+0.5*h*L_1); k_3 = F_xyz((x(i)+0.5*h),(y(i)+0.5*h*k_2),(z(i)+0.5*h*L_2)); L_3 = G_xyz((x(i)+0.5*h),(y(i)+0.5*h*k_2),(z(i)+0.5*h*L_2)); k_4 = F_xyz((x(i)+h),(y(i)+k_3*h),(z(i)+L_3*h)); % Corrected L_4 = G_xyz((x(i)+h),(y(i)+k_3*h),(z(i)+L_3*h)); y(i+1) = y(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h; % main equation z(i+1) = z(i) + (1/6)*(L_1+2*L_2+2*L_3+L_4)*h; % main equation end
$\endgroup$
0
$\begingroup$

A Fortran code shown below:

enter image description here

produces the following result

output

!Runge-Kutta Fourth Order Method !For 2nd Order Differentiation Equation !First you have to define the function F(x,y,z) = z !dy/dx G(x,y,z) = 6*y-z !dz/dx = d2y/dx2 INTEGER :: n,i REAL :: k1,l1,k2,l2,k3,l3,k4,l4 !Most Important Write (*,*) "Given Equation '(y2)-6(y1)+(y0)=0'" Write (*,*) "Xo=0, Yo=3, Zo=Y'o=1, Xn=1, n=?" Xo=0 !Given Condition Yo=3 !Given Condition Zo=1 !Given Condition Xn=1 !Given Condition read (*,*) n !n=number of Intercept h=(Xn-Xo)/n do i=1,n !you have to do the Calculation 'n' times k1 = h*F(Xo,Yo,Zo) l1 = h*G(Xo,Yo,Zo) k2 = h*F(Xo+h/2,Yo+k1/2,Zo+l1/2) l2 = h*G(Xo+h/2,Yo+k1/2,Zo+l1/2) k3 = h*F(Xo+h/2,Yo+k2/2,Zo+l2/2) l3 = h*G(Xo+h/2,Yo+k2/2,Zo+l2/2) k4 = h*F(Xo+h,Yo+k3,Zo+l3) l4 = h*G(Xo+h,Yo+k3,Zo+l3) !Sum Up Yn = Yo+(k1+2*k2+2*k3+k4)/6 Zn = Zo+(l1+2*l2+2*l3+l4)/6 !Operation for Next calculation Xo=Xo+h !(+h) than previous Term Yo=Yn !Now Yn becomes Yo Zo=Zn !Now Zn becomes Zo End Do Write (*,*) "Xn,Yn =",Xo,Yo Stop End
$\endgroup$
1
  • $\begingroup$ Can your clarify your question? $\endgroup$ Commented Jan 31, 2019 at 19:39
0
$\begingroup$

A Rust implementation is below:

const M: usize = 2; // Size of the system fn main() { let mut y: [f64; M] = [0.0; M]; let mut a: f64 = 0.0; let mut x: f64 = a; let mut b: f64 = 0.0; let n: usize = 10; let h: f64 = (b - a) / n as f64; // Initial conditions y[0] = 3.0; // y(0) y[1] = 1.0; // y'(0) for _i in 0..n { y = rk4_step(x, y, h, &derivatives); x += h; } println!("{:?}", y); } // Define the derivatives for the system of equations fn derivatives(x: f64, y: [f64; M]) -> [f64; M] { let mut dy_dx: [f64; M] = [0.0; M]; dy_dx[0] = y[1]; // dy/dx = z dy_dx[1] = 6.0 * y[0] - y[1]; // dz/dx = 6y - z dy_dx } // Perform one step of RK4 fn rk4_step(x: f64, y_n: [f64; M], h: f64, f: &dyn Fn(f64, [f64; M]) -> [f64; M]) -> [f64; M] { let mut k1: [f64; M] = [0.0; M]; let mut k2: [f64; M] = [0.0; M]; let mut k3: [f64; M] = [0.0; M]; let mut k4: [f64; M] = [0.0; M]; let mut y_n_plus_1: [f64; M] = [0.0; M]; k1 = scale_vec(&f(x, y_n), h); k2 = scale_vec(&f(x + h / 2.0, add_vec(&y_n, &scale_vec(&k1, 0.5))), h); k3 = scale_vec(&f(x + h / 2.0, add_vec(&y_n, &scale_vec(&k2, 0.5))), h); k4 = scale_vec(&f(x + h, add_vec(&y_n, &k3)), h); y_n_plus_1 = add_vec( &y_n, &scale_vec( &add_vec(&k1, &add_vec(&scale_vec(&k2, 2.0), &add_vec(&scale_vec(&k3, 2.0), &k4))), 1.0 / 6.0, ), ); y_n_plus_1 } // Helper function: element-wise vector addition fn add_vec(a: &[f64; M], b: &[f64; M]) -> [f64; M] { let mut result: [f64; M] = [0.0; M]; for i in 0..M { result[i] = a[i] + b[i]; } result } // Helper function: element-wise scalar multiplication of a vector fn scale_vec(a: &[f64; M], scalar: f64) -> [f64; M] { let mut result: [f64; M] = [0.0; M]; for i in 0..M { result[i] = a[i] * scalar; } result }
$\endgroup$
1
  • $\begingroup$ Your answer could be improved with additional supporting information. Please edit to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center. $\endgroup$ Commented Dec 4, 2023 at 13:25

You must log in to answer this question.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.