January 24, 2017 Class 5---Error Analysis for the Euler Method (II),
Introduction to Newton's 2nd Law
Activities:
Error estimates
Notes in error estimate may be more readable in the pdf file
~sg/euler_error_analysis.pdf
There is also a link on the class web page.
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 Euler 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.
You might like to plot the error divided by the step size to
see that it is not constant.
What are the main take aways? First, in each step the leading
error is order (step size)^2, but when you integrate over a
finite range, there are order 1/(step size) steps, so the leading
error is order (step size)^1, i.e., linear.
Second, in general there are higher order terms and it may not
be possible to negect them. Also, if you think you can estimate
the higher order contribution, you had better be very careful
about the leading order term.
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 will looked at the
code in some detail. One key issue of discussion was pointer
dereferencing as when a pointer variable x_ptr is used with the
* operator in front of it.
*x_ptr = 6.3
places the value 6.3 in the memory location that x_prt points to.
If desired, we can talk more about this later.
We started looking at energy conservation using awk and found
that it is not conserved. I made graphs that showed the
change in the energy is linear with time and that the slope
depends on the step size. With smaller step size, change in
energy is decreased.
At the end of class, students were challenged to write a
script to find the energy after 1 second as a function of
integration step size. We will look at ways to do that in
the next class.