FEM does not match exact solution

Question

I am trying to solve :

$$-u''(x) + u(x) = \sin(2\pi x)\, ,\quad 0<x<1\, ,$$ $t>0$, with $u(0) = u(1) = 0$. That has as exact solution

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

But the Forward Euler approximation solution does not match the exact solution.

Any help?

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

Comments above are right: it seems that you are also integrating in time (and indeed you also set the numer of points in time in your code), but the equation is only in variable $x$.

The following snippet produce the correct solution to your problem with linear elements in Python.

To compute $\int_0^1 \phi_i(x) f(x)dx$ I used integrate.quad from scipy, which performs Gaussian quadrature. That integral has been computed splitting the interval in $[x_{i-1},x_{i}]$ and $[x_{i},x_{i+1}]$, since the hat functions are non-differentiable at $x=x_i$. Btw, notice that this integral may be computed analytically, but quadrature is usually the choice in a fem solver.

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