| ||||||||||||||||||||||||||||||||||||||||||||||
| ||||||||||||||||||||||||||||||||||||||||||||||
The first few vectors $ c_{k}$, $d_{k}$ are given by
$
\begin{eqnarray}
c_{2}=1 \;\;\;\;\; && \;\; d_{2}=1 \\
&& \\
c_{3}=\left(\begin{array}{r} 3/2 \\ -1/2 \end{array} \right)
\;\;\; && \;\;
d_{3}=\left(\begin{array}{r} 4/3 \\-1/3 \end{array} \right) \\
&& \\
c_{4}=\left(\begin{array}{r} 23/12 \\-16/12 \\5/12
\end{array} \right) \; && \;\;
d_{4}=\left(\begin{array}{r} 19/12 \\-10/12 \\3/12
\end{array} \right) \\
&& \\
c_{5}=\left(\begin{array}{r} 55/24 \\-59/24 \\37/24 \\-9/24
\end{array}\right) \; && \;\;
d_{5}=\left(\begin{array}{r} 323/180 \\-264/180 \\159/180 \\-38/180
\end{array} \right)
\end{eqnarray}
$
Evaluation step: Insert the preliminary result $y_{n+1}^{P},\dot{y}_{n+1}^{P}$ in the physical law for $b[y,\dot{y}]$:
$
\begin{eqnarray}
b_{n+1}^{P} &\equiv& b\left[ y_{n+1}^{P},\dot{y}_{n+1}^{P}\right]
\end{eqnarray}
$
Corrector step:
$
\begin{eqnarray}
\dot{y}_{n+1}\Delta t - \dot{y}_{n}\Delta t& = &(\Delta t)^{2}
\left[ b_{n+1}^{P} - \frac{\textstyle 1}{\textstyle 2} \nabla b_{n+1}
- \frac{\textstyle 1}{ \textstyle12} \nabla^{2} b_{n+1}
- \frac{\textstyle 1}{\textstyle 24} \nabla^{3} b_{n+1}
- \dots \right]
\\
&& \\
y_{n+1} - y_{n} - \dot{y}_{n}\Delta t & = & \frac{\textstyle (\Delta t)^{2}}{\textstyle 2}
\left[ b_{n+1}^{P} - \frac{\textstyle 2}{\textstyle 3} \nabla b_{n+1}
- \frac{\textstyle 1}{\textstyle 12} \nabla^{2} b_{n+1}
- \frac{\textstyle 7}{\textstyle 180}\nabla^{3} b_{n+1}
- \dots \right]
\end{eqnarray}
$
Defining the vector
$
\begin{eqnarray}
b_{k}^{P} &\equiv& \{ b_{n+1}^{P},b_{n},\dots b_{n-k+3} \}^{T}
\end{eqnarray}
$
and coefficient vectors $ e_{k}$, $f_{k}$, we write
The first few coefficient vectors are
$
\begin{eqnarray}
e_{2}=1 \;\;\;\;\; && \;\; f_{2}=1 \\
&& \\
e_{3}=\left(\begin{array}{r} 1/2 \\1/2 \end{array}\right)
&& \;\;
f_{3}=\left(\begin{array}{r} 1/3 \\2/3 \end{array}\right) \\
&& \\
e_{4}=\left(\begin{array}{r} 5/12 \\8/12 \\-1/12
\end{array}\right) \; && \;\;
f_{4}=\left(\begin{array}{r} 3/12 \\10/12 \\-1/12
\end{array}\right) \\
&& \\
e_{5}=\left(\begin{array}{r} 9/24 \\19/24 \\-5/24 \\1/24
\end{array}\right) \; && \;\;
f_{5}=\left(\begin{array}{r} 38/180 \\171/180 \\-36/180 \\7/180
\end{array}\right)
\end{eqnarray}
$
4.3.3 Nordsieck Formulation of the PC MethodInstead of threading a NGB polynomial through preceding points, expand $y(t)$ about $t_{n}$. The resulting "Taylor predictor" may be written as
$
\begin{eqnarray}
y_{n+1}^{P} & = & y_{n}
+ y\,'_{n}\Delta t
+ y\,''_{n} \frac{\textstyle (\Delta t)^{2}}{\textstyle 2!}
+ y\,'''_{n} \frac{\textstyle (\Delta t)^{3}}{\textstyle 3!}
+ O[(\Delta t)^{4}]
\\
{y\,'}^{P}_{n+1} \Delta t & = &
\;\;\;\;\;\;\;\;\;
y\,'_{n} \Delta t
+ y\,''_{n} (\Delta t)^{2}
+ y\,'''_{n} \frac{\textstyle (\Delta t)^{3}}{\textstyle 2!}
+ O[(\Delta t)^{4}]
\\
{y\,''}^{P}_{n+1} \frac{\textstyle (\Delta t)^{2}}{\textstyle 2!} & = &
\;\;\;\;\;\;\;\; \;\;\;\;\;\;\;\; \;\;\;\;\;
y\,''_{n} \frac{\textstyle (\Delta t)^{2}}{\textstyle 2!}
+ y\,'''_{n} \frac{\textstyle (\Delta t)^{3}}{\textstyle 2!}
+ O[(\Delta t)^{4}]
\\
{y\,'''}^{P}_{n+1} \frac{\textstyle (\Delta t)^{3}}{\textstyle 3!} & = &
\;\;\;\;\;\;\;\; \;\;\;\;\;\;\;\; \;\;\;\;\;\;\;\; \;\;\;\;\;\;\;\;
\;\;\;\;\;\;\;
y\,'''_{n}\frac{\textstyle (\Delta t)^{3}}{\textstyle 3!}
+ O[(\Delta t)^{4}]
\end{eqnarray}
$
Defining the vector
$
\begin{eqnarray}
z_{n} \equiv \left(
\begin{array}{c}
y_{n} \\
y\,'_{n} \Delta t \\
y\,''_{n} \frac{\textstyle (\Delta t)^{2}}{\textstyle 2!} \\
\vdots
\end{array} \right)
\end{eqnarray}
$
and the (Pascal triangle) matrix
$
\begin{eqnarray}
A & \equiv &
\left(
\begin{array}{ccccc}
1 & 1 & 1 & 1 & \dots \\
0 & 1 & 2 & 3 & \dots \\
0 & 0 & 1 & 3 & \dots \\
& \ddots & \ddots & 1 & \ddots \\
& & & & \ddots
\end{array}
\right)
\end{eqnarray}
$
we have
$
\begin{eqnarray}
z_{n+1}^{P} = A \cdot z_{n}
\end{eqnarray}
$
Now evaluate: $\Longrightarrow$ Evaluation step: Insert $z_{n+1}^{P}$ in the force law: $b_{n+1}^{P} \equiv b[y_{n+1}^{P}, {y\,'}^{P}_{n+1}]$
Corrector step:
$
\begin{eqnarray}
z_{n+1} = z_{n+1}^{P} + \gamma c
\end{eqnarray}
$
with an optimized coefficient vector $c$. The first few vectors are
$
\begin{eqnarray}
c =
\left( \begin{array}{c} 1/6 \\5/6 \\1 \\1/3 \end{array} \right) \,,\;\;\;
\left( \begin{array}{c} 19/120 \\3/4 \\1 \\1/2 \\1/12
\end{array} \right)\,,\;\;\;
\left( \begin{array}{c} 3/20 \\251/360 \\1 \\11/18 \\1/6 \\1/60
\end{array} \right)\,,\dots
\end{eqnarray}
$
Advantages of Nordsieck PC: - Self-starting - Adjustable time steps
Stability: Again, between explicit (bad) and implicit (good).
4.3.4 Runge-Kutta Method for 2nd order ODE
Advantages of RK: - Self-starting - Adjustable time steps
But:
Repeated evaluation of the acceleration $b(y)$
in a single time step
4.2.5 Symplectic Algorithms [For some background information on Lagrange's and Hamilton's formulation of mechanics, seehere .]Hamiltonian dynamics is an important field of application for ODE algorithms. One simple performance test is the conservation of total mechanical energy. Another quantity that should be conserved in Hamiltonian dynamics is the symplectic form. A number of algorithms have been devised which by construction conserve this quantity.
In a classical system with $M$ degrees of freedom, let
the complete set of (generalized) coordinates be denoted by
$q$, and the conjugate momenta by $p$. Then
$
\frac{\textstyle d q}{\textstyle dt} = \nabla_{p} H( q, p) \;\;\;\;\;\;\;
\frac{\textstyle d p}{\textstyle dt} = -\nabla_{q} H( q, p)
$
where $H(q, p)$ is the Hamiltonian. Rewrite this by defining the phase space vector $z \equiv \{ q, p \}$:
$
\frac{\textstyle d z}{\textstyle dt} = J \cdot \nabla_{z} H(z)
$
where
$
J \equiv \left(\begin{array}{cc} 0 & I \\ - I & 0 \end{array} \right)
$
is the "symplectic matrix". ("symplectic" means "intertwined").
The canonical transformation
$
z(t_{0}) \; \Longrightarrow \; z(t)
$
conserves not only the energy (= numerical value of the Hamiltonian), but also the symplectic form
$
s(z_{1}, z_{2}) \equiv z_{1}^{T} \cdot J \cdot z_{2}
$
EXAMPLE: Harmonic oscillator:
$
H( z) \equiv H(q,p) = \frac{\textstyle k}{\textstyle 2} q^{2}
+ \frac{\textstyle p^{2}}{\textstyle 2m}
$
Canonical transformation:
$
z(t) = \left(\begin{array}{r} q \\ \\p \end{array} \right) =
\left( \begin{array}{cc}
\cos \omega t & \frac{\textstyle 1}{\textstyle m \omega } \sin \omega t \;\; \\
\\ -m \omega \sin \omega t&\cos \omega t \end{array} \right)
\cdot \left( \begin{array}{r} q(0) \\ \\ p(0) \end{array} \right)
\equiv A \cdot z(0)
$
Symplectic structure: Let $\{q_{1}(0), p_{1}(0) \}$ and $\{q_{2}(0), p_{2}(0) \}$ be two different initial conditions; then
$
\begin{eqnarray}
s(q_{1}(0) \dots p_{2}(0)) & \equiv & (q_{1}(0), p_{1}(0)) \cdot
\left(\begin{array}{cc} 0 & 1 \\ - 1 & 0 \\ \end{array} \right)
\cdot \left(\begin{array}{r} q_{2}(0) \\ p_{2}(0) \\ \end{array} \right) \\
&&
\\
& = & q_{1}(0) p_{2}(0) - p_{1}(0) q_{2}(0)
\end{eqnarray}
$
Geometric interpretation: $s$ is the area of a parallelogram defined by the two initial state vectors $ z_{1,2}$. $s$ is indeed constant:
$
\begin{eqnarray}
s(z_{1}(t),z_{2}(t)) & = & z_{1}^{T}(t) \cdot J \cdot z_{2}(t) \\
&& \\
& = & z_{1}^{T}(0) \cdot A^{T} \cdot J \cdot A \cdot z_{2}(0) \\
&& \\
& = & z_{1}^{T}(0) \cdot J \cdot z_{2}(0)
\end{eqnarray}
$
This means that $A^{T} \cdot J \cdot A = J$.
Let us check the performance of the Euler-Cauchy scheme:
$
z(t) = \left(\begin{array}{r} q \\ \\ p \end{array} \right) =
\left(\begin{array}{cc} 1 & \frac{\textstyle \Delta t}{\textstyle m} \; \\
\\-m \omega^{2} \Delta t & 1 \; \end{array} \right) \cdot
\left(\begin{array}{r} q(0) \\ \\ p(0) \end{array} \right)
\equiv E \cdot z(0)
$
The EC algorithm enhances the symplectic form by a factor $1+ (\omega \Delta t)^{2}$ at each time step. It is therefore not a symplectic algorithm.
Symplectic integrator of fourth order (Neri and Candy-Rozmus):
Third-order scheme (Ruth):
same structure as Candy's algorithm, but with coefficients
$
\begin{eqnarray}
(a_{1},a_{2},a_{3}) & = & (2/3,-2/3,1)
\\
(b_{1},b_{2},b_{3}) & = & (7/24, 3/4, -1/24)
\end{eqnarray}
$
Of the algorithms discussed previously only the Størmer-Verlet formula is symplectic.
A little theory:
Example: the first-order symplectic algorithm
$
q_{n+1}=q_{n}+P(p_{n})\Delta t
, \;\;\;\;
p_{n+1}=p_{n}+F(q_{n+1})\Delta t
$
exactly conserves
$
\tilde{H} \equiv H + H_{1}\Delta t + H_{2}(\Delta t)^{2}+
H_{3}(\Delta t)^{3}+ \dots
$
where
$
H_{1}=\frac{\textstyle 1}{\textstyle 2}H_{p}H_{q} , \;\;\;
H_{2}=\frac{\textstyle 1}{\textstyle 12}(H_{pp}H^{2}_{q}+H_{qq}H^{2}_{p}) , \;\;\;
H_{3}=\frac{\textstyle 1}{\textstyle 12}H_{pp}H_{qq}H_{p}H_{q}\;\dots
$
($H_{q} =\nabla_{q}H$ etc.) To take the harmonic oscillator, the perturbed Hamiltonian
$
\tilde{H}=H_{ho} + \frac{\textstyle \omega^{2} \Delta t}{\textstyle 2}pq
$
is conserved exactly.
A similar first order symplectic algorithm, known
as the Euler-Cromer formula,
$
p_{n+1}=p_{n}+F(q_{n})\Delta t , \;\;\;\;
q_{n+1}=q_{n}+P(p_{n+1})\Delta t
$
exactly conserves the perturbed Hamiltonian
$
\tilde{H}=H_{ho} - \frac{\textstyle \omega^{2} \Delta t}{\textstyle 2}pq
$
When applied to oscillator-like equations of motion it is a definite improvement over the (unstable) Euler-Cauchy method. EXERCISE: (See also here) Apply the (non-symplectic) RK method and the (symplectic) Størmer-Verlet algorithm (or the Candy procedure) to the one-body Kepler problem with elliptic orbit. Perform long runs to assess the long-time performance of the integrators. (For RK the orbit should eventually spiral down towards the central mass, while the symplectic procedures should only give rise to a gradual precession of the perihelion.) 4.2.6 Numerov's MethodMostly used for a certain kind of initial value problems (IVPs) that arise when boundary value problems (BVPs) are rewritten as IVPs. Many BVP have the general form
$
\frac{\textstyle d^{2}y}{\textstyle dx^{2}}=-g(x)y+s(x)
$
with given boundary values $y(x_{1})$ and $y(x_{2})$. Example: one-dimensional Poisson equation
$
\frac{\textstyle d^{2}\phi}{\textstyle dx^{2}}= -\rho(x)
$
with $\phi$ given at $x_{1}$ and $x_{2}$. Thus $g(x)=0$ and $s(x)=-\rho(x)$. Let us assume we have, instead of $y(x_{1})$ and $y(x_{2})$, the full set of initial values $y(x_{1})$ and $y'(x_{1})$.
Divide the interval $[x_{1},x_{2}]$ into sub-intervals of length $\Delta x$
and at each intermediate point $x_{n}$ expand $y(x)$
into a power series. Adding the Taylor formulae for $y_{n+1}$ and $y_{n-1}$
we find
$
y_{n+1}=2y_{n}-y_{n-1}+{y_{n}}''(\Delta x)^{2}+y_{n}^{(4)}
\frac{\textstyle (\Delta x)^{4}}{\textstyle 12}+O[(\Delta x)^{6}]
$
Inserting ${y_{n}}'' = -g_{n} y_{n}+s_{n}$ we have
$
y_{n+1}=2y_{n}-y_{n-1}+(\Delta x)^{2}[-g_{n}y_{n}+s_{n}]
+\frac{\textstyle (\Delta x)^{4}}{\textstyle 12}y_{n}^{(4)} +O[(\Delta x)^{6}]
$
The fourth derivative $y^{(4)}$ is approximated by
$
\begin{eqnarray}
y_{n}^{(4)} & = & \frac{\textstyle d^{2}y''}{\textstyle dx^{2}} |_{n} =
\frac{\textstyle d^{2}(-gy+s)}{\textstyle dx^{2}} |_{n}
\approx \frac{\textstyle 1}{\textstyle (\Delta x)^{2}}
\delta_{n}^{2}(-gy+s) = \\
&& \\
&=& \frac{\textstyle 1}{\textstyle (\Delta x)^{2}}
[-g_{n+1}y_{n+1}+2g_{n}y_{n}-g_{n-1}y_{n-1} + s_{n+1}-2s_{n}+s_{n-1}]
\end{eqnarray}
$
In this way we find Numerov's formula
$
\begin{eqnarray}
y_{n+1}[1+\frac{\textstyle (\Delta x)^{2}}{\textstyle 12}g_{n+1}] &=&
2y_{n}[1-\frac{\textstyle 5}{\textstyle 12}(\Delta x)^{2}g_{n}]-
y_{n-1}[1+\frac{\textstyle (\Delta x)^{2}}{\textstyle 12}g_{n-1}]+
\\
&& \\
&& +\frac{\textstyle (\Delta x)^{2}}{\textstyle 12}
[s_{n+1}+10s_{n}+s_{n-1}]+O[(\Delta x)^{6}]
\end{eqnarray}
$
EXERCISE: (See also here) Write a code that permits to solve a given second-order equation of motion by various algorithms. Apply the program to problems of point mechanics and explore the stabilities and accuracies of the diverse techniques. vesely 2005-10-10
|