计算科学 - FEM 不匹配精确解 - 吾爱随笔录

FEM 不匹配精确解

计算科学有限元 Python

2021-12-07 02:18:16

我正在尝试解决：

- u^{″} (x) + u (x) = \sin (2 π x), 0 < x < 1,

$-u''(x) + u(x) = \sin(2\pi x)\, ,\quad 0<x<1\, ,$

t > 0

$t>0$ , 其中 . 这具有精确的解决方案

u (0) = u (1) = 0

$u(0) = u(1) = 0$

u (x) = \frac{\sin (2 π x)}{1 + 4 π^{2}} .

$u(x) = \frac{\sin(2\pi x)}{1 + 4\pi^2}\, .$

但是前向欧拉近似解与精确解不匹配。

有什么帮助吗？

import numpy as np
import matplotlib.pyplot as plt
L   = 1
Nx  = 19
Nt  = 800 
T   = 0.1
x   = np.linspace(0, L, Nx+1)    # mesh points in space
dx  = x[1] - x[0]
t   = np.linspace(0, T, Nt+1)    # mesh points in time
dt  = t[1] - t[0]
a   = 1
F   = a*dt/dx**2
u   = np.zeros(Nx+1)
u_n = np.zeros(Nx+1)

def I(x):
    return(np.sin( 2*x*np.pi ))

# Set initial condition u(x,0) = I(x)  
for i in range(0, Nx+1):
    u_n[i] = I(x[i])  

for n in range(0, Nt):
    # Compute u at inner mesh points
    for i in range(1, Nx):
        u[i] = u_n[i] + F*(u_n[i-1] - 2*u_n[i] + u_n[i+1])

    # Insert boundary conditions
    u[0]  = 0  
    u[Nx] = 0

    # Update u_n before next step
    u_n[:]= u
exact = np.sin(2*np.pi *x ) / (1+4*np.pi**2)
plt.plot(x,u)
plt.plot(x,exact)
plt.show()
```

1个回答

上面的评论是正确的：您似乎也在时间积分（实际上您还在代码中设置了时间点的数量），但方程仅在变量中。 $x$

以下代码段为 Python 中的线性元素问题提供了正确的解决方案。

为了计算，我使用了 from ，它执行高斯求积。处不可微分，因此已计算出将和中的区间分开的积分。顺便说一句，请注意这个积分可以解析计算，但正交通常是 fem 求解器中的选择。 $\int_0^1 \phi_i(x) f(x)dx$ integrate.quadscipy $[x_{i-1},x_{i}]$ $[x_{i},x_{i+1}]$ $x=x_i$

    import numpy as np
    import matplotlib.pyplot as plt
    from scipy import integrate
    
    M = 10 #points in x
    L = 1 #endpoint
    x = np.linspace(0,L,M+1)
    h = x[1]-x[0]
    
    def uex(x):
        return(np.sin( 2*x*np.pi ))/(1+4*np.pi**2)
    
    
    def stiffassembly(M):
        diag = np.zeros(M-1) #x_1,...,x_{M-1} (M-1)
        offd = np.zeros(M-2) #off diagonal terms
        for i in range(1,M):
            diag[i-1] = 1/h +1/h
    
        for k in range(1,M-1):
            offd[k-1] = -1/h
    
        A = np.diag(offd,-1) + np.diag(diag,0) + np.diag(offd,+1)
        return A
    
    
    def massmatrix(N):
        diag = np.zeros(N-1) #x_1,...,x_M-1 (M-1)
        subd = np.zeros(N-2) 
        supr = np.zeros(N-2)
        
        for i in range(1,N):
            diag[i-1] = (h + h)/3
        for k in range(1,N-1):
            supr[k-1] = h/6
            subd[k-1] = h/6
    
        M = np.diag(subd,-1) + np.diag(diag,0) + np.diag(supr,+1)
        return M
    
    
    
    def f(x):
        return np.sin(2*np.pi*x)
    
    
    def load(M):
        load = np.zeros(M-1)
        for k in range(1,M):
            load[k-1] = integrate.quad(lambda w: f(w)*(1/h)*(w-x[k-1]),x[k-1],x[k])[0] + integrate.quad(lambda w: f(w)*(1/h)*(x[k+1]-w),x[k],x[k+1])[0]
        
        return load
    
    
    u = np.linalg.solve(+stiffassembly(M) + massmatrix(M),load(M))
    U = np.r_[0,u,0]
    plt.plot(x,U,'o',x,uex(x),'-')

其它你可能感兴趣的问题

上一篇基于状态方程的变系数隐式方法下一篇迭代稀疏数据集的有效方法