September 11, 2007 --- Class 5---Error Analysis for the Euler Method (II),
Introduction to Chapter 3, Energy (non)conservation
Information:
Activities:
Error estimates
Most of of this class was devoted to understanding the
error in the Euler algorithm for integrating a differential
equation. We did a careful estimate of the error for the
simple case
dy/dx = f(x).
This is simple because the right hand side depends on x,
not y(x). So this is just doing an integral.
If we use a step size dx, the error in a bin with left hand
edge at x is:
(1/2!) (dx)^2 f'(x) + (1/3!) (dx)^3 f''(x) + higher order terms
In our first example of the Euler program, f(x)=2x. Thus,
f'(x)=2 and f''(x)=0, so the second term above is 0. Why
is the error order dx, not (dx)^2? That is because the
above expression is for each bin. In our example, we
integrate x from 1 to 2, so the number of bins n=(1/dx).
Thus, the total error is dx for this case.
In our second example, f(x)=3x^2. f'(x)=6x and f''(x)=6,
so the second error term does not vanish. If you very
carefully evalute the two error terms you will find that
the first term is
3 (dx)^2 * n (1 + (n-1)*dx/2 )
where n=1/dx. Using the fact that n*dx=1, this first
term can be written as
3 dx (3/2 -dx/2) or 9/2 *dx - 3/2 * (dx)^2.
Actually, only this term requires care. The f'' term is
easily seen to be
(dx)^2.
Adding the two terms, our final result is
9/2 *dx - 1/2 * (dx)^2.
Before we did the careful estimate, we plotted the error
divided by dt and found that the coefficient of the second
term is negative, not positive and we would get from just
looking at the second term in the Taylor expansion.
Introduction to Chapter 3
Solving Newton's Law of Motion
We started Chapter 3 and discussed how to turn a 2nd order
differential equation into a series of coupled first order
differential equations. We also discussed why position,
velocity and acceleration are so important, and the next
derivative, called the jerk, is not well known.
The program fall_euler_ansi.c in ~sg contains code to
integrate the equations for a falling body.
It is set up for the Euler algorithm. We discussed the
code in detail.
If the lines for updating v and y are
reversed, you get the Euler-Cromer algorithm.
You will apply the Euler-Richardson algorithm, and
explore the step size dependence of the trajectory of a
falling body. We talked about energy conservation and
how to write an awk expression to calculate the energy.
Energy (non)conservation
We studied the energy of a falling body as caculated by the
Euler algorithm. We looked at the energy at t = 1 sec.
We also created a shell script to automate the procedure:
#!/bin/csh
# this is the best shell script I have ever written
foreach step ( 0.025 0.05 0.1 0.2)
fall_euler_ansi <