Runge-Kutta Methods

Ryan Sullenberger  2 May 2007 ME102

Runge-Kutta Methods Fourth-Order RK Method 25.8 : Solve the following problem with the fourth-order RK method where y(0) = 4 and y`(0) = 0. Solve from x = 0 to 5 with h = 0.5. Plot the results.

DiffEQ and Initial Conditions:

d

2

dx d dx

d

2

y

+ 0.6⋅  d y  + 8⋅ y = 0

 dx  

y = y1

2

dx

2

y = −0.6

d dx

y

− 8y

d

y( 0 ) = 4

dx

=

−0.6y 1 − 8y 0

y ( 0) = 0

System of first order ODE's and IC's:

 4  Y0 :=  0 

D( x , Y) :=

 

Y

1

 

 −0.6⋅ Y1 − 8⋅ Y0 

4th-order RK solution: N:= 100

xY := rkfixed( Y0 , 0 , 5 , N , D)

(

xY = x Y

0

Y

1

)

0

xY

=

1

2

0

0

4

0

1

0.05

3.96

-1.571

2

0.1

3.844

-3.065

3

0.15

3.656

-4.454

4

0.2

3.401

-5.714

5

0.25

3.087

-6.825

6

0.3

2.721

-7.769

7

0.35

2.313

-8.532

8

0.4

1.871

-9.105

9

0.45

1.406

-9.482

10

0.5

0.926

-9.661

11

0.55

0.443

-9.644

12

0.6

-0.035

-9.439

13

0.65

-0.498

-9.054

14

0.7

-0.938

-8.502

15

0.75

-1.346

...

〈0〉

X := xY

〈1〉

S1 := xY

〈2〉

S2 := xY

Graphed solution: 10

5

S1 S2

0

−5

− 10

0

1

2

3 X

4

5

Solving the ODE by writing a program that engages fourth-order RK method:

:= 0.5

N := 10

f1( z) := z

f2( y, z)

h

( x y z ) :=

y z

0

0

x

0

:= −0.6⋅ z − 8 ⋅ y

←4 ←0 ←0

for  i ∈ 0 .. N − 1 x

i+1

← ( i + 1) ⋅ h

k11

← f1( zi)

k12

← f2( yi , zi)

y_new

← yi +

z_new

← zi +

1 2 1 2

⋅ k11 ⋅ h ⋅ k12 ⋅ h

k21

← f1( z_new)

k22 k22

← f2( y_ne y_new w , z_ne z_new w)

y_new

← yi +

z_new

← zi +

1 2 1 2

⋅ k21 ⋅ h ⋅ k22 ⋅ h

k31

← f1( z_new)

k32 k32

← f2( y_ne y_new w , z_ne z_new w)

y_new

← yi + k31⋅ h

z_new

← zi + k32⋅ h

k41

← f1( z_new)

k42 k42

y_new w , z_ne z_new w) ← f2( y_ne

y

← yi +

z

i+1

← zi + i+ 1

(x y z)

1 6 1 6

⋅ ( k11 + 2⋅ k21 + 2 ⋅k31 + k41 ) ⋅ h ⋅ ( k12 + 2⋅ k22 + 2 ⋅k32 + k42 ) ⋅ h

Programs solutions: 0

x

=

0

0

0

0

0

4

0

0

1

0.5

1

1.037

1

-9.289

2

1

2

-2.428

2

-3.197

3

1.5

3

-1.557

3

5.365

4

2

4

1.154

4

4.072

5

2.5

5

1.481

5

-2.333

6

3

6

-0.293

6

-3.638

7

3.5

7

-1.132

7

0.372

8

4

8

-0.185

8

2.66

9

4.5

9

0.724

9

0.656

10

0.378

10

-1.626

10

y

5

=

z

Compare to MathCAD's rkfixed solutions:

Programs solutions graphed:

10

10 y

5

5

z

y z

=

0

S1 S2

−5 − 10

0

1

2

3 x

4

5

0

−5 − 10

0

1

2

3

x, x, X , X

4

5