finite-element-analysis - 初级有限元代码无法正确求解方程

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

离散我的空间:

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) 您正在使用弱公式，尽管您已经隐式地这样做了。一个“弱”的公式与所涉及的导数的顺序无关。它很弱，因为您没有在每个位置都完全满足微分方程。有限元最小化解的加权残差，在域上积分。函数 phi_j 实际上离散了加权函数。当您只有一阶导数时，不同之处在于您不必应用高斯散度定理(该定理简化为一维的分部积分)来消除二阶导数。您可以看出这没有完成，因为 phi_j 在 LHS 中没有区分。

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

其次，关于您的实现:

3) 您使用的集成点数量超出了必要的数量。您的元素使用线性插值函数，这意味着您只需要位于元素中心的一个积分点即可准确计算积分。查看高斯求积的详细信息以了解原因。此外，您还指定了集成点的数量作为节点数量的倍数。这应该以元素数量的倍数来完成(在您的例子中为 n_vertex-1)，因为元素是您要集成的域。

4) 您只需从公式中删除两个端节点即可构建系统。这不是指定边界条件的正确方法。我建议首先构建完整的系统，并使用应用狄利克雷边界条件的典型方法之一。另外，请考虑约束两个节点对于您要求解的微分方程意味着什么。存在什么函数满足 x' = 1、x(0) = 0、x(5) = 0？您尝试将 2 个边界条件应用于一阶微分方程，从而对系统产生了过度约束。

不幸的是，无法进行任何小调整来使代码正常工作，但我希望上面的评论可以帮助您重新考虑您的方法。

编辑以解决您的更改:

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 的方程。积分后执行此操作。