初学者有限元代码不能正确求解方程

Question

我正在尝试编写求解极难微分方程的代码： x' = 1 用有限元法。据我了解，我可以获得解决方案 u

与基函数phi_i(x)，而我可以获得u_i作为线性方程组的解：

用微分算子D（这里只有一阶导数）。作为基础，我正在使用帐篷功能：

    def tent(l, r, x):
    m = (l + r) / 2 
    if x >= l and x <= m:
        return (x - l) / (m - l)
    elif x < r and x > m:
        return (r - x) / (r - m)
    else:
        return 0

def tent_half_down(l,r,x):
    if x >= l and x <= r:
        return (r - x) / (r - l)
    else:
        return 0

def tent_half_up(l,r,x):
    if x >= l and x <= r:
        return (x - l) / (r - l)
    else:
        return 0

def tent_prime(l, r, x):
    m = (l + r) / 2 
    if x >= l and x <= m:
        return 1 / (m - l)
    elif x < r and x > m:
        return 1 / (m - r)
    else:
        return 0

def tent_half_prime_down(l,r,x):
    if x >= l and x <= r:
        return - 1 / (r - l)
    else:
        return 0

def tent_half_prime_up(l, r, x):
    if x >= l and x <= r:
        return 1 / (r - l)
    else:
        return 0

def sources(x):
    return 1

离散化我的 space：

n_vertex = 30
n_points = (n_vertex-1) * 40
space = (0,5)
x_space = np.linspace(space[0],space[1],n_points)
vertx_list = np.linspace(space[0],space[1], n_vertex)
tent_list = np.zeros((n_vertex, n_points))
tent_prime_list = np.zeros((n_vertex, n_points))

tent_list[0,:] = [tent_half_down(vertx_list[0],vertx_list[1],x) for x in x_space]
tent_list[-1,:] = [tent_half_up(vertx_list[-2],vertx_list[-1],x) for x in x_space]
tent_prime_list[0,:] = [tent_half_prime_down(vertx_list[0],vertx_list[1],x) for x in x_space]
tent_prime_list[-1,:] = [tent_half_prime_up(vertx_list[-2],vertx_list[-1],x) for x in x_space]
for i in range(1,n_vertex-1):
    tent_list[i, :] = [tent(vertx_list[i-1],vertx_list[i+1],x) for x in x_space]
    tent_prime_list[i, :] = [tent_prime(vertx_list[i-1],vertx_list[i+1],x) for x in x_space]

计算线性方程组：

b = np.zeros((n_vertex))
A = np.zeros((n_vertex,n_vertex))
for i in range(n_vertex):
    b[i] = np.trapz(tent_list[i,:]*sources(x_space))
    for j in range(n_vertex):
        A[j, i] = np.trapz(tent_prime_list[j] * tent_list[i])

然后求解重构

u = np.linalg.solve(A,b)
sol = tent_list.T.dot(u)

但这不起作用，我只得到一些上下图案。我究竟做错了什么？

Answer 1

首先，关于术语和符号的几点评论：

1) 你正在使用弱公式，尽管你已经隐式地这样做了。 "weak" 的公式与所涉及的导数顺序无关。它很弱，因为您不能在每个位置都精确地满足微分方程。 FE 最小化解决方案的加权残差，在域上集成。函数 phi_j 实际上将加权函数离散化。当您只有 first-order 个导数时，不同之处在于您不必应用高斯散度定理（它简化为一维的部分积分）来消除 second-order 个导数。你可以看出这没有完成，因为 phi_j 在 LHS 中没有区分。

2) 我建议不要使用 "A" 作为微分运算符。您还将此符号用于全局系统矩阵，因此您的符号不一致。人们经常使用"D"，因为这更符合它用于区分的想法。

其次，关于您的实施：

3) 您正在使用比必要更多的集成点。您的元素使用线性插值函数，这意味着您只需要一个位于元素中心的积分点即可准确计算积分。查看高斯正交的详细信息以了解原因。此外，您已将集成点数指定为节点数的倍数。这应该作为元素数量的倍数来完成（在您的情况下，n_vertex-1），因为元素是您要集成的域。

4) 您已经通过简单地从公式中删除两个末端节点来构建您的系统。这不是指定边界条件的正确方法。我建议先构建完整的系统，然后使用一种典型的方法来应用 Dirichlet 边界条件。此外，考虑约束两个节点对您要求解的微分方程意味着什么。存在什么函数满足 x' = 1, x(0) = 0, x(5) = 0？您尝试将 2 个边界条件应用于 first-order 微分方程，从而过度约束了系统。

遗憾的是，没有任何小的调整可以让代码正常工作，但我希望上面的评论能帮助您重新考虑您的方法。

编辑以解决您的更改：

1) 假设矩阵 A 是用 A[row,col] 寻址的，那么你的索引是向后的。您应该与 A[i,j] = ...

整合

2) 应用约束的一种简单方法是用所需的约束替换一行。例如，如果您希望 x(0) = 0，则为所有 j 设置 A[0,j] = 0，然后设置 A[0,0] = 1 并设置 b[0] = 0。这将替换其中一个u_0 = 0 的方程。积分后执行此操作。

初学者有限元代码不能正确求解方程

Beginner Finite Elemente Code does not solve equation properly

finite-element-analysis