3.38. Ordinary Differential Equations¶

3.38.1. Finite Difference Formulas¶

We define the backward difference operator by:

Repeated application gives:

We can also derive a formula for where is any real number, independent of :

Now we can express the following general integral using the function value from either left () or right () hand side of the interval :

Code:

>>> from sympy import var, simplify, integrate
>>> var("nabla t h")
(nabla, t, h)
>>> s = integrate((1-nabla)**(-t/h), (t, 0, h))
>>> simplify(s)
h*nabla/(-log(1 - nabla) + nabla*log(1 - nabla))
>>> s.series(nabla, 0, 5)
h + h*nabla/2 + 5*h*nabla**2/12 + 3*h*nabla**3/8 + 251*h*nabla**4/720 + O(nabla**5)
>>> s2 = s*(1-nabla)
>>> simplify(s2)
-h*nabla/log(1 - nabla)
>>> s2.series(nabla, 0, 5)
h - h*nabla/2 - h*nabla**2/12 - h*nabla**3/24 - 19*h*nabla**4/720 + O(nabla**5)

Keeping terms only to third-order, we obtain:

Similarly:

Code:

>>> from sympy import var
>>> var("f0 f1 f2 f3")
(f0, f1, f2, f3)
>>> nabla1 = f0 - f1
>>> nabla2 = f0 - 2*f1 + f2
>>> nabla3 = f0 - 3*f1 + 3*f2 - f3
>>> 24*(f0 + nabla1/2 + 5*nabla2/12 + 3*nabla3/8)
-59*f1 - 9*f3 + 37*f2 + 55*f0
>>> 24*(f0 - nabla1/2 - nabla2/12 - nabla3/24)
f3 - 5*f2 + 9*f0 + 19*f1

3.38.2. Integrating ODE¶

Set of linear ODEs can be written in the form:

(3.38.2.1)

For example for the Schrödinger we have

Now we need to choose a grid , where is some uniform grid. For example :

where . We also need the derivative, for the exampe above we get:

Now we substitute this into (3.38.2.1):

We can integrate this system from to on a uniform grid :

where and we use some method to approximate the integral, see the previous section.

(3.38.3.1)

The left hand side can be written as:

So the Poisson equation can also be written as:

(3.38.3.2)

Now we determine the values of , and the behavior of and . The equation determines up to an arbitrary constant, so we set and now the potential is determined uniquely.

The 3D integral of the (number) density is equal to the total (numeric) charge, which is equal to (number of electrons). We can then use the Poisson equation to rewrite the integral in terms of :

Let

Then around we get and (for some constant ). As such, is a point charge (delta function) at the origin. From now on, we will assume no point charge, i.e. .

In the limit , we get the equation:

by integrating (and requiring that vanished in infinity to get rid of the integration constant), we get for :

Integrating (3.38.3.2) directly, we get:

We already know that behaves like in infinity, so vanishes. Requiring itself to vanish in infinity, the left hand side simplifies to and we get:

Last thing to determine is . To do that, we expand the charge density and potential (and it’s derivatives) into a series around the origin:

And substitute into the equation (3.38.3.1):

We now multiply the whole equation by and then set . We get , so . We put it back into the equation to get:

This must hold for all , so we get the following set of equations for :

from which we express for all . We already know the values for and from earlier, so overall we get:

in particular:

So we get the following series expansion for and :

Examples¶

It is useful to have analytic solutions to test the numerical solvers. Here we present a few.

Gaussian Charge¶

The Gaussian charge is simply a Gaussian, normalized in such a way that the total charge is :

(3.38.3.3)

Let us verify the normalization by calculating the total charge :

So the total charge is , as expected. Code:

>>> from sympy import var, integrate, exp, Symbol, oo
>>> var("r")
r
>>> alpha=Symbol("alpha", positive=True)
>>> integrate(exp(-alpha**2*r**2)*r**2, (r, 0, oo))
sqrt(pi)/(4*alpha**3)

Now we calculate the potential from the Poisson equation (3.38.3.2):

We have two integration constants and . We fix the potential using the condition , which implies . The other constant is a point charge at the origin, which in our case (3.38.3.3) is zero, so .

We finally obtain the potential:

We can calculate the electrostatic self-energy, i.e. the energy of interaction of the charge with the potential generated by this charge :

Code:

>>> from sympy import var, integrate, exp, Symbol, oo, erf
>>> var("r")
r
>>> alpha=Symbol("alpha", positive=True)
>>> integrate(exp(-alpha**2*r**2)*erf(alpha*r)*r, (r, 0, oo))
sqrt(2)/(4*alpha**2)

Exponential Charge¶

The exponential charge is simply an exponential, normalized in such a way that the total charge is :

(3.38.3.4)

Let us verify the normalization by calculating the total charge :

So the total charge is , as expected.

Now we calculate the potential from the Poisson equation (3.38.3.2):

Similarly as for the Gaussian charge, we require the potential to vanish at infinity, which implies . Then we calculate the point charge at the origin:

We do not have any point charge at the origin, so , from which it follows . We finally obtain:

Let us calculate the self-energy:

Code:

>>> from sympy import var, integrate, exp, Symbol, oo
>>> var("r Z B")
(r, Z, B)
>>> alpha=Symbol("alpha", positive=True)
>>> integrate(exp(-alpha*r)*r**2, (r, 0, oo))
2/alpha**3
>>> V = integrate(-Z*alpha**3/2 * r * exp(-alpha*r), r, r)/r
>>> V.simplify()
-Z*(alpha*r + 2)*exp(-alpha*r)/(2*r)
>>> ((V+B/r).diff(r)*r**2).simplify()
(-2*B*exp(alpha*r) + Z*alpha*r*(alpha*r + 1) + Z*(alpha*r +
2))*exp(-alpha*r)/2
>>> ((V+B/r).diff(r)*r**2).limit(r, 0)
-B + Z
>>> integrate(exp(-alpha*r)*((1-exp(-alpha*r))/r-alpha*exp(-alpha*r)/2)*r**2, (r, 0, oo))
5/(8*alpha**2)

Piecewise Polynomial Charge¶

We will use a second-derivative continuous piecewise polynomial for , normalized in such a way that the total charge is :

(3.38.3.5)

Let us verify the normalization by calculating the total charge :

So the total charge is , as expected.

Now we calculate the potential from the Poisson equation (3.38.3.2):

Similarly as for the Gaussian charge, we require the potential to vanish at infinity, which implies . Then we calculate the point charge at the origin:

We do not have any point charge at the origin, so , from which it follows . Then is calculated from the condition of a continuous first derivative at :

So . Finally, is calculated from the continuous values of :

which implies . We finally obtain:

Let us calculate the self-energy:

Let us also calculate the following integral:

Which agrees with [Pask2012], equation (10c). The following integral over the sphere of radius :

Again in agreement with [Pask2012], the paragraph after equation (17).

Code:

>>> from sympy import var, pi, integrate, solve
>>> var("r r_c Z A B")
(r, r_c, Z, A, B)
>>> n = -21*Z*(r-r_c)**3*(6*r**2+3*r*r_c+r_c**2)/(5*pi* r_c**8)
>>> 4*pi*integrate(n*r**2, (r, 0, r_c))
Z
>>> V = integrate(-4*pi*r*n, r, r)/r
>>> V.simplify()
Z*r**2*(9*r**5 - 30*r**4*r_c + 28*r**3*r_c**2 - 14*r_c**5)/(5*r_c**8)
>>> ((V+A+B/r).diff(r)*r**2).simplify()
-B + 63*Z*r**8/(5*r_c**8) - 36*Z*r**7/r_c**7 + 28*Z*r**6/r_c**6 - 28*Z*r**3/(5*r_c**3)
>>> (V+A).diff(r).subs(r, r_c)
-Z/r_c**2
>>> (V+A).subs(r, r_c)
A - 7*Z/(5*r_c)
>>> A = solve((V+A).subs(r, r_c)-Z/r_c, A)[0]
>>> A
12*Z/(5*r_c)
>>> V = V + A
>>> V.simplify()
Z*(r**2*(9*r**5 - 30*r**4*r_c + 28*r**3*r_c**2 - 14*r_c**5) + 12*r_c**7)/(5*r_c**8)
>>> 2*pi*integrate(n*V*r**2, (r, 0, r_c))
15962*Z**2/(17875*r_c)
>>> 4*pi*integrate(n*(Z/r-V)*r**2, (r, 0, r_c))
10976*Z**2/(17875*r_c)
>>> 4*pi*integrate((Z/r-V)*r**2, (r, 0, r_c))
14*pi*Z*r_c**2/75

Alternatively, one can also calculate this using a Piecewise function:

>>> from sympy import var, pi, integrate, solve, Piecewise, oo, Symbol
>>> var("r Z A B")
(r, Z, A, B)
>>> r_c = Symbol("r_c", positive=True)
>>> n = Piecewise((-21*Z*(r - r_c)**3*(6*r**2 + 3*r*r_c + r_c**2)/(5*pi*r_c**8), r <= r_c), (0, True))
>>> 4*pi*integrate(n*r**2, (r, 0, oo))
Z
>>> V = integrate(-4*pi*r*n, r, r)/r
>>> V.simplify()
Piecewise((Z*r**2*(9*r**5 - 30*r**4*r_c + 28*r**3*r_c**2 - 14*r_c**5)/(5*r_c**8), r <= r_c), (0, True))
>>> ((V+A+B/r).diff(r)*r**2).simplify()
Piecewise((-B + 63*Z*r**8/(5*r_c**8) - 36*Z*r**7/r_c**7 + 28*Z*r**6/r_c**6 - 28*Z*r**3/(5*r_c**3), r <= r_c), (-B, True))
>>> (V+A).diff(r).subs(r, r_c)
-Z/r_c**2
>>> (V+A).subs(r, r_c)
A - 7*Z/(5*r_c)
>>> A = solve((V+A).subs(r, r_c)-Z/r_c, A)[0]
>>> A
12*Z/(5*r_c)
>>> V = V + Piecewise((A, r <= r_c), (0, True))
>>> V.simplify()
Piecewise((Z*(r**2*(9*r**5 - 30*r**4*r_c + 28*r**3*r_c**2 - 14*r_c**5)/r_c**7 + 12)/(5*r_c), r <= r_c), (0, True))
>>> 2*pi*integrate(n*V*r**2, (r, 0, oo))
15962*Z**2/(17875*r_c)
>>> 4*pi*integrate(n*(Z/r-V)*r**2, (r, 0, oo))
10976*Z**2/(17875*r_c)
 [Pask2012] (1, 2) Pask, J. E., Sukumar, N., Mousavi, S. E. (2012). Linear scaling solution of the all-electron Coulomb problem in solids. International Journal for Multiscale Computational Engineering, 10(1), 83–99. doi:10.1615/IntJMultCompEng.2011002201