计算科学 - 交错网格的二阶导数算子矩阵的有限差分逼近 - 吾爱随笔录

我正在做一项计算物理作业，我正在寻求帮助，因为我遇到了困难！

问题是：

编写一个函数来创建交错网格的二阶导数算子矩阵的有限差分逼近。给定输入 $N$ （矩阵的大小）和 $\delta x$ （网格间距），函数应该以三个数组的形式返回三对角矩阵 $(a,b,c)$ . 包括 Neumann 边界条件 $u0(−1) = u0(1) = 0$ 通过修改特定元素 $b$ .

作业的目的是编写一些函数来求解这个非线性椭圆方程：

u^{″} + 500 e^{- 100 x^{2}} - u^{3} = 0 .

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

在练习课上，我写了一个更简单的脚本来解决

u^{″} (x) = f (x) = - e^{x} (- 2 + 2 x + 5 x^{2} + x^{3}) .

$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-')

因此，Newton-Raphson 方法使用导数的有限差分逼近进行离散化

u^{″} = (u [i + 1] + u [i - 1] - 2 u [i]) / d x^{2} .

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

所以，我要做的第一件事是为设置某种数组，以便在 for 循环中敲击它并使用上面的等式填充它，我不知道 'u' 应该在这里，我可以使用一个零数组，但我需要一些值才能为获取一些东西。 $u$ $u''$

在上面的示例中，我rhs使用我生成的交错网格创建数组，称为pts.

我是否需要为制作另一个交错网格，然后使用它来制作 RHS 数组： $u$

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

那行不通，但是沿着这些思路...

我也觉得我应该使用边界条件来为 $u$

u^{'} (1) = u^{'} (- 1) = 0

$u'(1) = u'(-1) = 0$ 所以

u^{″} (0) = 0

$u''(0) = 0$

任何指针都将不胜感激，我真的不确定如何正确输入这个额外的变量。 $u$