差分Python的有限差分法方程

Question

我必須解決歐拉伯努利差動梁方程是：

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

希望這是有幫助的。 （由於這是我第一次發佈答案）

這個埃爾米特正定帶狀矩陣的係數是由於應用鬼節點方法。它是處理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() 

歡呼