有限元法求解微分方程

有限元法是一种数值解微分方程的方法，它将求解区域分割成许多小的单元，每个单元内部的微分方程可以用简单的代数方程来近似表示。通过将所有单元的代数方程组合起来，可以得到整个求解区域的微分方程的数值解。以下是有限元法求解微分方程的一般步骤：

1. 将求解区域分割成许多小的单元，每个单元内部的微分方程可以用简单的代数方程来近似表示。

2. 对于每个单元，选择一个适当的插值函数来近似表示未知函数。

3. 将每个单元的代数方程组合起来，得到整个求解区域的微分方程的数值解。

4. 对于边界条件，可以采用类似插值函数的方法来近似表示。

5. 通过求解代数方程组得到数值解。

以下是一个使用有限元法求解二阶常微分方程的Python代码示例：

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

def fem_solve(a, b, c, f, xa, xb, n):
    h = (xb - xa) / n
    x = np.linspace(xa, xb, n+1)
    A = lil_matrix((n+1, n+1))
    b = np.zeros(n+1)
    for i in range(1, n):
        A[i, i-1] = a(x[i]) / h**2 - b(x[i]) / (2*h)
        A[i, i] = -2 * a(x[i]) / h**2 + c(x[i])
        A[i, i+1] = a(x[i]) / h**2 + b(x[i]) / (2*h)
        b[i] = f(x[i])
    A[0, 0] = 1
    A[n, n] = 1
    b[0] = xa
    b[n] = xb
    u = spsolve(A.tocsr(), b)
    return x, u

# 示例：求解 y'' + y = x, y(0) = 0, y(1) = 1
a = lambda x: 1
b = lambda x: 0
c = lambda x: 1
f = lambda x: x
xa = 0
xb = 1
n = 10
x, u = fem_solve(a, b, c, f, xa, xb, n)
print(u)