Ordinary Differential Equations

Verlet

The Verlet method takes in: a function a[x,t] that computes the acceleration given the current position and time, the initial position x0 and velocity v0, the time step dt and number of steps to take, Ns.  As output, it returns ordered pairs of {t, x} where x is the numerical approximation to the solution of (((d^2) )/(dt^2))x(t)= a(x,t) at time t.

Verlet[a_, x0_, v0_, dt_, Ns_] := Module[{xout,xnext, xcur,xprev,t,index},
xout = Table[0.0,{j,1,Ns}];
xcur = x0;
xprev = x0 - dt v0 +(1/2) a[x0,0.0] dt^2;
t = 0.0;
xout[[1]] = {t,x0};
For[index = 2, index ≤ Ns, index=index+1,
xnext = 2.0 xcur - xprev + dt^2 a[xcur, t];
t = t + dt;
xout[[index]] = {t, xnext};
xprev = xcur;
xcur = xnext;
];
Return[xout];
]

Example

We’ll take a linear acceleration with frequency (2 π):

a[x_,t_] :=-(2 Pi)^2 x

Use 100 steps to go from 0 ⟶ 1 in time, starting from rest at x0 =1:

Nsteps = 100;
X = Verlet[a, 1.0,0.0,1.0/(Nsteps-1),Nsteps];

A plot of the numerical solution:

G0 = ListPlot[X]

and we can plot the continuous exact solution together with our discrete, numerical solution:

G1 = Plot[Cos[2 Pi t],{t,0,1}];
Show[G0,G1]

The residual (difference between exact and numerical solution) is also available

Y =Table[{t,Cos[2 Pi t]},{t,0,1,1.0/(Nsteps-1)}];
res = Table[{X[[j,1]], X[[j,2]] - Y[[j,2]]},{j,1,Length[X]}];
ListPlot[res]

Midpoint RK

The RKODE function generates numerical solutions to ((df(x))/dx) = G(x,f(x)) given: the RHS function G, the starting value of x (called x0 as input), and the initial condition f(x0) (called f0), the step size dx and number of steps to take, Ns.  RKODE works for scalars or vectors, and returns a list with the pairs {xn, fn} where fn is the approximation to the solution evaluated at xn: fn ~ f(xn).

RKODE[G_, x0_,f0_, dx_, Ns_] := Module[{k1,k2,retvals,index,xn,fn},
retvals = Table[{0.0,0.0},{i,1,Ns}];
xn = x0;
fn = f0;
retvals[[1]] = {x0,f0};
For[index = 2, index ≤ Ns, index=index+1,
k1 = dx G[xn,fn];
k2 = dx G[xn + dx/2, fn + k1/2];
fn = fn + k2;
xn = xn + dx;
retvals[[index]] = {xn,fn};
];
Return[retvals];
]

Example

For an acceleration a = -(2 π)^2 x, as above, we define f = {x[t], x’[t]} (so the first component of f is the position, the second, the velocity) and then G is:

G[x_,f_] := {f[[2]],-(2 Pi)^2 f[[1]]}

Again run for 100 steps from x =0 ⟶ 1 (x is the generic name for the independent variable here) starting from f0 = {1,0} (i.e. initial position at 1, initial velocity zero):

Nsteps = 100;
fvals = RKODE[G,0.0,{1,0}, 1.0/(Nsteps-1),Nsteps];
X = Table[{fvals[[j,1]],fvals[[j,2,1]]},{j,1,Length[fvals]}];

Plot the numerical and exact solutions together,

G0 = ListPlot[X];
G1 = Plot[Cos[2 Pi t],{t,0,1}];
Show[G0,G1]

and the residual:

Y =Table[{t,Cos[2 Pi t]},{t,0,1,1.0/(Nsteps-1)}];
res = Table[{X[[j,1]], X[[j,2]] - Y[[j,2]]},{j,1,Length[X]}];
ListPlot[res]

Unstable Euler

Define the matrix from (2.86) for the explicit Euler update:

UD = {{1,dt},{-w^2 dt,1}};

Set the parameters ω and dt

w = 2 Pi;
dt = .01;

Use the command “MatrixPower” to raise the matrix UD to a power and generate the solutions in (2.87):

f = Table[MatrixPower[UD,j].{1,0},{j,0,1000}];
X = Table[{j dt, f[[j,1]]},{j,1,Length[f]}];

We can see the growth behavior in the position variable X:

ListPlot[X]

To make the point, we’ll connect the plot, avoiding the points with their irrelevant “bunching” - this is Figure 2.3:

ListLinePlot[X,PlotStyle->GrayLevel[0.5]]

Stable Euler

Define the (inverse of the) matrix from (2.88) for the implicit Euler update:

UD = Inverse[{{1,-dt},{w^2 dt,1}}];

Set the parameters ω  and dt

w = 2 Pi;
dt = .01;

Again using “MatrixPower” to generate the sequence:

f = Table[MatrixPower[UD,j].{1,0},{j,0,1000}];
X = Table[{j dt, f[[j,1]]},{j,1,Length[f]}];

This time, the amplitude decays over time:

ListPlot[X]

Or, connecting the points (Figure 2.4):

ListLinePlot[X,PlotStyle->GrayLevel[0.5]]