# FEM does not match exact solution

I am trying to solve :

$$-u''(x) + u(x) = \sin(2\pi x)\, ,\quad 0 $$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 - x
t   = np.linspace(0, T, Nt+1)    # mesh points in time
dt  = t - t
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
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()
$$$$

• The code refers to initial conditions but the PDE is not time-dependent. Why?
– knl
Sep 20 '20 at 19:38
• I typed your equations, please check that they are OK. Also, please use MathJax for future posts. Sep 20 '20 at 23:02
• You mention that your FE solution does not match the analytic solution. You could describe how they are different or show a plot. Sep 20 '20 at 23:03
• i want to to take different values of like t =0,1,10,100 and see how it the temperature distribution diffuses on these values.
– user37062
Sep 20 '20 at 23:52
• At the least then, the question is poorly explained. You have a code that references a time variable, your description shows something else, and you don't tell us how the solution doesn't match your expectations. How about reducing the code to its minimal form that only solves the problem as described, and show us a picture with what it yields? Sep 21 '20 at 3:56

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-x

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)

` 