一般的常微分方程边值问题的有限元方法代码

以下是一个使用有限元方法求解常微分方程边值问题的示例代码：

import numpy as np
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import spsolve

# 定义常微分方程边值问题
def f(x):
    return np.sin(x)

def p(x):
    return 1.0

def q(x):
    return 1.0

def r(x):
    return 1.0

def exact(x):
    return np.sin(x)

# 定义有限元网格
a = 0.0
b = np.pi
n = 10
h = (b - a) / n

x = np.linspace(a, b, n+1)

# 定义局部刚度矩阵和局部载荷向量
def local_stiffness_matrix(i):
    A = np.zeros((2,2))
    A[0,0] = p(x[i]) / h + q(x[i]) * h / 3.0
    A[0,1] = -p(x[i]) / h + q(x[i]) * h / 6.0
    A[1,0] = -p(x[i]) / h + q(x[i]) * h / 6.0
    A[1,1] = p(x[i]) / h + q(x[i]) * h / 3.0
    return A

def local_load_vector(i):
    b = np.zeros(2)
    b[0] = f(x[i]) * h / 2.0 + r(x[i]) * h**3 / 12.0
    b[1] = f(x[i+1]) * h / 2.0 + r(x[i+1]) * h**3 / 12.0
    return b

# 组装全局刚度矩阵和载荷向量
K = csr_matrix((n+1, n+1))
F = np.zeros(n+1)

for i in range(n):
    A = local_stiffness_matrix(i)
    B = local_load_vector(i)
    K[i:i+2, i:i+2] += csr_matrix(A)
    F[i:i+2] += B

# 处理边界条件
K = K[1:n, 1:n]
F = F[1:n] - exact(a+h) * K[:,0]

# 求解线性方程组
U = np.zeros(n+1)
U[1:n] = spsolve(K, F)

# 输出数值解和精确解
print("Numerical\tExact")
for i in range(n+1):
    print("{:.6f}\t{:.6f}".format(U[i], exact(x[i])))

这个代码使用了线性有限元方法，将区间 $[0,\pi]$ 等分成 $10$ 个小区间，每个小区间使用一个线性函数逼近原方程。通过组装局部刚度矩阵和局部载荷向量，得到全局刚度矩阵和载荷向量。然后，通过处理边界条件和求解线性方程组，得到数值解。最后，将数值解和精确解进行比较。