Error analysis of the Euler algorithm: Local error proportional
to h^2, global error is proportional to h.
Wikipedia: Euler's Algorithm
Euler's original paper on his algorithm
- Runge-Kutta (4) algorithm: Local error is proportional
to h^5, global error is proportional to h^4.
- Dormand-Prince algorithm:
Read the entry in PHASER Help and also see
- How to compute e: comparison of several algorithms.
- What to do in practice: In real life applications, it is impossible
to be certain if numerically generated solutions of differential equations
are reliable, or how large the errors might be. The golden rule
of practice is
to compute the same solution using several algorithms with various
step sizes. If the variations in these computations are large (more than
you can accept for a particular application) your problem is most likely
a dangerous one. If the variations in these computations are small
enough, most likely you get a reliable answer out of the machine.
- Largest safe step size:
Consider the differential equation x'= -12x.
Show, using the Linearization Theorem, that x=0 is an
asymptotically stable equilibrium point.
Now, in the XiVsTime View of Phaser, compute several
solutions starting near the equilibrium point
using Dormand-Prince 5(4) with the default settings.
Assume that this is a correct picture (errors are very small).
Your picture should indeed show the asymptotic stability of the origin.
- Next, in the XiVsTime View of Phaser
algorithm with various step sizes compute the same solutions above. Use a large graph point size and connect points.
Determine the largest
step size that results in a picture where the origin
looks like an asymptotically stable point as you obtain with Dormand-Prince 5(4).
Hint: h = 0.3 does not give a correct picture.
- Impossible computation:
Consider the "explosion" problem x'= x^2, x(0)=1, and try to compute
x(1) using several different algorithms with various step sizes.
Verify by differentiation that the unique
solution of this initial value problem is x(t) = 1/(1-t).
Conclude that x(1) is infinite.
Now in Phaser, compute this number, x(1),
using steps h = 0.01 and h = 0.005 and algorithms Euler,
Improved Euler, Runge-Kutta4,
Dormand-Prince5(4), and Dormand-Prince8(5,3).
Following "What to do in practice" rule above, conclude that this problem is
indeed dangerous to compute.
Next, In the Dormand-Prince 5(4), set the step h = 0.01 and decrease
the Absolute error and Relative error from 1.0E-7 until Phaser
refuses to compute the value x(1). What is the smallest error tolerance
for which PHASER refuses to compute infinity? Monitoring the
Console in the Numericas Editor, report at which t value
PHASER stops computations for your set tolerances.
- Two steps of Improved Euler(2):
Consider the initial value problem x1' = x, x(0)=1.
With step size h = 0.5, compute, by hand using paper and pencil,
two steps of Improved Euler (2) algorithm to
obtain the approximate value of x(1). You can look up the formula for Improved Euler (2) in Phaser Help;
make sure you show all the intermediate numbers Ks.
Now compare your answer with the one you get from Phaser.
When the right-hand-side of a differential equation
contains time t explicitly, the equation is called nonautonomous.
Now, consider a differential equation of the form
dx/dt = f (t, x) with x(t_0) = x_0.
Then Euler becomes
x_(n+1) = x_n + h*f(x_n, t_n)
t_(n+1) = t_n + h .
Now, in the Gompertz equation
x' = a*(exp(-b*t))*x,
set a = 3 and b = 2.1, and x(0) = 4.5; and take step size h = 0.2.
By hand (with the help of a calculator)
compute two steps of Nonutonomous Euler on Gompertz and compare your numbers
to the ones you get from PHASER.
- Reading Euler: Try to read the original paper of Euler listed above.
Write a short commentary on this paper. Is his algorithm the same as the one
we derived in class? What does he have to say about errors?