矩阵向量乘法，其中向量已经进行了插值 - Python

Question

矩阵向量乘法，其中向量已经进行了插值 - Python

pythonscipyinterpolationfinite-element-analysis

3

我使用有限元方法来逼近Laplace方程 $-u_{xx} = sin(\pi*x)$ ，将其转换为矩阵系统AU=F，其中A是刚度向量，解出U（对我的问题不是非常重要）。现在我已经得到了近似的U，在计算AU时应该得到类似F的向量，其中F为：

AU在x=0到x=1范围内（例如20个节点）给出以下绘图：

然后我需要将U插值为更长的向量，并找到AU（对于更大的A也是如此，但不插值）。我通过以下方式进行插值：

U_inter = interp1d(x,U)
U_rich = U_inter(longer_x)

直到我将其与更长的A矩阵相乘，看起来都运行正常：

看来每个峰值都出现在x的一个节点上（即原始U的节点）。有人知道这是什么原因吗？以下是我的查找A、U和F的代码。

import numpy as np
import math
import scipy
from scipy.sparse import diags
import scipy.sparse.linalg
from scipy.interpolate import interp1d
import matplotlib
import matplotlib.pyplot as plt

def Poisson_Stiffness(x0):
    """Finds the Poisson equation stiffness matrix with any non uniform mesh x0"""

    x0 = np.array(x0)
    N = len(x0) - 1 # The amount of elements; x0, x1, ..., xN

    h = x0[1:] - x0[:-1]

    a = np.zeros(N+1)
    a[0] = 1 #BOUNDARY CONDITIONS
    a[1:-1] = 1/h[1:] + 1/h[:-1]
    a[-1] = 1/h[-1]
    a[N] = 1 #BOUNDARY CONDITIONS

    b = -1/h
    b[0] = 0 #BOUNDARY CONDITIONS

    c = -1/h
    c[N-1] = 0 #BOUNDARY CONDITIONS: DIRICHLET

    data = [a.tolist(), b.tolist(), c.tolist()]
    Positions = [0, 1, -1]
    Stiffness_Matrix = diags(data, Positions, (N+1,N+1))

    return Stiffness_Matrix

def NodalQuadrature(x0):
    """Finds the Nodal Quadrature Approximation of sin(pi x)"""

    x0 = np.array(x0)
    h = x0[1:] - x0[:-1]
    N = len(x0) - 1

    approx = np.zeros(len(x0))
    approx[0] = 0 #BOUNDARY CONDITIONS

    for i in range(1,N):
        approx[i] = math.sin(math.pi*x0[i])
        approx[i] = (approx[i]*h[i-1] + approx[i]*h[i])/2

    approx[N] = 0 #BOUNDARY CONDITIONS

    return approx

def Solver(x0):

    Stiff_Matrix = Poisson_Stiffness(x0)

    NodalApproximation = NodalQuadrature(x0)
    NodalApproximation[0] = 0

    U = scipy.sparse.linalg.spsolve(Stiff_Matrix, NodalApproximation)

    return U

x = np.linspace(0,1,10)
rich_x = np.linspace(0,1,50)
U = Solver(x)
A_rich = Poisson_Stiffness(rich_x)
U_inter = interp1d(x,U)
U_rich = U_inter(rich_x)
AUrich = A_rich.dot(U_rich)
plt.plot(rich_x,AUrich)
plt.show()

- James

你能修改你的代码，使其可以独立运行吗？另外，我在你的代码中没有看到任何插值部分（interp1d），那在哪里呢？ - user707650

嗨，我现在已经包含了我的代码和用于找到图形的程序（包括我使用的interp1d），谢谢。 - James

1

通过微小的缩进更正，此代码便是完整的 - 它可以通过简单地复制和粘贴来运行。 - hpaulj

1

如果您使用不同的插值方法，它似乎可以工作。尝试使用“kind='cubic'”。 - Moritz

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- hpaulj · Accepted Answer

注释1:

我添加了Stiffness_Matrix = Stiffness_Matrix.tocsr()语句，以避免效率警告。有限元计算已经足够复杂，我需要打印一些中间值才能确定发生了什么。

注释2:

plt.plot(rich_x,A_rich.dot(Solver(rich_x)))绘制得很好。你得到的噪音是由于插值U_rich和真实解U_rich-Solver(rich_x)之间的差异造成的。

注释3:

我认为你的代码没有问题。问题在于你可以用这种方式测试插值。我对有限元理论有点陌生，但我认为你需要使用形状函数来进行插值，而不是简单的线性函数。

注释4:

直觉上，A_rich.dot(U_rich)询问的是，产生U_rich的强制项F是什么样子。与Solver(rich_x)相比，U_rich具有平坦区域，其值小于真实解。那么，什么样的F会产生这种情况呢？一个尖锐的、在点处具有NodalQuadrature(x)的函数，但在中间接近零。这就是你的图表所显示的内容。

更高阶的插值将消除平坦区域，并产生一个更平滑的反向计算F。但你真的需要重新审视有限元理论。

你可能会发现查看以下内容很有启示性：

plt.plot(x,NodalQuadrature(x))
plt.plot(rich_x, NodalQuadrature(rich_x))

第二个图形更加平滑，但高度只有第一个图形的五分之一。

更好的比较请看：

plt.plot(rich_x,AUrich,'-*')  # the spikes
plt.plot(x,NodalQuadrature(x),'o')  # original forcing
plt.plot(rich_x, NodalQuadrature(rich_x),'+') # new forcing

在该模型中，强制作用并不是连续的，而是每个节点上的一个值。随着节点数量的增加（rich_x），每个节点上的大小会变小。