FEM ne correspond pas à la solution exacte

Question

J'essaye de résoudre:

$$-u''(x) + u(x) = \sin(2\pi x)\, ,\quad 0<x<1\, ,$$ $t>0$, avec $u(0) = u(1) = 0$. Cela a comme solution exacte

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

Mais la solution d'approximation Forward Euler ne correspond pas à la solution exacte.

De l'aide?

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()
```

Answer 1

Les commentaires ci-dessus sont justes: il semble que vous intégrez également dans le temps (et en effet vous définissez également le nombre de points dans le temps dans votre code), mais l'équation est uniquement en variable $x$.

L'extrait de code suivant produit la solution correcte à votre problème avec les éléments linéaires en Python.

Pour calculer $\int_0^1 \phi_i(x) f(x)dx$J'ai utilisé integrate.quadde scipy, qui effectue une quadrature gaussienne. Cette intégrale a été calculée en divisant l'intervalle en$[x_{i-1},x_{i}]$ et $[x_{i},x_{i+1}]$, puisque les fonctions de chapeau ne sont pas différentiables à $x=x_i$. Btw, notez que cette intégrale peut être calculée analytiquement, mais la quadrature est généralement le choix dans un solveur fem.

    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),'-')