微分方程的Python有限差分法

Question

我必须求解 Euler Bernoulli 微分梁方程：

u’’’’(x) = f(x) ; 

(x为光束轴点坐标)

和边界条件：

u(0)=0,    u’(0)=0,    u’’(1)=0,    u’’’(1)=a 

我研究了数值有限差分理论，它将一系列推导表示为：

U’k = (1/2*h)(Uk+1 - Uk-1) 

U’’k = (1/h2)(Uk+1 - 2 Uk + Uk-1) 

U’’’k = (1/2h3)(Uk+2 - 2 Uk+1 + 2 Uk-1 + Uk-2)

U’’’’k = (1/h4)(Uk+2 - 4 Uk+1 + 6 Uk - 4 Uk-1 + Uk-2)

（k+1、k+2等都是下标）

我找到了一个表达如下的脚本：

import numpy as np
from scipy.linalg import solveh_banded
import matplotlib.pyplot as plt

def beam4(n,ffun,a):
    x = np.linspace(0,1,n+1)

    h = 1.0/n

    stencil = np.array((1,-4,6))

    B = np.outer(stencil,np.ones(n))

    B[1,-1] = -2;    B[2,0] = 7;    B[2,-1] = 1;    B[2,-2] = 5

    f = ffun(x)

    f *= h** 4;     f[-1] *= 0.5;     f[-1] -= a*h**3

    u = np.zeros(n+1) 

    u[1:] = solveh_banded(B,f[1:],lower=False)

    return x,u

但我不明白为什么要这样构建系数矩阵：

stencil = np.array((1,-4,6))

B = np.outer(stencil,np.ones(n))

B[1,-1] = -2;    B[2,0] = 7;    B[2,-1] = 1;    B[2,-2] = 5    

f = ffun(x)

f *= h**4;    f[-1] *= 0.5;    f[-1] -= a*h**3  "

提前致谢！！！！

Answer 1

希望这会有所帮助。 （由于这是我第一次发布答案）

这个Hermitian正定带状矩阵的系数是由于应用了鬼节点方法。它是在不损失精度的情况下处理 FDM 边界条件的最有效和最流行的方法之一（这里这些系数通常会给出二阶收敛率）。如果您无法看到矩阵，请检查我下面代码中的“K”矩阵：

from numpy import linspace,zeros,array,eye,dot
from numpy.linalg import solve
from pylab import plot,xlabel,ylabel,legend,show
a       = 0.2          ;b     = 0.0
LX,dx   = 1.0,0.05     ;nx    = int(LX/dx)+1
X       = linspace(0.0,LX,nx)
Fs      = X**2.0
""""""""""""""""""""""""""""""""""""""""""""""""""""""
def calcB(l,b):
    if   b==0:    return [0.0,0.0]
    elif b==1:    return 1.0/l/l/l/l*array([-4.0,7.0,-4.0,1.0])
    elif b==nx-2: return 1.0/l/l/l/l*array([1.0,-4.0,5.0,-2.0])
    elif b==nx-1: return 1.0/l/l/l/l*array([2.0,-4.0,2.0])
    else:         return 1.0/l/l/l/l*array([1.0,-4.0,6.0,-4.0,1.0])
U       = zeros(nx)     ;V       = zeros(nx)
M       = eye(nx)       ;K       = zeros((nx,nx))
F       = zeros(nx)     ;F[nx-2] = -b/dx/dx
F[nx-1] = 2.0*b/dx/dx-2.0*a/dx
for i in range(nx):
    if i == 0:      I   = [i,i+1]
    elif i == 1:    I   = [i-1,i,i+1,i+2]
    elif i == nx-2: I   = [i-2,i-1,i,i+1]
    elif i == nx-1: I   = [i-2,i-1,i]
    else:           I   = [i-2,i-1,i,i+1,i+2]
    for k,j in enumerate(I):
        K[i,j]  += calcB(dx,i)[k]
""""""""""""""""""""""""""""""""""""""""""""""""""""""
pn      = [0]        ;eq2     = pn
eq1     = [i for i in range(nx) if i not in pn]

MM1_    = K[eq1,:]; MM11,MM12 = MM1_[:,eq1],MM1_[:,eq2]
RR      = F+Fs

U[eq2]  = [0.0]
U[eq1]  = solve(MM11,(RR[eq1]-dot(MM12,U[eq2]))) 
######################Plotting#########################
Us = lambda x: x**6.0/360.0+x**3.0/6.0*(a-1.0/3.0)+x**2.0/2.0*(1.0/4.0-a)
plot(X,U,'bo',label='FDM') ;plot(X,Us(X),'g-',label='solution')
xlabel('X'); ylabel('U'); legend(loc='best')
show()

干杯