Jan. 23, 2013 Class 5---Error Analysis for the Euler Method (II),
Introduction to Newton's 2nd Law
Activities:
We would like to automate the procedure of finding the final
value y(2) as a function of step size so that we can easily
plot the final result as we vary the step size.
The shell script ~sg/chap2/ans_vs_step.csh solves
this problem. Here is ans_vs_step.csh with extra comments:
#!/bin/csh # anything after the hash mark or sharp is a comment
# the first comment is special in that it tells the shell
# to use /bin/csh to interpret or run the commands in this file
# such a comment starts with #!
# the word foreach starts a loop. The loop ends with the "end" line
foreach step ( 0.5 0.2 .1 .05 .02 .01 .005 )
echo -n $step " " # -n tells echo not to end with a new line
euler $step |tail -1|awk '{print $2}'
end
The script above uses the tail command to look at the last
line of output. There is also a head command to look at the
beginning of a file or standard input. head -9 will look at the
first 9 lines. head +3 will look at all of the file starting
with line 3. Consult the man pages for head and tail
to learn more.
We plotted the result as a function of step size. It is
clear that as the step size approaches zero we approach the
analytic answer. Let's examine how we approach the
correct answer.
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 used a program polyfit
to fit the result as a function of step size with a second
order polynomial. That result agrees with the error calculation
above. It is also possible to plot the error
divided by dt and find that the coefficient of the second
term is negative, not positive as we would get from just
looking at the second term in the Taylor expansion.
Solving Newton's Law of Motion
We 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.
In our next class, we will talk about energy conservation.