33.4 Unit Tests

33.4.1 Runge Kutta FLASH Test for a 2D Elliptical Path

33.4.1.1 Derivation of the 2D elliptical path ODE

Consider the following ODE in 2D, operating on a 2D vector ${\bf r}$:

$$\dfrac{d{\bf r}}{dt} = \left(\begin{array}{cc} \lambda_1 & \omega \\ -\omega & \lambda_2 \end{array}\right) {\bf r} = {\bf A}{\bf r}.$$ (33.22)

This ODE describes a rotational ($\omega$) and sheer ( $\lambda_1,\lambda_2$) force on a particle located at ${\bf r}$. If $\lambda_1,\lambda_2=0$, then the ODE describes a particle revolving in a circle of radius ${\bf r}(0)$, the initial radius of the particle. The general solution to the above ODE is:

$${\bf r}(t) = e^{{\bf A}t}{\bf r}(0).$$ (33.23)

The matrix exponential is given by (using ${\bf A}={\bf U}{\bf D}{\bf U}^{-1}$ and $e^{\bf A}={\bf U}e^{\bf D}{\bf U}^{-1}$, where ${\bf D}$ is diagonal):

$$e^{{\bf A}t} = \dfrac{1}{\Delta} \left(\begin{array}{cc} (e^{{\bf A}t})_{11} & (... ...bf A}t})_{12} \\ (e^{{\bf A}t})_{21} & (e^{{\bf A}t})_{22} \end{array}\right),$$ (33.24)

where

$$(e^{{\bf A}t})_{11} = \dfrac{1}{2}e^{t(\lambda_1+\lambda_2+\Delta)/2}(\lambda_1-\lambda... ... + \dfrac{1}{2}e^{t(\lambda_1+\lambda_2-\Delta)/2}(-\lambda_1+\lambda_2+\Delta)$$ (33.25)

$$(e^{{\bf A}t})_{12} = \omega(e^{t(\lambda_1+\lambda_2+\Delta)/2}-e^{t(\lambda_1+\lambda_2-\Delta)/2}$$ (33.26)

$$(e^{{\bf A}t})_{21} = - (e^{{\bf A}t})_{12}$$ (33.27)

$$(e^{{\bf A}t})_{22} = \dfrac{1}{2}e^{t(\lambda_1+\lambda_2+\Delta)/2}(-\lambda_1+\lambd... ...) + \dfrac{1}{2}e^{t(\lambda_1+\lambda_2-\Delta)/2}(\lambda_1-\lambda_2+\Delta)$$ (33.28)

and

$$\Delta = \sqrt{(\lambda_1-\lambda_2)^2-4\omega^2}.$$ (33.29)

For general values of $\lambda_1,\lambda_2,\omega$, the path of the particle is a spiral around the origin. Closed paths are possible, if at certain times $t\neq 0$ we have $e^{{\bf A}t}={\bf I}$. A necessary condition is thus that the offdiagonal elements of $e^{{\bf A}t}$ have to become zero for certain $t\neq 0$. Since each of these elements has a factor of the structure $e^B-e^C$, then either both $B$ and $C$ are real and equal or $B$ and $C$ are imaginary with oposite sign. Both situations can only happen, if $\lambda_2=-\lambda_1$. The first case ($B,C$ real and equal) is only possible if $\Delta=0$ and the second case ($B,C$ imaginary and of opposite sign) only if $\Delta$ is imaginary. Hence, necessary conditions for a closed path are:

$$\lambda_2 = - \lambda_1$$ (33.30)

$$\omega^2 \geq \lambda_1^2.$$ (33.31)

Let us now introduce the parameter $\vert k\vert>1$ such that:

$$\omega = k\lambda_1.$$ (33.32)

Then our closed path ODE becomes:

$$\dfrac{d{\bf r}}{dt} = \lambda_1 \left(\begin{array}{cc} 1 & k \\ -k & -1 \end{array}\right) {\bf r}$$ (33.33)

and, since $\lambda_1$ is now just merely a scaling factor, we can further condense the closed path ODE to:

$$\dfrac{d{\bf r}}{dt} = \left(\begin{array}{cc} 1 & k \\ -k & -1 \end{array}\right) {\bf r}.$$ (33.34)

This is the 2D ellipse ODE we are going to use for our Runge Kutta test. Inserting $\lambda_1=1$, $\lambda_2=-1$ and $\omega=k$ into the above general solution, we obtain, after inserting the identity $e^{i\phi}=\cos(\phi)+i\sin(\phi)$, the following real path solution:

$$\left(\begin{array}{c} x(t) \\ y(t) \end{array}\right) = \left(\begin{array}{cc} \dfrac{\sin(\phi)}{\sqrt{k^2-1}} + \cos(\... ...\phi) \end{array}\right) \left(\begin{array}{c} x_0 \\ y_0 \end{array}\right),$$ (33.35)

where

$$\phi = t\sqrt{k^2-1}.$$ (33.36)

For $\phi=n\pi; n=1,3,5,\ldots$, the 2x2 matrix is equal to $-{\bf I}$, i.e. the particle made a half revolution. For $\phi=n\pi; n=2,4,6,\ldots$, the 2x2 matrix is equal to ${\bf I}$ and the particle made a complete revolution. The time $T$ it takes for the particle to make one complete revolution (its time period) is thus:

$$T = \dfrac{2\pi}{\sqrt{k^2-1}}.$$ (33.37)

The square of the distance $d^2=[x(t)]^2+[y(t)]^2$ the particle is found from the origin has a minimum and a maximum at the following angles:

$$\phi_{max} (k>1) = \arctan\left[\sqrt{\dfrac{k+1}{k-1}}\dfrac{x_0+y_0}{x_0-y_0}\right]$$ (33.38)

$$\phi_{min} (k>1) = -\arctan\left[\sqrt{\dfrac{k-1}{k+1}}\dfrac{x_0-y_0}{x_0+y_0}\right]$$ (33.39)

$$\phi_{max} (k<1) = - \phi_{min} (k>1)$$ (33.40)

$$\phi_{min} (k<1) = - \phi_{max} (k>1).$$ (33.41)

These angles (range $[-\pi/2,+\pi/2]$) give only the angles to the nearest maximum/minimum point from the original position $(x_0,y_0)$. The other two angles corresponding to the remaining maximum/minimum distance points are obtained by adding $\pi$. The corresponding times are:

$$t_{max} = \phi_{max}/\sqrt{k^2-1}$$ (33.42)

$$t_{min} = \phi_{min}/\sqrt{k^2-1}.$$ (33.43)

The corresponding minimum and maximum square distances are:

$$d^2_{min} = \dfrac{k(x_0^2+y_0^2)+2x_0y_0}{k+\mbox{sgn}[k]}$$ (33.44)

$$d^2_{max} = \dfrac{k(x_0^2+y_0^2)+2x_0y_0}{k-\mbox{sgn}[k]},$$ (33.45)

where the signum function sgn is:

$$k>0 \longrightarrow [k] = +1$$ (33.46)

$$k<0 \longrightarrow [k] = -1.$$ (33.47)

From these formulas we see that

• if $x_0=y_0$ and $k>1$ , then the first $t_{min}$ is equal to zero, i.e. the particle is at its minimum distance $\sqrt{2}\vert x_0\vert$ from the origin.
• if $x_0=y_0$ and $k<-1$ , then the first $t_{max}$ is equal to zero, i.e. the particle is at its maximum distance $\sqrt{2}\vert x_0\vert$ from the origin.
• if $x_0=-y_0$ and $k>1$, then the first $t_{max}$ is equal to zero, i.e. the particle is at its maximum distance $\sqrt{2}\vert x_0\vert$ from the origin.
• if $x_0=-y_0$ and $k<-1$, then the first $t_{min}$ is equal to zero, i.e. the particle is at its minimum distance $\sqrt{2}\vert x_0\vert$ from the origin.
• the aspect ratio $A$ of the elliptical path is equal to $A=\sqrt{(k+\mbox{sgn}[k])/(k-\mbox{sgn}[k])}$.

We next develop the approach and formulas to find the minimum/maximum distance from a general point in a 2D domain to a general ellipse inside this 2D domain.

33.4.1.2 Minimum and maximum distance from a point to a general ellipse in a 2D cartesian domain

We will first get the general parametric equation of an ellipse in 2D cartesian space. We start by placing a 2D ellipse with its center at the cartesian origin and with its minor/major axis aligned along the cartesian y/x axis. Let the semi-major axis have length $a$ and the semi-minor axis have length $b$. The parametric equation for this ellipse is:

$$\left(\begin{array}{c} x \\ y \end{array}\right) = \left(\begin{array}{c} a\sin\theta \\ b\cos\theta \end{array}\right),$$ (33.48)

where the angle parameter is in the range $0\leq\theta\leq 2\pi$. Next we clockwise rotate the entire ellipse such that the minor axis makes an angle $\alpha$ with the cartesian y-axis. The range of this rotation angle is $0\leq\alpha\leq \pi$. The parametric equation for the rotated ellipse becomes

$$\left(\begin{array}{c} x \\ y \end{array}\right) = {\bf R}\left(\begin{array}{c} a\sin\theta \\ b\cos\theta \end{ar... ...}\right) \left(\begin{array}{c} a\sin\theta \\ b\cos\theta \end{array}\right),$$ (33.49)

where ${\bf R}$ is the rotation matrix that sends every point to its new rotated point. Finally we add the center translation vector ${\bf T}=(T_x,T_y)$ to the equation for ellipses that are not located at the cartesian origin:

$$\left(\begin{array}{c} x \\ y \end{array}\right) = \left(\begin{array}{cc} \cos\alpha & \sin\alpha \\ -\sin\alpha &... ...heta \end{array}\right)+ \left(\begin{array}{c} T_x \\ T_y \end{array}\right).$$ (33.50)

This is the general parametric equation of any ellipse in 2D cartesian space. Consider a specific point $(x,y)$ on the ellipse. The tangent line at this ellipse point contains the tiny derivative vector

$${\bf v} = \left(\begin{array}{c} dx/d\theta \\ dy/d\theta \end{array}\righ... ...\sin\theta \\ -a\sin\alpha\cos\theta-b\cos\alpha\sin\theta \end{array}\right),$$ (33.51)

with its origin at the ellipse point $(x,y)$. The other vector we consider is the vector

$${\bf w} = \left(\begin{array}{c} x_p - x \\ y_p - y \end{array}\right) = \... ... y_p + a\sin\alpha\sin\theta - b\cos\alpha\cos\theta - T_y \end{array}\right),$$ (33.52)

which also has its origin at the ellipse point $(x,y)$. The minimum distance occurs for those $\theta$ for which ${\bf v}\cdot {\bf w}=0$ (perpendicular vectors). This gives the following equation for $\theta$:

$$R\cos\theta + S\sin\theta + U\cos\theta\sin\theta = 0,$$ (33.53)

where

$$R = a[(y_p-T_y)\sin\alpha - (x_p-T_x)\cos\alpha]$$ (33.54)

$$S = b[(x_p-T_x)\sin\alpha + (y_p-T_y)\cos\alpha]$$ (33.55)

$$U = a^2-b^2.$$ (33.56)

Note, that $U\neq 0$, whenever $a\neq b$. For $a=b$ (a circle) we will have $U=0$ and the equation becomes easily solvable. The elliptical $U$ term makes the equation much harder to solve. Eliminating the $\sin\theta$ term, we arrive at the following quartic equation in $\cos\theta$:

$$\cos^4\theta + 2\left(\dfrac{S}{U}\right)\cos^3\theta + \left[\le... ...os^2\theta - 2\left(\dfrac{S}{U}\right)\cos\theta - \left(\dfrac{S}{U}\right)^2 = 0.$$

Its solutions give us up to 4 different possible angles $\theta$. However, to each $\cos\theta$ solution there correspond two $\sin\theta$ solutions via $\sin\theta = \pm \sqrt{1-\cos^2\theta}$, increasing the total possible solutions to 8. Each of these 8 solutions must be checked for the minimum and maximum distance given by $\vert{\bf w}\vert$. As far as the number of real solutions for the above quartic is concerned, when the point is inside or close to the ellipse, we can visualize that there are 4 possible directions from the point to the ellipse for which we have ${\bf v}\cdot {\bf w}=0$ (i.e. we hit the ellipse curve at right angles). Hence for this point location we expect 4 real quartic solutions. On the other hand, if we are far away from the ellipse, only 2 such points on the elliptical curve will exist and we expect 2 real and 2 complex solutions. We should never get 4 complex solutions of the above quartic.

33.4.1.3 The FLASH test problem for the Runge Kutta 2D ellipse test

We will now set up the unit test for the Runge Kutta 2D ellipse test. We start by stating the initial conditions and the requested shape of the ellipse:

1. The particle is initially at position $(x_0,y_0) = (1,1)$.
2. We want the particle to follow an elliptical path with given aspect ratio $A$.

From this requested aspect ratio we calculate our needed $k$ as:

$$k = \dfrac{A^2+1}{A^2-1}$$ (33.57)

and we take $k$ to be positive from now on. This means that the particle will start moving along the elliptical curve in the clockwise direction. By looking at the general formulas derived in the previous sections, we can state the following properties and equations for the ellipse. The square distance extrema and locations will be:

$$d^2_{min} = 2 \mbox{located at $(1,1)$ and $(-1,-1)$}$$ (33.58)

$$d^2_{max} = 2A^2 \mbox{located at $(A,-A)$ and $(-A,A)$} .$$ (33.59)

From this we see that the major and minor semi-axis will have lengths:

$$a = \sqrt{2}A$$ (33.60)

$$b = \sqrt{2}.$$ (33.61)

From the particle's initial position at $(1,1)$, the minor and major semi-axis will be reached at the following times ( $n=0,1,2,\ldots$):

$$t_{minor} = \dfrac{(A^2-1)}{2A}\pi n$$ (33.62)

$$t_{major} = \dfrac{(A^2-1)}{2A}\pi (n+1/2).$$ (33.63)

Hence the period of one complete revolution (returning to $(1,1)$) is given by:

$$T = \dfrac{(A^2-1)}{A}\pi.$$ (33.64)

The ellipse is represented by the following two (equivalent) forms:

$$\dfrac{(x-y)^2}{4A^2}+\dfrac{(x+y)^2}{4} = 1$$ (33.65)

and the parametric time equations:

$$x(t) = A\sin\left(\dfrac{2A}{A^2-1}t\right) + \cos\left(\dfrac{2A}{A^2-1}t\right)$$ (33.66)

$$y(t) = - A\sin\left(\dfrac{2A}{A^2-1}t\right) + \cos\left(\dfrac{2A}{A^2-1}t\right).$$ (33.67)

Each of the following Runge Kutta runs follows the elliptical curve through a full ellipse revolution, starting at the initial at position $(x_0,y_0) = (1,1)$. The calculations are stopped as soon as the $y=1$ line is crossed once, with the last crossing point still included in the calculation. Chosing an ellipse with an aspect ration of $A = 2$, we report the following data using fractional errors of $e_{frac}=10^{-8}$ and $10^{-12}$ and a $(x,y)$ coordinate error base vector of ${\bf e}_{base}=(1,1)$ (Eq.(33.15) for each Runge Kutta run in separate tables:

• The embedded RK scheme used and its order.
• The number of RK steps used to trace the full elliptical curve.
• The overall largest RK step coordinate error $max({\bf e}_{step})$, which should be $\leq e_{frac}\cdot{\bf e}_{base}$.
• The maximum time distance error $\max (e_{time})$ from the elliptical curve, where
  

  $$e_{time} = \sqrt{[x-x(t)]^2+[y-y(t)]^2}$$ (33.68)

  
  is a measure of how far the RK curve deviates from the true elliptical curve at time $t$. The quantity $e_{time}$ can be viewed as a global error of the RK scheme.
• The maximum of the closest distance error $\max (e_{closest})$ to the true elliptical curve. The closest elliptical distance error $e_{closest}$ is calculated using the approach outlined in section 33.4.1.2 and is a measure of how closely the RK scheme traces out the true elliptical curve.

Table 33.1 shows the obtained FLASH results. The higher order Fehlberg (5) and Cash-Karp (5) RK schemes outperform the lower order ones. The Cash-Karp RK scheme gives the best performance and is thus the set default RK scheme for the Runge Kutta unit.

Table 33.1: Number of RK steps, largest RK step coordinate error $max({\bf e}_{step})$, maximum time distance error $\max (e_{time})$ and maximum closest distance error $\max (e_{closest})$ values for a complete RK elliptical curve tracing with amplification $A = 2$, starting at the point $(x_0,y_0) = (1,1)$. Results shown as a function of embedded RK scheme used with fractional errors $e_{frac}=10^{-8}$ and $e_{frac}=10^{-12}$.

RK scheme (order)	RK steps	$max({\bf e}_{step})$	$\max (e_{time})$	$\max (e_{closest})$
Fractional error $e_{frac}=10^{-8}$
Euler-Heun (2)	67887	$8.10\times 10^{-9}$	$2.63\times 10^{-8}$	$1.55\times 10^{-11}$
Bogacki-Shampine (3)	1102	$7.39\times 10^{-9}$	$1.12\times 10^{-7}$	$1.09\times 10^{-7}$
Fehlberg (4)	228	$6.97\times 10^{-9}$	$1.78\times 10^{-6}$	$1.73\times 10^{-6}$
Fehlberg (5)	84	$6.84\times 10^{-9}$	$3.23\times 10^{-8}$	$3.02\times 10^{-8}$
Cash-Karp (5)	60	$7.14\times 10^{-9}$	$2.72\times 10^{-8}$	$2.61\times 10^{-8}$
Fractional error $e_{frac}=10^{-12}$
Euler-Heun (2)	6788695	$8.10\times 10^{-13}$	$3.85\times 10^{-12}$	$4.14\times 10^{-10}$
Bogacki-Shampine (3)	23729	$7.29\times 10^{-13}$	$1.12\times 10^{-11}$	$1.09\times 10^{-11}$
Fehlberg (4)	2272	$6.60\times 10^{-13}$	$1.75\times 10^{-9}$	$1.72\times 10^{-9}$
Fehlberg (5)	526	$6.06\times 10^{-13}$	$3.22\times 10^{-12}$	$3.10\times 10^{-12}$
Cash-Karp (5)	372	$6.13\times 10^{-13}$	$2.79\times 10^{-12}$	$2.69\times 10^{-12}$

33.4.2 Runge Kutta FLASH Test for Binomial Homogeneous Higher Order ODE

A binomial homogeneous ODE is an ODE of the form ($n\geq 1$):

$$\sum_{i=0}^{n}\binom{n}{i}f^{(i)}(x) = 0,$$ (33.69)

whose general solution is

$$f(x) = e^{-x}\left(\sum_{k=0}^{n-1}c_k x^k\right).$$ (33.70)

We are interested in the particular solution with $\forall c_k=1$:

$$f(x) = e^{-x}P(x)$$ (33.71)

$$P(x) = \sum_{k=0}^{n-1}x^k.$$ (33.72)

The derivatives of this particular solution are given by:

$$f^{(i)}(x) = e^{-x}\sum_{k=0}^i(-1)^{i+k}\binom{i}{k}P^{(k)}(x).$$ (33.73)

The values of the $k$-th derivative of $P(x)$ at $x = 0$ are the factorials of $k$, i.e. we have $P^{(k)}(0)=k!$. Thus the inital set of conditions at $x = 0$ for the particular solution are:

$$f^{(i)}(0) = i!\sum_{k=0}^i\frac{(-1)^k}{k!}$$

$$= !i,\;\;\;0\leq i\leq n-1$$ (33.74)

where $!i$ denotes the number of out-of-place permutations (derangements, no numbers in original places) of $i$. A recursive formula for $!i$ can be given: $!i=(i-1)[!(i-1)+!(i-2)]$, with starting values $!0=1$ and $!1=0$. Note that the given set of initial condition implies

$$f^{(n)}(0) = !n\;-\;n!,$$ (33.75)

due to the definition of the binomial ODE equation.

To solve the binomial ODE using the Runge Kutta scheme, the usual transformation to a set of coupled 1st order ODE's is done. Define the $n$-dimensional vector

$${\bf y} = \left(f,f^{(1)},f^{(2)},\ldots,f^{(n-1)}\right)^T,$$ (33.76)

the coupled set of 1st order ODE's are in matrix form:

$$\frac{d{\bf y}}{dx} = \left( \begin{array}{cccccc} 0 & 1 & 0 & 0 & \cdots & 0 \\ 0... ...n}{2} & -\binom{n}{3} & \cdots & -n \end{array} \right){\bf y}.$$ (33.77)

We compare the Runge Kutta vector ${\bf y}$ with the vector ${\bf y}$ evaluated using Eq.(33.73).