November 6, 2007 --- Class 20 --- Simpson's Rule Rederived, AHO, NDSolve,
Hit or Miss Method
Simpson's Rule Rederived
Starting from the trapezoid rule, we used the formulae below
for correcting the leading 1/N^2 errors.
S_N = A + E/N^2
S_2N = A + E/(2N)^2
so A = 4/3 S_2N -1/3 S_N
Applying this to a one-dimensional integral, we will easily find that
the weighting of the sum of two expressions on the right hand side
exactly coincides with the Simpson's rule formula. We previously
found that Simpson's Rule has a leading error of order 1/N^4. Since
the approach above removes the leading 1/N^2 error and since the
error of the trapezoid rule only has even powers of 1/N (asserted, but
not proved), this is not surprising.
Anharmonic Oscillator: Analytic Considerations
For homework you were asked to study the period of the
harmonic oscillator and how it depends on the force constant.
It turns out that a lot can be learned about this problem by
analytic means. We discussed how to use time translation to
get simple initial conditions. We then considered how we
could rescale x and t to simplify the equation of motion.
When we had done so, we found that we could determine the
dependence of the period on k, m and the initial amplitude,
up to a single unknown constant. This constant was the period
with unit mass, k and amplitude. We just found that period using
Mathematica.
Solving Differential Equations with Mathematica: Anharmonic Oscillator
If you have a simple differential equation to solve, you might want
to use Mathematica rather than a numerical code you develop yourself.
To solve a differential equation numerically, we use NDSolve. To solve
one analytically, use DSolve. For example:
NDSolve[ {x''[t] == - x[t]^3, x[0] == 1, x'[0] == 0}, x, {t, 0, 10}]
is the command required to solve the anharmonic oscillator with k/m=1.
x''[t] stands for the second derivative at time t. x'[0] stands for
the first derivative, or velocity, at t=0. Within the list that is
the first argument to NDSolve, we have the differential equation
itself and the initial conditions on position and velocity that are
required to specify the solution. The solution will be an
"interpolating function" for x[t] over the range
0 <= t <= 10.
If your next command is
Plot[ Evaluate[x[t] /. %], {t, 0, 10}]
Mathematica will display a plot of the solution. The command
Evaluate[x[t] /. %] tells Mathematica to use the output of the
previous command as a substitution rule and then to evaluate the
function x. Assuming that In[1] was the command to solve the
differntial equation, Out[1] is the rule for the interpolating function
that defines the solution. You can then find the solution at a
particular time with a command like:
x[7.5] /. Out[1]
I do not find the following commands intuitive, but they can be used
to find the period:
x /. First[Out[1]]
FindRoot[ %[t] == 1, {t, 7.5}]
The output of the first command is the interpolating function itself.
In the second command, % stands for the interpolating function, we
make it a function of t and then ask Mathematica to find the value of
t that makes the function 1. We have the search start near 7.5 because
looking at the plot of the solution to the anharmonic oscillator, we
see the the period is close to 7.5. The result from Mathematica is
7.4163.
Higher Dimensional Integrals
Monte Carlo methods like hit or miss and the sample mean methods
have errors that fall like n^{-1/2}. This is not very good. To
cut the error in half you need to do 4 times a much work. For
Simpson's rule the error goes like n^{-4}, so if you do twice a
much work, the error falls by a factor of 16. However, for
multidimensional integrals, a method with an error like h^a,
where h is the spacing between grid points, will have an error like
n^{-a/d} because the number of points on the grid goes like h^{-d}.
As an example, if a=4 and d=10, a/d=0.4 so the error is falling
more slowly than n^{-1/2}. So, for high dimension integrals,
Monte Carlo methods may be bad, but other methods are worse.
Hit-or-Miss Method
The Hit-or-Miss method and Sample Mean methods are described in CSM
Sec. 11.2. I introduced the hit-or-miss method. In the next class,
we will test this method.