Monday, September 9, 2013

Differential Equations #11: Numerical Solutions with the Runge-Kutta 4th Order Method


Update:  One of the numerical example (example 1) has incorrect numerical results.  This blog has been edited.  The results posted are now correct. Apologizes. -  Eddie
1/4/2015

Introduction

In this blog we will explore one of the ways to solve the differential equation y' = f(x,y) by numerical methods. One of the common methods is the Runge-Kutta 4th Order Method, probably the most common methods.

With the initial condition (x_n, y_n) and step size h, the next point (x_n+1, y_n+1) is calculated by the following equations:
x_n+1 = x_n + h
y_n+1 = y_n + (k1 + 2*k2 + 2*k3 + k4)/6 with

k1 = h * f(x_n, y_n)
k2 = h * f(x_n + h/2, y_n + k1/2)
k3 = h * f(x_n + h/2, y_n + k2/2)
k4 = h * f(x_n + h, y_n + k3)

(Almost) obviously use either a calculator or mathematical software for this process.

Examples:
1. dy/dx = -0.23 * (y + 1)
Initial Condition: (0, 24)
We want y1 when x1 = 1 and y2 when x2 = 2.


Observe that h = 1 with x0 = 0 and y0 = 24. Note the lack of x term in f(x,y).

Then:
x1 = 0 + 1 = 1
k1 ≈ -5.75
k2 ≈ -5.08875
k3 ≈ -5.16479375
k4 ≈ -4.562097438
y1 ≈ 18.86346918

x2 = 1 + 1 = 2
k1 ≈ -4.568597911
k2 ≈ -4.043209151
k3 ≈ -4.103628858
k4 ≈ -3.624763273
y2 ≈ 14.78229631

So our points are:
(0, 24)
(1, 18.86346918)
(2, 14.78229631)


2. y' = y - x with I.C. y(0) = 2
Find y(0.1) and y(0.2).


Note h = 0.1 with x0 = 0, y0 = 2, and f(x,y) = -x + y.

x1 = 0.1
k1 = 0.1 * f(0, 2) = 0.2
k2 = 0.1 * f(0.5 * 0.1, 2 + 0.5 * 0.2) = 0.205
k3 = 0.1 * f(0.5 * 0.1, 2 + 0.5 * 0.205) = 0.20525
k4 = 0.1 * f(0.1, 2 + 0.20525) = 0.210525
y1 = 2 + (0.2 + 2 * 0.205 + 2 * 0.20525 + 0.210525)/6 ≈ 2.205171

x2 = 0.2
k1 ≈ 0.210517
k2 ≈ 0.216043
k3 ≈ 0.216319
k4 ≈ 0.222149
y2 ≈ 2.421203

So our points are:
(0, 2)
(0.1, 2.205171)
(0.2, 2.421203)


Talk to you all next time!
Eddie

This blog is property of Edward Shore. 2013

No comments:

Post a Comment

HP Prime and Casio fx-CG50: Trapezoid Midsegment, Height, Area

HP Prime and Casio fx-CG50:  Trapezoid Midsegment, Height, Area The program TRAPEZ calculates the following: Midsegment leng...