Consistency and Stability Issues in the Numerical Integration of the First and Second Order Initial Value Problem

Author(s)
Isaac Fried

Abstract

In this note we consider some basic, yet unusual, issues pertaining to the accuracy and stability of numerical integration methods to follow the solution of first order and second order initial value problems (IVP). Included are remarks on multiple solutions, multi-step methods, effect of initial value perturbations, as well as slowing and advancing the computed motion in second order problems.

In this note we consider some basic, yet unusual, issues pertaining to the accuracy and stability of numerical integration methods to follow the solution of first order and second order initial value problems (IVP). Included are remarks on multiple solutions, multi-step methods, effect of initial value perturbations, as well as slowing and advancing the computed motion in second order problems.

Keywords

Initial Value Problems, Numerical Integration, Consistency, Stability, Multiple Solutions, Sensitivity to Initial Conditions, Slowing and Advancing the Computed Motion

Initial Value Problems, Numerical Integration, Consistency, Stability, Multiple Solutions, Sensitivity to Initial Conditions, Slowing and Advancing the Computed Motion

1. Introduction

Numerically following the solution of an evolution equation in time is one of the central tasks of numerical analysis and has received over the years vast and repeated attention, see [1] [2] . In this note we consider some essential and interesting, but possibly not widely known, aspects of the numerical solution of initial value problems. The rationale behind the procedures advocated here for the accurate and stable solution of the initial value problem, is often to greatly facilitate their introduction in the classroom. For the first order problem, we consider the usefulness of implicit methods, see [3] , for capturing multiple solutions emanating from a point of bifurcations, and also simple procedures for determining the accuracy and stability of multistep methods. For the second order problem, see [4] , we consider the means of slowing and advancing the computed speed of the numerical solution of the equation of motion.

2. The Generalized Mean Value Theorem—Taylor’s Theorem

The decisive significance of Taylor’s theorem (as can be looked up in any elementary calculus textbook) to applied mathematics in general, and to numerical analysis in particular, is that it ascertains that every differential function looks locally like a polynomial. Polynomials having the advantage of being easily computed, differentiated and integrated, relieving us by their use of the burden of possibly high complications in the symbolic manipulation of functions, often only implicitly given as the solution of an initial value problem (IVP) or a boundary value problem (BVP). We look at this theorem here from an unusual angle.

The mean value theorem (MVT) states that if function f(x), $f\left(0\right)=0$ , ${f}^{\prime}\left(0\right)\ne 0$ , is continuous in the closed interval [0, x] and differentiable on the open interval (0, x), then point ξ exists, strictly inside the interval, $0<\xi <x$ , at which the slope of the chord equals to the slope of the tangent line to f(x) at point ξ, or

$\frac{f\left(x\right)-f\left(0\right)}{x}={f}^{\prime}\left(\xi \right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{or}\text{\hspace{0.17em}}\text{\hspace{0.17em}}f\left(x\right)=x{f}^{\prime}\left(\xi \right),\text{\hspace{0.17em}}0<\xi <x$ (1)

implying that f(x) is nearly linear near x = 0.

The MVT theorem is a direct result, of the geometrically intuitively plausible, Rolle’s theorem. The generalized mean value theorem (GMVT) is a result of the application of the MVT (or Rolle’s theorem) to the higher order derivative functions of f(x). Here it is in its most concise form: Let function f(x) be continuous at x = 0, and such that

$f\left(0\right)=0,{f}^{\prime}\left(0\right)=0,{f}^{\u2033}\left(0\right)=0,\cdots ,{f}^{\left(n-1\right)}\left(0\right)=0,{f}^{\left(n\right)}\left(0\right)\ne 0.$ (2)

Then, function f(x) may be expressed in the form

$f\left(x\right)=\frac{1}{n!}{x}^{n}{f}^{\left(n\right)}\left(\xi \right),\text{\hspace{0.17em}}0<\xi <x$ (3)

implying that if f^{(n)}(x) is bounded near x = 0, then f(x), like x^{n}, is small if
$\left|x\right|\ll \text{1}$ .

3. Where Is ξ

We consider first the simplest case of n = 1, for which Equation (3) is

$f\left(x\right)=x{f}^{\prime}\left(\xi \right),\text{\hspace{0.17em}}f\left(0\right)=0,\text{\hspace{0.17em}}{f}^{\prime}\left(0\right)\ne 0,\text{\hspace{0.17em}}0<\xi <x$ (4)

and take

$f\left(x\right)=Ax+B{x}^{2}+C{x}^{3},\text{\hspace{0.17em}}{f}^{\prime}\left(x\right)=A+2Bx+3C{x}^{2}.$ (5)

such that $A={f}^{\prime}\left(0\right),B={f}^{\u2033}\left(0\right)/2!,C={f}^{\u2034}\left(0\right)/3!$ .

We further assume that, approximately

$\xi =kx+m{x}^{2}$ (6)

and have from this that

$$ (68)

For ${\alpha}_{0}=-3/4$ , ${\alpha}_{1}=7/4$ , and with ${y}^{\prime}\left(t\right)=0$ , the two-step method becomes

${y}_{2}=\frac{1}{4}\left(-3{y}_{0}+7{y}_{1}\right)$ (69)

for which the characteristic equation is

$4{z}^{2}-7z+3=0,\text{\hspace{0.17em}}{z}_{1}=1,\text{\hspace{0.17em}}{z}_{2}=\frac{3}{4}$ (70)

and hence

${y}_{n}={c}_{1}{\left(1\right)}^{n}+{c}_{2}{\left(\frac{3}{4}\right)}^{n}.$ (71)

Say we start the method with ${y}_{0}=\text{1}$ , ${y}_{1}=1+\epsilon $ , such that

${c}_{1}+{c}_{2}=1,\text{\hspace{0.17em}}\text{\hspace{0.17em}}{c}_{1}+\frac{3}{4}{c}_{2}=1+\epsilon $ (72)

resulting in ${c}_{1}=1+4\epsilon $ , ${c}_{2}=-4\epsilon $ , and

${y}_{n}=1+4\epsilon -4\epsilon {\left(\frac{3}{4}\right)}^{n}$ (73)

and

${y}_{\infty}=1+4\epsilon .$ (74)

17. A Three Step Method

The integration method

${y}_{3}={y}_{2}+\frac{1}{12}\tau \left(23{{y}^{\prime}}_{2}-16{{y}^{\prime}}_{1}+5{{y}^{\prime}}_{0}\right)$ (75)

is correct for y = 1, y = t, y = t^{2} and y = t^{3}. For

$y\left(t\right)=\frac{1}{24}{M}_{4}{t}^{4}$ (76)

we have from Equation (75) that

$y\left(3\tau \right)-{y}_{3}=\frac{3}{8}{M}_{4}{\tau}^{3}.$ (77)

For ${y}^{\prime}=-y$ , $y\left(0\right)={y}_{0}$ the characteristic equation of the three step method becomes

$12{z}^{3}+\left(-12+23\tau \right){z}^{2}-16\tau z+5\tau =0$ (78)

and at τ = 0 it reduces to

$12{z}^{3}-12{z}^{2}=0$ , of roots ${z}_{1}=1,{z}_{2}={z}_{3}=0$ . (79)

Implicit differentiation of Equation (78) with respect to τ yields

$36{z}^{2}{z}^{\prime}+23{z}^{2}+\left(-12+23\tau \right)2z{z}^{\prime}-16z-16\tau {z}^{\prime}+5=0$ (80)

and at τ = 0, z = 1, we have that ${z}^{\prime}=-1$ , so that near τ = 0

$z=1-\tau $ (81)

and the method is stable.

We further have from Equation (79) that at τ = 6/11, z_{3} = −1 and z_{1} = z_{2} = 1/2.

At the repeating root z = 0, the derivative function ${z}^{\prime}$ does not exist. Instead we write τ in terms of z as

$\tau =\frac{12\left(1-z\right)}{23{z}^{2}-16z+5}{z}^{2}$ (82)

and if z = 0, nearly, then

$\tau =\frac{12}{5}{z}^{2}$ (83)

and

$z=\pm \sqrt{\frac{5}{12}\tau}.$ (84)

18. Advancing and Retarding the Computed Motion

Next, we turn our attention to the second order IVP, see [6] and [7] , as typified by the model second order problem

${y}^{\u2033}+y=0,\text{\hspace{0.17em}}y\left(0\right)={y}_{0},\text{\hspace{0.17em}}{y}^{\prime}\left(0\right)={{y}^{\prime}}_{0}$ (85)

which we propose to approximate as

$\frac{{y}_{0}-2{y}_{1}+{y}_{2}}{{\tau}^{2}}+\left(1\text{\hspace{0.17em}}+\in \right){y}_{1}=0,\text{\hspace{0.17em}}\in \text{\hspace{0.17em}}=\alpha \tau +\beta {\tau}^{2}.$ (86)

The characteristic equation of this method is

${z}^{2}+\left(-2+{\tau}^{2}\left(1\text{\hspace{0.17em}}+\in \right)\right)z+1=0.$ (87)

For a sufficiently small τ, the roots of the characteristic equation are complex and $\left|z\right|=\text{1}$ . Hence the closed-form prediction of the computed y at step n

${y}_{n}={c}_{1}\mathrm{cos}\left(n\theta \right)+{c}_{2}\mathrm{sin}\left(n\theta \right)$ (88)

where c_{1} and c_{2} are determined by the initial conditions.

From the characteristic Equation (87) we have that

$\mathrm{cos}\left(\theta \right)=1-\frac{1}{2}{\tau}^{2}\left(1\text{\hspace{0.17em}}+\in \right),\text{\hspace{0.17em}}\mathrm{sin}\left(\theta \right)=\frac{1}{2}\sqrt{4-{\left(-2+{\tau}^{2}\left(1\text{\hspace{0.17em}}+\in \right)\right)}^{2}}$ (89)

or

$\theta =\tau +\frac{1}{2}\alpha {\tau}^{2}+\frac{1}{24}\left(1-3{\alpha}^{2}+12\beta \right){\tau}^{3}+O\left({\tau}^{5}\right),\text{\hspace{0.17em}}\in \text{\hspace{0.17em}}=2-{\tau}^{2}-2\mathrm{cos}\left(\tau \right).$ (90)

To drop the τ^{2} term in the above equation we take α = 0, and are left with

$\frac{\theta}{\tau}=1+\frac{1}{24}\left(1+12\beta \right){\tau}^{2}+O\left({\tau}^{4}\right)$ (91)

suggesting that θ may be advanced or retarded relative to τ with a proper choice of β. For instance, for β = −1/12 we have that θ = τ, nearly.

19. Integration of the Equation of Motion Linear Prediction

Inasmuch as the initial conditions to the second order equation of motion are usually given in terms of initial position and initial velocity, we prefer to reduce the second order problem into a coupled system of first order equations for position and velocity, and follow them in tandem.

To remain both concise and specific, we consider the model initial value problem

${x}^{\prime}=-y,\text{\hspace{0.17em}}{y}^{\prime}=x,\text{\hspace{0.17em}}{x}_{0}=x\left(0\right)=1,\text{\hspace{0.17em}}{y}_{0}=y\left(0\right)=0$ (92)

where x = x(t), y = y(t), t > 0, and where ${()}^{\prime}$ denotes differentiation with respect to time t. This initial value problem, that coincides with a single second order problem, is solved by $x=\mathrm{cos}\left(t\right)$ , $y=\text{sin}\left(t\right)$ representing a constant circular motion of period $T=2\text{\pi}$ .

We propose to follow the initial value problem with the explicit scheme

$\begin{array}{l}{x}_{1}={x}_{0}+\tau {{x}^{\prime}}_{0}\\ {y}_{1}={y}_{0}+\tau \left({\alpha}_{0}{{y}^{\prime}}_{0}+{\alpha}_{1}{{y}^{\prime}}_{1}\right)\end{array}$ (93)

in which τ is the time step, where x_{1} = x(τ) and y_{1} = y(τ), approximately, and with the coefficients α_{0}, α_{1} to be presently determined by stability and accuracy considerations.

With ${x}^{\prime}=-y$ , ${y}^{\prime}=x$ , system (93) becomes the system of recursions

$\begin{array}{l}{x}_{1}={x}_{0}-\tau {y}_{0}\\ {y}_{1}={y}_{0}+\tau \left({\alpha}_{0}{x}_{0}+{\alpha}_{1}{x}_{1}\right)\end{array}$ (94)

that explicitly produces x_{1} and y_{1} out of x_{0} and y_{0}, and then x_{2} and y_{2} out of x_{1} and y_{1}, and so on up to x_{n} and y_{n}. System (94) is solved by
${x}_{n}={z}^{n}{x}_{0}$ ,
${y}_{n}={z}^{n}{y}_{0}$ for magnification factor z that satisfies the pair of linear equations

$\begin{array}{l}z{x}_{0}={x}_{0}-\tau {y}_{0}\\ z{y}_{0}={y}_{0}+\tau \left({\alpha}_{0}{x}_{0}+{\alpha}_{1}z{x}_{0}\right)\end{array}$ (95)

for any x_{0} and y_{0}. Equation (95) is recast in the matrix-vector form to assume the form

$\left[\begin{array}{cc}z-1& \tau \\ -\tau \left({a}_{0}+{a}_{1}z\right)& z-1\end{array}\right][\begin{array}{c}{x}_{0}\\ {y}_{0}\end{array}]=0$ (96)

and the condition that it have a nontrivial solution is that the determinant of the system’s matrix of coefficients be zero, leading to the characteristic equation

$\mathrm{det}\left[\begin{array}{cc}z-1& \tau \\ -\tau \left({a}_{0}+{a}_{1}z\right)& z-1\end{array}\right]=0,\text{\hspace{0.17em}}{z}^{2}+2\left(-1+\frac{1}{2}{a}_{1}{\tau}^{2}\right)z+1+{a}_{0}{\tau}^{2}=0$ (97)

for z. The periodic nature of the solution to this initial value problems dictates that z be complex. Let
$\left|z\right|$ be the modulus of complex z. If
$\left|z\right|<\text{1}$ , then
${\left|z\right|}^{n}\to 0$ as
$n\to \infty $ , and if
$\left|z\right|>\text{1}$ , then
${\left|z\right|}^{n}\to \infty $ as
$n\to \infty $ . To avoid these undesirable eventualities of an artificial energy sink and an artificial, numerically induced, energy source, we select α_{0} = 0 in Equation (94), and are left with the reduced characteristic equation

${z}^{2}+2\left(-1+\frac{1}{2}{\alpha}_{1}{\tau}^{2}\right)z+1=0,\text{\hspace{0.17em}}a{z}^{2}+bz+c=0$ (98)

that possesses two complex roots z_{1} and z_{2} such that
${z}_{1}{z}_{2}=c/a={\left|z\right|}^{2}=1$ .

In fact,

$z=1-\frac{1}{2}{\alpha}_{1}{\tau}^{2}\pm i\tau \sqrt{{\alpha}_{1}-\frac{1}{4}{\alpha}_{1}^{2}{\tau}^{2}}$ (99)

where i^{2} = −1, and z is complex if

${\alpha}_{1}>0,\text{\hspace{0.17em}}4-{\alpha}_{1}{\tau}^{2}>0.$ (100)

By the fact that $\left|z\right|=\text{1}$ , the complex solution to Equation (98) is energy conserving, and may be written as

${z}_{1}=\mathrm{cos}\left(\theta \right)-i\mathrm{sin}\left(\theta \right),\text{\hspace{0.17em}}{z}_{2}=\mathrm{cos}\left(\theta \right)+i\mathrm{sin}\left(\theta \right)$ (101)

with

$\mathrm{cos}\left(\theta \right)=1-\frac{1}{2}{\alpha}_{1}{\tau}^{2},\text{\hspace{0.17em}}\mathrm{sin}\left(\theta \right)=\tau \sqrt{{\alpha}_{1}-\frac{1}{4}{\alpha}_{1}^{2}{\tau}^{2}}.$ (102)

Now, x_{n} and y_{n} are generally written as

${x}_{n}={c}_{1}{z}_{1}^{n}+{c}_{2}{z}_{2}^{n},\text{\hspace{0.17em}}{y}_{n}={{c}^{\prime}}_{1}{z}_{1}^{n}+{{c}^{\prime}}_{2}{z}_{2}^{n}$ (103)

with constants
${c}_{1},{c}_{2},{{c}^{\prime}}_{1},{{c}^{\prime}}_{2}$ determined by the initial conditions. Given x_{0} = 1, y_{0} = 0 we get from Equation (94), x_{1} = 1, y_{1} = α_{1}τ. Writing x_{n} and y_{n} in Equation (102) for n = 0 and n = 1 we obtain the two systems of linear equations

$\left[\begin{array}{cc}1& 1\\ {z}_{1}& {z}_{2}\end{array}\right][\begin{array}{c}{c}_{1}\\ {c}_{2}\end{array}]=[\begin{array}{c}1\\ 1\end{array}],\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left[\begin{array}{cc}1& 1\\ {z}_{1}& {z}_{2}\end{array}\right][\begin{array}{c}{{c}^{\prime}}_{1}\\ {{c}^{\prime}}_{2}\end{array}]={a}_{1}\tau [\begin{array}{c}0\\ 1\end{array}]$ (104)

readily solved for ${c}_{1},{c}_{2},{{c}^{\prime}}_{1},{{c}^{\prime}}_{2}$ as

$[\begin{array}{c}{c}_{1}\\ {c}_{2}\end{array}]=\frac{1}{{z}_{2}-{z}_{1}}[\begin{array}{c}{z}_{2}-1\\ -{z}_{1}+1\end{array}],\text{\hspace{0.17em}}\text{\hspace{0.17em}}[\begin{array}{c}{{c}^{\prime}}_{1}\\ {{c}^{\prime}}_{2}\end{array}]=\frac{{a}_{1}\tau}{{z}_{2}-{z}_{1}}[\begin{array}{c}-1\\ 1\end{array}]$ (105)

in which
${z}_{1}=\mathrm{cos}\left(\theta \right)-i\mathrm{sin}\left(\theta \right)$ ,
${z}_{2}=\mathrm{cos}\left(\theta \right)+i\mathrm{sin}\left(\theta \right)$ ,
${z}_{2}-{z}_{1}=2i\mathrm{sin}\left(\theta \right)$ . Writing z_{1} and z_{2} in terms of θ recasts Equation (103) into the form

$\begin{array}{l}{x}_{n}=\left({c}_{1}+{c}_{2}\right)\mathrm{cos}\left(n\theta \right)+i\left({c}_{2}-{c}_{1}\right)\mathrm{sin}\left(n\theta \right),\\ {y}_{n}=\left({{c}^{\prime}}_{1}+{{c}^{\prime}}_{2}\right)\mathrm{cos}\left(n\theta \right)+i\left({{c}^{\prime}}_{2}-{{c}^{\prime}}_{1}\right)\mathrm{sin}\left(n\theta \right)\end{array}$ (106)

and we have from Equation (105) that

$\begin{array}{l}{c}_{1}+{c}_{2}=1,\text{\hspace{0.17em}}{c}_{2}-{c}_{1}=-i{\left(\mathrm{sin}\left(\theta \right)\right)}^{-1}{\alpha}_{1}{\tau}^{2}\\ {{c}^{\prime}}_{1}+{{c}^{\prime}}_{2}=0,\text{\hspace{0.17em}}{{c}^{\prime}}_{2}-{{c}^{\prime}}_{1}=-i{\left(\mathrm{sin}\left(\theta \right)\right)}^{-1}{\alpha}_{1}\tau \end{array}$ (107)

with which we finally get

${x}_{n}=\mathrm{cos}\left(n\theta \right)+\frac{1}{2}{\alpha}_{1}{\tau}^{2}{\left(\mathrm{sin}\left(\theta \right)\right)}^{-1}\mathrm{sin}\left(n\theta \right),\text{\hspace{0.17em}}{y}_{n}={\alpha}_{1}\tau {\left(\mathrm{sin}\left(\theta \right)\right)}^{-1}\mathrm{sin}\left(n\theta \right)$ (108)

as the general numerical solution to our initial value problem.

20. Period Control

A cycle is completed when
$\mathrm{sin}\left(n\theta \right)=0$ or
$n\theta =2\text{\pi}$ . Then, according to Equation (108) y_{n} = 0 and x_{n} = 1. From
$\tau n\theta =2\text{\pi}\tau $ and T = nτ we obtain the computed period as

$T=2\text{\pi}\frac{\tau}{\theta}$ (109)

and to retain
$T=2\text{\pi}$ we select α_{1} in Equation (93) so as to guarantee τ = θ or
$\mathrm{sin}\left(\tau \right)=\mathrm{sin}\left(\theta \right)$ . This condition becomes, in view of Equation (102),

$\mathrm{sin}\left(\tau \right)=\tau \sqrt{{\alpha}_{1}-\frac{1}{4}{\alpha}_{1}^{2}{\tau}^{2}}$ (110)

leading to the quadratic equation

$\frac{1}{4}{\tau}^{2}{\alpha}_{1}^{2}-{\alpha}_{1}+\frac{{\mathrm{sin}}^{2}\left(\tau \right)}{{\tau}^{2}}=0$ (111)

for α_{1}, and resulting in

${\alpha}_{1}=2\left(1-\mathrm{cos}\left(\tau \right)\right)/{\tau}^{2}$ (112)

or

${\alpha}_{1}=1-\frac{1}{12}{\tau}^{2}+\frac{1}{360}{\tau}^{4}.$ (113)

if τ is small.

21. Quadratic Prediction

Inclusion of the acceleration in the prediction of x_{1} suggests the higher order scheme

$\begin{array}{l}{x}_{1}={x}_{0}+\tau {{x}^{\prime}}_{0}+\frac{1}{2}{\tau}^{2}{{x}^{\u2033}}_{0}\\ {y}_{1}={y}_{0}+\frac{1}{2}\tau \left({\alpha}_{0}{{y}^{\prime}}_{0}+{\alpha}_{1}{{y}^{\prime}}_{1}\right)\end{array}$ (114)

that becomes for

$\begin{array}{l}{x}^{\prime}=-y,\text{\hspace{0.17em}}{y}^{\prime}=x,\text{\hspace{0.17em}}{x}^{\u2033}=-x\\ {x}_{1}={x}_{0}-\tau {y}_{0}-\frac{1}{2}{\tau}^{2}{x}_{0},\\ {y}_{1}={y}_{0}+\frac{1}{2}\tau \left({\alpha}_{0}{x}_{0}+{\alpha}_{0}{x}_{1}\right).\end{array}$ (115)

Substitution of x_{1} = zx_{0}, y_{1} = zy_{0} in Equation (115) results in the system

$\left[\begin{array}{cc}z-1+\frac{1}{2}{\tau}^{2}& \tau \\ -\frac{1}{2}\tau \left({a}_{0}+{a}_{1}z\right)& z-1\end{array}\right][\begin{array}{c}{x}_{0}\\ {y}_{0}\end{array}]=0$ (116)

from which we obtain the quadratic characteristic equation

$\begin{array}{l}\mathrm{det}\left[\begin{array}{cc}z-1+\frac{1}{2}{\tau}^{2}& \tau \\ -\frac{1}{2}\tau \left({a}_{0}+{a}_{1}z\right)& z-1\end{array}\right]=0,\\ {z}^{2}+2\left(-1+\frac{1}{4}{\tau}^{2}+\frac{1}{4}{a}_{1}{\tau}^{2}\right)z+1+\frac{1}{2}{\tau}^{2}\left({a}_{0}-1\right)=0\end{array}$ (117)

for magnification factor z. To assure
$\left|z\right|=\text{1}$ for the complex roots of Equation (117) we set α_{0} = 1 and are left with

${z}^{2}+2\left(-1+\frac{1}{4}{\tau}^{2}\beta \right)z+1=0$ (118)

where β = 1 + α_{1}. The two roots of Equation (118) are

$z=1-\frac{1}{4}{\tau}^{2}\beta \pm i\tau \sqrt{\frac{1}{2}\beta -\frac{1}{16}{\tau}^{2}{\beta}^{2}}$ (119)

and z is complex if

$\beta >0,\text{\hspace{0.17em}}8-{\tau}^{2}\beta >0.$ (120)

Because $\left|z\right|=\text{1}$ we may write the complex roots of Equation (118) as

$z=\mathrm{cos}\left(\theta \right)\pm i\mathrm{sin}\left(\theta \right),\text{\hspace{0.17em}}\mathrm{cos}\left(\theta \right)=1-\frac{1}{4}{\tau}^{2}\beta ,\text{\hspace{0.17em}}\mathrm{sin}\left(\theta \right)=\tau \sqrt{\frac{1}{2}\beta -\frac{1}{16}{\tau}^{2}{\beta}^{2}}.$ (121)

The numerical scheme is period conserving if τ = θ, or $\mathrm{sin}\theta =\mathrm{sin}\tau $ . This is assured, according to Equation (121), by β such that

$\mathrm{sin}\left(\tau \right)=\tau \sqrt{\frac{1}{2}\beta -\frac{1}{16}{\tau}^{2}{\beta}^{2}}$ (122)

or

$\frac{1}{16}{\tau}^{2}\beta -\frac{1}{2}\beta +{\left(\frac{\mathrm{sin}\left(\tau \right)}{\tau}\right)}^{2}=0$ (123)

resulting in

$\beta =4\left(1-\mathrm{cos}\left(\tau \right)\right)/{\tau}^{2}$ (124)

or

${\alpha}_{1}=1-\frac{1}{6}{\tau}^{2}+\frac{1}{180}{\tau}^{4}$ (125)

if τ is small.

22. Conclusion

We accomplished here showing how to routinely determine the consistency and stability of any multistep method, explicit as well as implicit, for the stepwise integration of the first order initial value problem. We have also demonstrated here the advantage of using implicit methods to capture different solutions emanating from a branch-off point. For the integration of the second order equation of motion we have shown how to advance and retard the motion of the computed solution.

Cite this paper

Fried, I. (2019) Consistency and Stability Issues in the Numerical Integration of the First and Second Order Initial Value Problem.*Applied Mathematics*, **10**, 676-690. doi: 10.4236/am.2019.108048.

Fried, I. (2019) Consistency and Stability Issues in the Numerical Integration of the First and Second Order Initial Value Problem.

References

[1] Hairer, E., Norset, S.P. and Wanner, G. (1980) Solving Ordinary Differential Equations. Springer-Verlag, Berlin, Heidelberg.

[2] Gautschi, W. (2010) Numerical Analysis. 2nd Edition, Birkhauser, Berlin.

[3] Fried, I. (2000) The Advantage of an Implicit Correction at Bifurcation. Journal of Sound and Vibration, 235, 153-155.

https://doi.org/10.1006/jsvi.1999.2829

[4] Fried, I. and Waldrop, L. (2001) Kira Is Energy Consuming. Journal of Sound and Vibration, 244, 169-172.

https://doi.org/10.1006/jsvi.2000.3467

[5] Dahlquist, G. (1963) A Special Stability Problem for Linear Multistep Methods. BIT Numerical Mathematics, 3, 27-43.

https://doi.org/10.1007/BF01963532

[6] Lambert, J.D. and Watson, I.A. (1976) Symmetric Multistep Methods for Periodic Initial Value Problems. IMA Journal of Applied Mathematics, 18, 189-202.

https://doi.org/10.1093/imamat/18.2.189

[7] Quinlan, G.D. and Tremaine, S. (1990) Symmetric Multistep Methods for the Numerical Integration of Planetary Orbits. The Astronomical Journal, 100, 1694-1700.

https://doi.org/10.1086/115629