交错网格二阶导数算子矩阵的有限差分逼近

Question

我正在做一个计算物理作业，我正在寻找一些帮助，因为我卡住了！在

问题是：

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

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

u’’+500*exp（−100*x^2）−u^3=0

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

u’’（x）=f（x）=-exp（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]-2u[i]）/dx2

所以，我要做的第一件事是为u建立一个数组，以便在for循环中敲打它，用上面的等式填充它，我不知道这里应该是什么，我可以用一个零数组，但是我需要一些值来得到u的值

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

是否需要为u创建另一个交错网格，然后使用该网格生成RHS阵列：

^{pr2}$

那不太管用，但有些事情。。。在

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

谢谢， 利昂

Answer 1

牛顿-拉斐逊法的问题是，正如你所说，它需要一些初始值。当我在uni做这个任务时，我把最终用户的输入作为最初的猜测。对于您的特定情况，您可能希望硬编码最初的猜测。在