我必须求解欧拉-伯努利微分梁方程,它是:
u’’’’(x) = f(x) ;
(x是梁轴点的坐标)
边界条件:
^{pr2}$我研究了数值有限差分理论,它将一系列导数表示为:
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 "
提前谢谢!!!!在
希望这对你有帮助。(由于这是我第一次发布答案)
这种厄米正定带状矩阵的系数是由于采用了鬼节点法。它是在不损失精度的情况下处理FDM边界条件的最有效和最流行的方法之一(这里这些系数通常给出二阶收敛速度)。如果您无法查看矩阵,请检查我下面代码中的“K”矩阵:
干杯
相关问题 更多 >
编程相关推荐