Discretization methods for ODE

Quantum Mechanics
                      K. Konishi, G. Paffuti

Explore in Mathematica the discretization technique for ordinary differential equations.

ClearAll["Global`*"];
Off[General::spell, General::spell1];

Derivatives

Operator δ

The finite difference operator δ is defined by its action on a function f, defined on some discrete set of points, labeled with an index:

δ[f [j] ] = f[ j+ 1/2] - f[ j - 1/2]

In Mathematica this operator can be defined with

δ[f_[x_]]:= f[x+(1/2)] - f[x-(1/2)];

We can enforce linearity with

δ[a_.f_[x_]]:= a δ[f[x]];
δ[a_.f_[x_] + b_]:=  Expand[a δ[f[x]] + δ[b]];
δ[0] = 0;

This operator clearly approximates the derivative of f[x] at point x=j.

Let us take for instance a uniform grid with spacing h. f[j] is the value of f at the j-th point of the grid. Using Taylor expansion

δ[f [j] ] ≃ f'[j]h + 𝒪((h^3))

Let us note that 𝒪((h^2))corrections are absent due to the symmetric definition.

Applying two times the operator we have an approximation for second derivative

δ[δ[f[j]]]

We will call this operator δ2 or (δ^2):

(δ^2) f[j]= f[j+1]- 2 f[j] + f[j-1]

Using δ it is possible to write an arbitrary approximation to second derivatives. One has an expansion of the form

(δ^2)y[j]= (h^2)(1 + (1/12)(δ^2) - (1/240)(δ^4) + ...)y''[j]

Let us check this expansion. It is evidently enough to check it for polynomials (i.e. Taylor expansions around  x[j] ).

fh = Sum[a[n] j^n h^n,{n,1,8}];
dfh = Sum[n(n-1) a[n] j^(n-2) h^(n-2),{n,2,8}];

The l.h.s. of () is (at the origin)

lhs=δ[δ[f[j]]]/.f[j_]-> fh/.j-> 0

The first to terms inthe r.h.s. of () give an approximation up to order 𝒪( (h^6))

Expand[(h^2)(dfh + (1/12) δ[δ[df[j]]] )/.df[j_]-> dfh/.j-> 0]

while next order approximate lhs to order 𝒪( (h^8))

Expand[(h^2)(dfh + (1/12) δ[δ[df[j]]] - (1/240)δ[δ[δ[δ[df[j]]]]] )/.df[j_]-> dfh/.j-> 0]

Differential equations

Let us consider a differential equation of the form

y''= F * y + G;

Using the expansion () one get, at the second order approximation

(δ^2)y[j]= (h^2)(1 + (1/12)(δ^2) )(F[j]y[j] + G[j])

This is a system of linear equations for y[j] and can be solved by several procedures. In practice a function is defined on a grid of points and can be represented by a vector. For this reason it is convenient to define the operator δ2 directly as a matrix, in order to allow the use of all Linear Algebra machinery in a simple way.

δ2[v_]:= Module[{nL,d2},
nL = Length[v];
d2 = ListConvolve[{1,-2,1},v];
d2=Append[d2=Prepend[d2,v[[2]]-2v[[1]]],-2v[[nL]]+v[[nL-1]]];
d2];

The Mathematica function ListConvolve (see Help) make the second difference when used with the list {1,-2,1}:

ListConvolve[{1,-2,1},{v[1],v[2],v[3],v[4],v[5]}]

On a grid with nL points there are only nL-2 second order differences as the ending points miss one of the neighbours. We will often consider differential equations with homogeneous Dirichlet conditions, i.e. y[0] = y[L] = 0.  The grid of points will consist of internal points:

x[j] = j * h ;          j = 1, . . . nL.

Then we can define the second difference operator by assigning a zero value to f[x] computed at external grid points. This is exactly what the two instructions Prepend and Append realize in the function δ2. As a check

δ2[{v[1],v[2],v[3],v[4],v[5]}]

Let us now show explicitely that applying the operator δ2 we get a system of linear equations

np = 5;

yv = Table[y[i],{i,1,np}]; 
Fv= Table[F[i],{i,1,np}];
Gv= Table[G[i],{i,1,np}];

eq1=Table[δ2[yv][[i]] == Expand[(h^2)(yv*Fv + Gv + (1/12)δ2[yv*Fv + Gv])[[i]]],{i,1,np}];
TableForm[eq1]

Matrix Form

As already noticed it is convenient to write operator in matrix form. It is easy to see from above computations that δ2 matrix has the form

delta2 = (({{-2, 1, 0, 0, 0}, {1, -2, 1, 0, 0}, {0, 1, -2, 1, 0}, {0, 0, 1, -2, 1}, {0, 0, 0, 1, -2}}))

We can use the function SparseArray to get a tridiagonal matrix. The instruction Normal transforms the sparsearray structure in a normal matrix for readibility

delta2= Normal[SparseArray[{{i_,i_} ->  -2,{i_,j_}/;j-i==1-> 1,
{i_,j_}/; i-j== 1->  1},{np,np}]];

delta2//MatrixForm

The matrix form of the rhs of the above equations is generated through

(matF=Normal[SparseArray[{Band[{1,1}]->(10/12) Fv,
Band[{1,2}]-> (1/12)Fv[[2;;]],Band[{2,1}]->(1/12) Fv[[1;;np-1]]}
]])//MatrixForm

In the upper diagonal the index start from 2 while in the lower diagonal it starts from 1. That everything is correct can be checked by comparing the matrix expression

matF.yv//MatrixForm

with the F dependent terms in the above equations

TableForm[eq1]

The terms containing the vector G are written with the matrix 1 + δ2/12, called matS below

(matS=Normal[SparseArray[{Band[{1,1}]->10/12,
Band[{1,2}]-> 1/12,Band[{2,1}]->1/12},{np,np}]])//MatrixForm

In matrix form the above equations take the form

δ2 y = (h^2)matF·y + (h^2)matS·Gv

eq2= Table[(delta2.yv )[[i]] ==Expand[(h^2 matF.yv)[[i]] +( h^2 matS.Gv) [[i]]],{i,1,np}];
TableForm[eq2]

This form is identical to the previous one :

eq1== eq2

Note

The matrix matF is not hermitian.  This is correct as the self-adjoint operator is the second derivative, which is formally multiplied by 1+δ2/12 = matS. To have explicit hermitian matrices we have to multiply the operators by (matS^(-1)). This can be explicitly checked:

MatrixForm/@{Inverse[matS].delta2,Inverse[matS].matF}

This transformation has to be avoided as the difference matrix ceases to be sparse.

Integrals

The operator δ2 can also be used to write approximate integration formulas.The classic Simpson rule

(∫_(Subscript[x, j-1])^(Subscript[x, j+1]))f[x]dx = (h/3)(f[j-1] +4 f[j] + f[j+1])+O((h^5));

can be written

(∫_(Subscript[x, j-1])^(Subscript[x, j+1]))f[x]dx = 2h(f[j] +(1/6)(δ^2)f[j])+O((h^5));

In fact

2h(Simplify[f[j] + (1/6)δ[δ[f[j]]]])

To test the accuracy is enough to consider a polynomial with a zero in (Subscript[x, j]). With our previous example

exactIntegral=Integrate[Sum[a[n] x^n,{n,1,8}],{x,-h,h}]

Expand[Simplify[2h*(f[j] + (1/6)δ[δ[f[j]]])]/.f[j_]-> fh/.j-> 0]

We can have two order more in accuracy by using

(∫_(Subscript[x, j-1])^(Subscript[x, j+1]))f[x]dx = 2h(f[j] +(1/6)(δ^2)f[j]-(1/180)(δ^4)f[j])+O((h^7))=
(h/90)(-f[j+2] + 34 f[j+1] + 114 f[j] + 34 f[j-] - f[j-2])+O((h^7));

Let us first verify the equality in the previous equation:

Simplify[2h(f[j] + (1/6) δ[δ[f[j]]] - (1/180)δ[δ[δ[δ[f[j]]]]] )]

Now we can verify that this approximation coincides with the exact result up to order  𝒪( (h^7))

approxIntegral=Expand[Simplify[2h*(f[j] + (1/6)δ[δ[f[j]]]- (1/180)δ[δ[δ[δ[f[j]]]]])]/.f[j_]-> fh/.j-> 0];

approxIntegral-exactIntegral

Usually the integral is on a finite grid. At the end points the usual Simpson rule can be used

(∫_(Subscript[x, 0])^(Subscript[x, 2]))f[x]dx = (h/3)(f[0]+4 f[1] + f[2])

If the function vanishes at the boundaries the integral can be written in terms of internal points only.