- Logistic diferential equation
- Interacting populations
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.
Bacterial growth experiment:
A culture of 8 bacteria are placed in a test tube, and their
numbers is counted daily. At first, when their numbers was small
they grew at a rate of %185 per day. After a few days, their numbers
stabilized at about 410.
We will assume that the bacteria grow according
to the logistic model x'=ax - bx^2.
Determine the numerical values of the parameters a and b for our bacterial growth.
Hint: It is reasonable to assume that when their numbers is small, bacteria can
grow "exponentially;" that is a = 1.85 is a good approximation. Now using the
information about the equilibria, find b.
With the parameter values from above, using Phaser, determine how long it takes for the bacteria culture
to grow to within %10 of the carrying capacity. The default algorithm Dormand-Prince 5(4) is fine.
- Predator-prey model with logistic growth: Recall that in class we considered
the Predator-prey system ODE from Phaser
x1' = x1(a - bx2 - mx1)
x2' = x2(-c + dx1 -nx2)
Set a=2, b=1.1, c=1.5, d=0.5, m=0, and n=0. With these parameter values we assume that
in the absence of the other, each species grows or dies exponentially.
Take several (at least two) initial conditions and describe the behavior of the
populations in biological terms.
Now, increase the parameter m from 0 to 2.5 while keeping the other parameters fixed.
That is, assume that the prey grows according to the logistic model.
Describe the behavior of the populations as the parameter m is varied.
- Error vs. step size in Euler:
We mentioned in class that the global error bound in Euler's algorithm is
proportional to the step size. Now,
in PHASER solve
the initial-value problem x' = x, x(0) = 1 to compute x(1) = e = 2.7182818284590452354,
using Euler's algorithm with six different step sizes.
Using your favorite plotting program, plot the errors against the step sizes.
Do you get a linear relationship? Be careful of the scales on your graph.
- 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.
- Reading Euler: Try to read the original paper of Euler listed above.