RKMethods.mw
Runge-Kutta Methods
An IVP in Standard Form
The standard form for an IVP is the system of
-order ODEs, with initial conditions:
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_2.gif) |
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_3.gif) |
![diff(Y(t), t) = F(t, Y(t))](images/RKMethods_4.gif) |
![Y(t[0]) = Y0](images/RKMethods_5.gif) |
(1) |
where
is a vector of functions and
is a vector function of the scalar
and the vector
.
A 4th-order IVP
Let us use the following problem to test our numerical methods.
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_10.gif) |
![`+`(diff(diff(diff(diff(y(t), t), t), t), t), `*`(2, `*`(t, `*`(sin(t), `*`(diff(diff(diff(y(t), t), t), t))))), `-`(`/`(`*`(diff(y(t), t)), `*`(t)))) = `*`(t, `*`(exp(`+`(`-`(`*`(4, `*`(t)))))))](images/RKMethods_11.gif) |
(2) |
with initial conditions
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_12.gif) |
![y(1) = 0, (D(y))(1) = .5, ((`@@`(D, 2))(y))(1) = .1, ((`@@`(D, 3))(y))(1) = 0](images/RKMethods_13.gif) |
(3) |
The above inital value problem can be converted to a first-order system as follows. Introduce a new set of dependent variables:
> |
![Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi(](images/RKMethods_14.gif) |
In terms of the new functions
the original
-order ODE becomes the following first-order system.
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_21.gif)
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_22.gif) |
with initial conditions
> |
![Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi(](images/RKMethods_27.gif) |
![y[1](1) = 0, y[2](1) = .5, y[3](1) = .1, y[4](1) = 0](images/RKMethods_28.gif) |
(6) |
This now takes the form of a standard first-order system:
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_29.gif) |
![diff(Y(t), t) = F(t, Y)](images/RKMethods_30.gif) |
![Y(1) = Y1](images/RKMethods_31.gif) |
(7) |
where
is a vector of dimension 4 and
is a vector function of dimension 4 taking as arguments
and the vector
. Specifically,
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_36.gif) |
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_37.gif) |
![Typesetting:-mrow(Typesetting:-mverbatim(](images/RKMethods_38.gif) |
(8) |
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_39.gif) |
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_40.gif) |
![Typesetting:-mrow(Typesetting:-mverbatim(](images/RKMethods_41.gif) |
(9) |
The initial conditions are specified by the vector
:
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_43.gif) |
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_44.gif) |
![Typesetting:-mrow(Typesetting:-mverbatim(](images/RKMethods_45.gif) |
(10) |
Comparison Values
For comparison values when testing our numerical methods, apply the high-order method dverk78 to the original
-order ODE.
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_47.gif) |
![`+`(diff(diff(diff(diff(y(t), t), t), t), t), `*`(2, `*`(t, `*`(sin(t), `*`(diff(diff(diff(y(t), t), t), t))))), `-`(`/`(`*`(diff(y(t), t)), `*`(t)))) = `*`(t, `*`(exp(`+`(`-`(`*`(4, `*`(t)))))))](images/RKMethods_48.gif) |
(11) |
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_49.gif) |
![y(1) = 0, (D(y))(1) = .5, ((`@@`(D, 2))(y))(1) = .1, ((`@@`(D, 3))(y))(1) = 0](images/RKMethods_50.gif) |
(12) |
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_51.gif) |
![proc (x_dverk78) local res, data, vars, solnproc, outpoint, ndsol, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; `:=`(_EnvDSNumericSaveDigits, Digits); `:=`(Digits, 14); ...](images/RKMethods_52.gif) |
(13) |
For example:
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_53.gif) |
![4.2](images/RKMethods_54.gif) |
(14) |
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_55.gif) |
![[t = 4.20000000000000016, y(t) = 2.69988538740093586, diff(y(t), t) = 1.72666950368022687, diff(diff(y(t), t), t) = 2.79161521800559286, diff(diff(diff(y(t), t), t), t) = 15.5902337896826674]](images/RKMethods_56.gif)
![[t = 4.20000000000000016, y(t) = 2.69988538740093586, diff(y(t), t) = 1.72666950368022687, diff(diff(y(t), t), t) = 2.79161521800559286, diff(diff(diff(y(t), t), t), t) = 15.5902337896826674]](images/RKMethods_57.gif) |
(15) |
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_58.gif) |
![2.69988538740093586](images/RKMethods_59.gif) |
(16) |
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_60.gif) |
Runge-Kutta Method of Order 1
We have already seen a basic numerical time-stepping method of order 1, namely the Forward Euler Method:
> |
![Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi(](images/RKMethods_62.gif) |
![Y[`+`(n, 1)] = `+`(Y[n], `*`(h, `*`(F(t[n], Y[n]))))](images/RKMethods_63.gif) |
(17) |
In the class of Runge-Kutta methods discussed below, the Forward Euler Method is a Runge-Kutta method of order 1.
Recall that the above formula is derived from the Taylor series expansion
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_64.gif) |
![`≈`(Y(`+`(t[n], h)), series(`+`(Y(t[n]), `*`((D(Y))(t[n]), `*`(h)))+O(`^`(h, 2)),h,2))](images/RKMethods_65.gif) |
(18) |
and by dropping the term
we arrive at the Forward Euler Method.
It follows that this method has local truncation error
which leads to a global truncation error
. Hence the method is order 1.
Test Order of Accuracy: RK1
The following procedure defines the method RK1 (otherwise known as the Forward Euler Method).
The IVP to be solved is
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_105.gif) |
![proc (t, Y) options operator, arrow; Vector([Y[2], Y[3], Y[4], `+`(`-`(`*`(2, `*`(t, `*`(sin(t), `*`(Y[4]))))), `/`(`*`(Y[2]), `*`(t)), `*`(t, `*`(exp(`+`(`-`(`*`(4, `*`(t))))))))]) end proc](images/RKMethods_106.gif) |
(19) |
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_107.gif) |
![Vector[column](%id = 151434492)](images/RKMethods_108.gif) |
(20) |
Compute the solution a number of times specified by
, dividing the stepsize
by 2 each time.
For each stepsize
, compare the computed result at
of the
component of
(which corresponds to the original
), with the value
.
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_117.gif) |
![10](images/RKMethods_118.gif) |
(21) |
> |
![Typesetting:-mrow(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_119.gif) |
![1.0, 2.0](images/RKMethods_120.gif) |
(22) |
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_121.gif) |
![.4](images/RKMethods_122.gif) |
(23) |
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_123.gif) |
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_131.gif) |
![0.3906250000e-3](images/RKMethods_132.gif) |
(24) |
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_133.gif) |
![2560](images/RKMethods_134.gif) |
(25) |
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_135.gif) |
![Vector[column](%id = 150408840)](images/RKMethods_136.gif) |
(26) |
Compute successive ratios
> |
![Typesetting:-mrow(Typesetting:-mo(](images/RKMethods_137.gif) |
The ratio is approximately
.
Conclusion: Since
with
we conclude that RK1 is of order 1.
A Runge-Kutta Method of Order 2
We have previously developed a numerical time-stepping formula with a higher order, namely the Modified Euler Method which has order 2.
The Modified Euler Method can be stated in the following form, which is a common form for presenting Runge-Kutta formulas.
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_150.gif) |
> |
![Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi(](images/RKMethods_151.gif) |
Recall that the above formula is based on the Taylor series expansion
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_155.gif) |
![`≈`(Y(`+`(t[n], h)), series(`+`(Y(t[n]), `*`((D(Y))(t[n]), `*`(h)), `*`(`*`(`/`(1, 2), `*`(((`@@`(D, 2))(Y))(t[n]))), `*`(`^`(h, 2))))+O(`^`(h, 3)),h,3))](images/RKMethods_156.gif) |
(29) |
by dropping the term
It follows that this method has local truncation error
which leads to a global truncation error
. Hence the method is order 2.
Test Order of Accuracy: RK2
The following procedure defines the method RK2 (otherwise known as the Modified Euler Method).
Compute the solution a number of times specified by
, dividing the stepsize
by 2 each time.
For each stepsize
, compare the computed result at
of the
component of
(which corresponds to the original
), with the value
.
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_205.gif) |
![10](images/RKMethods_206.gif) |
(30) |
> |
![Typesetting:-mrow(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_207.gif) |
![1.0, 2.0](images/RKMethods_208.gif) |
(31) |
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_209.gif) |
![.4](images/RKMethods_210.gif) |
(32) |
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_211.gif) |
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_219.gif) |
![0.3906250000e-3](images/RKMethods_220.gif) |
(33) |
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_221.gif) |
![2560](images/RKMethods_222.gif) |
(34) |
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_223.gif) |
![Vector[column](%id = 150401424)](images/RKMethods_224.gif) |
(35) |
Compute successive ratios
> |
![Typesetting:-mrow(Typesetting:-mo(](images/RKMethods_225.gif) |
The ratio is approximately
.
Conclusion: Since
with
we conclude that RK2 is of order 2.
A Runge-Kutta Method of Order 3
For each order greater than
there is more than one possible Runge-Kutta method of the specified order.
The following is a Runge-Kutta method of order 3.
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_239.gif) |
> |
![Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi(](images/RKMethods_240.gif)
![Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi(](images/RKMethods_241.gif) |
Note that for this particular choice of parameters, the value
does not appear in the final formula although it is used in defining
.
The parameters have been chosen to ensure that the local truncation error is
and therefore the global truncation error is
.
Test Order of Accuracy: RK3
The following procedure defines the method RK3.
Compute the solution a number of times specified by
, dividing the stepsize
by 2 each time.
For each stepsize
, compare the computed result at
of the
component of
(which corresponds to the original
), with the value
.
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_296.gif) |
![10](images/RKMethods_297.gif) |
(38) |
> |
![Typesetting:-mrow(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_298.gif) |
![1.0, 2.0](images/RKMethods_299.gif) |
(39) |
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_300.gif) |
![.4](images/RKMethods_301.gif) |
(40) |
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_302.gif) |
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_310.gif) |
![0.3906250000e-3](images/RKMethods_311.gif) |
(41) |
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_312.gif) |
![2560](images/RKMethods_313.gif) |
(42) |
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_314.gif) |
![Vector[column](%id = 151466412)](images/RKMethods_315.gif) |
(43) |
Compute successive ratios
> |
![Typesetting:-mrow(Typesetting:-mo(](images/RKMethods_316.gif) |
The ratio is approximately
.
(Note that with a higher order method, the computed results stop improving because we reach the maximum accuracy for the given precision of the floating point arithmetic.)
Conclusion: Since
with
we conclude that RK3 is of order 3.
A Runge-Kutta Method of Order 4
The following is the most common Runge-Kutta method of order 4, as stated in the CS 370 Course Notes.
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_325.gif) |
> |
![Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi(](images/RKMethods_326.gif)
![Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi(](images/RKMethods_327.gif) |
The parameters have been chosen to ensure that the local truncation error is
and therefore the global truncation error is
.
Test Order of Accuracy: RK4
The following procedure defines the method RK4.
Compute the solution a number of times specified by
, dividing the stepsize
by 2 each time.
For each stepsize
, compare the computed result at
of the
component of
(which corresponds to the original
), with the value
.
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_382.gif) |
![10](images/RKMethods_383.gif) |
(46) |
> |
![Typesetting:-mrow(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_384.gif) |
![1.0, 2.0](images/RKMethods_385.gif) |
(47) |
Start with a larger value of
because the errors get very small quickly with this
-order method.
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_388.gif) |
![4.0](images/RKMethods_389.gif) |
(48) |
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_390.gif) |
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_398.gif) |
![0.3906250000e-2](images/RKMethods_399.gif) |
(49) |
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_400.gif) |
![256](images/RKMethods_401.gif) |
(50) |
> |
![Typesetting:-mrow(Typesetting:-mi(](images/RKMethods_402.gif) |
![Vector[column](%id = 150401568)](images/RKMethods_403.gif) |
(51) |
Compute successive ratios
> |
![Typesetting:-mrow(Typesetting:-mo(](images/RKMethods_404.gif) |
The ratio is converging to approximately
.
(Note that to get a better measure of the ratios, one would need to compute with higher-precision floating point.)
Conclusion: Since
with
we conclude that RK4 is of order 4.