Finite-difference approximation of the 2nd derivative operator matrix for a staggered grid

Question

I'm working on a computational physics assignment and I was looking for some help as I've got stuck!

The question is:

Write a function to create the finite-difference approximation of the 2nd derivative operator matrix for a staggered grid. Given the inputs $N$ (the size of the matrix) and $\delta x$ (the grid spacing), the function should return the tridiagonal matrix in the form of three arrays $(a,b,c)$. Include the Neumann boundary conditions $u0(−1) = u0(1) = 0$ by modifiying specific elements of $b$.

The aim of the assignment is to write a few functions in order to solve this non-linear elliptic equation:

$$u'' + 500e^{−100 x^2} − u^3 = 0\,.$$

In a practice class I wrote a more simple script for solving

$$u''(x) = f(x) = - e^{x}(-2 + 2x + 5x^2 + x^3)\,.$$

from pylab import *
from scipy import interpolate

def EllipticSolver(Nx):
    dx = 1.0/Nx
    pts = linspace(-dx/2.0,1.+dx/2.0,Nx+2)
    soln = zeros(Nx+2)
    rhs = zeros(Nx+2)
    rhs = -exp(pts)*(-2.0 + 2.*pts + 5.*pts**2 + pts**3)
    soln = tridiagonal2(soln,rhs,dx)
    return pts[1:-1], soln[1:-1]

def tridiagonal2(dat,d,dx):
    N = len(dat)
    print(N)
    dx2 = 1.0/(dx*dx)
    a = zeros(N)
    b = zeros(N)
    c = zeros(N)
    a[:] = dx2
    c[:] = dx2
    b[:] = -2.0*dx2

    b[1] = -1.0*dx2
    b[N-2] = -3.0*dx2

    c[1] = c[1]/b[1]
    d[1] = d[1]/b[1]
    for j in range(2,N-1):
        c[j] = c[j]/(b[j] - c[j-1]*a[j])
        d[j] = (d[j] - d[j-1]*a[j])/(b[j] - c[j-1]*a[j])

    print(a)
    print(b)
    print(c)
    print(d)

    dat[N-2] = d[N-2]
    for j in range(N-3,0,-1):
        dat[j] = d[j] - c[j]*dat[j+1]

    return dat

x1, soln1 = EllipticSolver(100)
x2, soln2 = EllipticSolver(200)
x3, soln3 = EllipticSolver(400)

figure(1)
clf()
plot(x1,soln1,'r-')
plot(x2,soln2,'g-')

s2 = interpolate.interp1d(x2,soln2,'cubic')
soln2a = s2(x1)

s3 = interpolate.interp1d(x3,soln3,'cubic')
soln3a = s3(x2)

diff1 = soln1 - soln2a
diff2 = soln2 - soln3a

figure(2)
clf()
plot(x1,diff1,'r-')
plot(x2,4.*diff2,'g-')

So, Newton-Raphson method, discretise using finite-difference approximations of the derivatives

$$u'' = (u[i+1] + u[i-1] - 2u[i]) / dx^2 \,.$$

So, the first thing I'm gonna have to do is set up some kind of array for $u$ in order to whack it in a for loop and populate it using the above equation I don't know what 'u' should be here, I could use an array of zeros, but I need some values in order to get something for $u''$.

In the example above, I create the rhs array using the staggered grid I produced called pts.

Do I need to make another staggered grid for $u$, and then use that to make a RHS array:

rhs = zeros(N) 
rhs = (ui+1 + ui-1 - 2ui) / dx2 

That wouldn't quite work but something along those lines...

I also feel like I should be using the boundary conditions to find some initial conditions for $u$ also

$$u'(1) = u'(-1) = 0$$ so then $u''(0) = 0$

Any pointers would be appreciated, I'm really unsure how to input this extra variable $u$ properly.

Answer 1

On an irregular grid a possible finite difference operator is the tridiagonal matrix $A$ which gives $$[AU]_i = \frac{\dfrac{U_{i+1}-U_i}{\delta_{i}} - \dfrac{U_{i}-U_{i-1}}{\delta_{i-1}}}{1/2(\delta_i+\delta_{i-1})}\,,$$ where $\delta_i = x_{i+1} - x_i$ and $U$ is a vector with $U_i \approx u(x_i)$.

$U$ is the numerical solution you're looking for, so you wont be able to plug it in. You're going to want some sort of function to solve for $U$, so $f(U) = 0$.