0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Those are Dirichlet boundary conditions you've written there. The script you have is for a regular grid. $\endgroup$ – Steve Dec 8 '15 at 9:57
2
$\begingroup$

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$.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.