Mathematics 22: Lecture 11
Runge-Kutta
Dan Sloughter
Furman University
January 25, 2008
Dan Sloughter (Furman University) Mathematics 22: Lecture 11 January 25, 2008 1 / 11
Order of approximations
I One may show that the error in Eulers method is bounded by the
step-size h times a constant. We call Eulers method a first-order
method.
Dan Sloughter (Furman University) Mathematics 22: Lecture 11 January 25, 2008 2 / 11
Order of approximations
I One may show that the error in Eulers method is bounded by the
step-size h times a constant. We call Eulers method a first-order
method.
I The modified Euler method is a second-order method: the error is
bounded by a constant times h2 .
Dan Sloughter (Furman University) Mathematics 22: Lecture 11 January 25, 2008 2 / 11
Runge-Kutta method
I Consider the initial-value problem
du
= f (t, u), u(t0 ) = u0 .
dt
Dan Sloughter (Furman University) Mathematics 22: Lecture 11 January 25, 2008 3 / 11
Runge-Kutta method
I Consider the initial-value problem
du
= f (t, u), u(t0 ) = u0 .
dt
T
I Divide [t0 , t0 + T ] into N equal intervals of length h = N.
Dan Sloughter (Furman University) Mathematics 22: Lecture 11 January 25, 2008 3 / 11
Runge-Kutta method
I Consider the initial-value problem
du
= f (t, u), u(t0 ) = u0 .
dt
T
I Divide [t0 , t0 + T ] into N equal intervals of length h = N.
I Let ti = t0 + ih, i = 0, 1, 2, . . . , N.
Dan Sloughter (Furman University) Mathematics 22: Lecture 11 January 25, 2008 3 / 11
Runge-Kutta method
I Consider the initial-value problem
du
= f (t, u), u(t0 ) = u0 .
dt
T
I Divide [t0 , t0 + T ] into N equal intervals of length h = N.
I Let ti = t0 + ih, i = 0, 1, 2, . . . , N.
I Having computed u0 , u1 , . . . , ui , let
k1 = f (ti , ui )
h h
k2 = f ti + , ui + k1
2 2
h h
k3 = f ti + , ui + k2
2 2
k4 = f (ti + h, ui + hk3 ).
Dan Sloughter (Furman University) Mathematics 22: Lecture 11 January 25, 2008 3 / 11
Runge-Kutta (contd)
I Let
h
ui+1 = ui + (k1 + 2k2 + 2k3 + k4 ) .
6
Dan Sloughter (Furman University) Mathematics 22: Lecture 11 January 25, 2008 4 / 11
Runge-Kutta (contd)
I Let
h
ui+1 = ui + (k1 + 2k2 + 2k3 + k4 ) .
6
I Runge-Kutta is a fourth-order method.
Dan Sloughter (Furman University) Mathematics 22: Lecture 11 January 25, 2008 4 / 11
Example
I Consider the initial-value problem
du
= u cos(t), u(0) = 1.
dt
on the interval [0, 6].
Dan Sloughter (Furman University) Mathematics 22: Lecture 11 January 25, 2008 5 / 11
Example
I Consider the initial-value problem
du
= u cos(t), u(0) = 1.
dt
on the interval [0, 6].
I Let h = 0.1 as before.
Dan Sloughter (Furman University) Mathematics 22: Lecture 11 January 25, 2008 5 / 11
Example
I Consider the initial-value problem
du
= u cos(t), u(0) = 1.
dt
on the interval [0, 6].
I Let h = 0.1 as before.
I For the first step, we have
k1 = (1.0) cos(0) = 1.0
k2 = (1.0 + (0.05)(1.0)) cos(0.05) = 1.0486878
k3 = (1.0 + (0.05)(1.0486878)) cos(0.05) = 1.0511191
k4 = (1.0 + (0.1)(1.0511191)) cos(0.1) = 1.0995910.
Dan Sloughter (Furman University) Mathematics 22: Lecture 11 January 25, 2008 5 / 11
Example (contd)
I And so
0.1
u1 = 1.0 + (1.0 + (2)(1.0486878) + (2)(1.0511191) + 1.0995910)
6
= 1.1049867.
Dan Sloughter (Furman University) Mathematics 22: Lecture 11 January 25, 2008 6 / 11
Example (contd)
I And so
0.1
u1 = 1.0 + (1.0 + (2)(1.0486878) + (2)(1.0511191) + 1.0995910)
6
= 1.1049867.
I Note: the exact value is u(0.1) = e sin(0.1) = 1.1049868.
Dan Sloughter (Furman University) Mathematics 22: Lecture 11 January 25, 2008 6 / 11
Example (contd)
I And so
0.1
u1 = 1.0 + (1.0 + (2)(1.0486878) + (2)(1.0511191) + 1.0995910)
6
= 1.1049867.
I Note: the exact value is u(0.1) = e sin(0.1) = 1.1049868.
I Recall: with Eulers method we had u1 = 1.1 and with the modified
Euler method we had u1 = 1.1047252.
Dan Sloughter (Furman University) Mathematics 22: Lecture 11 January 25, 2008 6 / 11
Using Octave
I Runge-Kutta in Octave:
octave:1> function w = f(t,u)
> w = u*cos(t);
> endfunction
octave:2> t = [0:0.1:6];
octave:3> u(1) = 1.0;
octave:4> for i = 1:60
> k1 = f(t(i),u(i));
> k2 = f(t(i)+0.05,u(i)+0.05*k1);
> k3 = f(t(i)+0.05,u(i)+0.05*k2);
> k4 = f(t(i)+0.1,u(i)+0.1*k3);
> u(i+1) = u(i) + (0.1/6)*(k1 + 2*k2 + 2*k3 + k4);
> endfor
Dan Sloughter (Furman University) Mathematics 22: Lecture 11 January 25, 2008 7 / 11
Using Octave
I Runge-Kutta in Octave:
octave:1> function w = f(t,u)
> w = u*cos(t);
> endfunction
octave:2> t = [0:0.1:6];
octave:3> u(1) = 1.0;
octave:4> for i = 1:60
> k1 = f(t(i),u(i));
> k2 = f(t(i)+0.05,u(i)+0.05*k1);
> k3 = f(t(i)+0.05,u(i)+0.05*k2);
> k4 = f(t(i)+0.1,u(i)+0.1*k3);
> u(i+1) = u(i) + (0.1/6)*(k1 + 2*k2 + 2*k3 + k4);
> endfor
I Note: u(6) u60 = 0.75623, which is exact to 5 decimal places.
Dan Sloughter (Furman University) Mathematics 22: Lecture 11 January 25, 2008 7 / 11
Using Octave (contd)
I Comparison of exact (green) and approximate (red) solutions:
Dan Sloughter (Furman University) Mathematics 22: Lecture 11 January 25, 2008 8 / 11
Octave: lsode
I The following commands use the built-in Octave function lsode to
solve our equation:
octave:1> function w = f(u, t)
> w = u*cos(t);
> endfunction
octave:2> t = [0:0.1:6];
octave:3> u = lsode("f",1.0,t);
Dan Sloughter (Furman University) Mathematics 22: Lecture 11 January 25, 2008 9 / 11
Octave: lsode
I The following commands use the built-in Octave function lsode to
solve our equation:
octave:1> function w = f(u, t)
> w = u*cos(t);
> endfunction
octave:2> t = [0:0.1:6];
octave:3> u = lsode("f",1.0,t);
I Note: u and t are reversed in the definition of f from our notation.
Dan Sloughter (Furman University) Mathematics 22: Lecture 11 January 25, 2008 9 / 11
lsode (contd)
I The method used by Octave is an adaptive step-size method.
Dan Sloughter (Furman University) Mathematics 22: Lecture 11 January 25, 2008 10 / 11
lsode (contd)
I The method used by Octave is an adaptive step-size method.
I That is, the actual step-size (value of h) used varies as the integration
proceeds based on the behavior of the function.
Dan Sloughter (Furman University) Mathematics 22: Lecture 11 January 25, 2008 10 / 11
lsode (contd)
I The method used by Octave is an adaptive step-size method.
I That is, the actual step-size (value of h) used varies as the integration
proceeds based on the behavior of the function.
I In particular, the values in the t vector do not determine the step size,
but are there only for evaluation and plotting purposes.
Dan Sloughter (Furman University) Mathematics 22: Lecture 11 January 25, 2008 10 / 11
lsode (contd)
I The method used by Octave is an adaptive step-size method.
I That is, the actual step-size (value of h) used varies as the integration
proceeds based on the behavior of the function.
I In particular, the values in the t vector do not determine the step size,
but are there only for evaluation and plotting purposes.
I In particular, if one only wanted to know u(6), t could be specified by
t = [0:6:6], in which case u(2) is the approximation to u(6).
Dan Sloughter (Furman University) Mathematics 22: Lecture 11 January 25, 2008 10 / 11
lsode (contd)
I The method used by Octave is an adaptive step-size method.
I That is, the actual step-size (value of h) used varies as the integration
proceeds based on the behavior of the function.
I In particular, the values in the t vector do not determine the step size,
but are there only for evaluation and plotting purposes.
I In particular, if one only wanted to know u(6), t could be specified by
t = [0:6:6], in which case u(2) is the approximation to u(6).
I Or, we could just use u = lsode("f",1.0,[0 6]); to perform the
evaluation.
Dan Sloughter (Furman University) Mathematics 22: Lecture 11 January 25, 2008 10 / 11
Runge-Kutta in Maxima
I To approximate a solution in Maxima using Runge-Kutta:
Dan Sloughter (Furman University) Mathematics 22: Lecture 11 January 25, 2008 11 / 11
Runge-Kutta in Maxima
I To approximate a solution in Maxima using Runge-Kutta:
I load("dynamics")
Dan Sloughter (Furman University) Mathematics 22: Lecture 11 January 25, 2008 11 / 11
Runge-Kutta in Maxima
I To approximate a solution in Maxima using Runge-Kutta:
I load("dynamics")
I u:rk(u*cos(t),u,1.0,[t,0,6,0.1])$
Dan Sloughter (Furman University) Mathematics 22: Lecture 11 January 25, 2008 11 / 11
Runge-Kutta in Maxima
I To approximate a solution in Maxima using Runge-Kutta:
I load("dynamics")
I u:rk(u*cos(t),u,1.0,[t,0,6,0.1])$
I The resulting ordered pairs are in the variable u.
Dan Sloughter (Furman University) Mathematics 22: Lecture 11 January 25, 2008 11 / 11
Runge-Kutta in Maxima
I To approximate a solution in Maxima using Runge-Kutta:
I load("dynamics")
I u:rk(u*cos(t),u,1.0,[t,0,6,0.1])$
I The resulting ordered pairs are in the variable u.
I To plot the result: wxplot2d([discrete,u])
I Note: In the lab use, plot2d([discrete,u])
Dan Sloughter (Furman University) Mathematics 22: Lecture 11 January 25, 2008 11 / 11
