Is my Finite Element Model for this problem correct?

Question

The equation I'm trying to solve is

$$ \frac{d^2}{dx^2}u = - \frac{1}{1+x^2} $$

on the domain $0 < x < 1$, with boundary conditions $u(0) = 0$ and $\frac{du}{dx}(1)=0$.

I got the weak form equation as:

$$\int_0^1 \frac{du}{dx} \frac{d\psi}{dx} \ dx =\int_0^1\left( \frac 1 {1+x^2} \right)(\psi) \, dx$$

I assumed a linear element.My element stiffness matrix is:

$$ k^e_{ij} = \frac{1}{h} \left[ \begin{matrix} 1 & -1 \\ -1 & 1 \end{matrix} \right] $$

and

$$ f_i^e = \left( \begin{matrix} \arctan h - \frac 1 {2h}\ln(1+h^2) \\ \frac 1 {2h} \ln(1+h^2) \end{matrix} \right) $$

Here are the results for a 10 element model of a python program that I did:

Exact   FEM
0.07355 0.09468
0.13721 0.17939
0.19127 0.25413
0.23617 0.31891
0.27245 0.37372
0.30073 0.41856
0.32166 0.45343
0.33587 0.47834
0.34399 0.49329
0.34657 0.49826

Are these results to be expected? Or is there a mistake in the model/matrices above?

Answer 1

Your finite element solution is not correct. For this kind of problem (1 -dimensional problem, constant material parameter) the finite element solution should be exact at the vertices. So the floating point numbers in your lists should be equal.

It's hard to say where the error is without seeing your full code. Some possible errors:

• Wrong load vector. I don't know how you evaluated yours, I usually do it using numerical quadrature so I don't use any analytical expressions.

• Error in the boundary condition.

Edit: I added the following test I made to confirm my claim that the finite element solution should be exact at vertices. The following code uses sp.fem python library to compute the solution with linear elements. The solution is compared to the exact solution as given by Mathematica.

import numpy as np
import matplotlib.pyplot as plt 
from spfem.mesh import MeshLine
from spfem.asm import *
from spfem.element import *
from spfem.utils import direct

N = 11
m = MeshLine(np.array([np.linspace(0,1,N)]), np.array([range(N-1),range(1,N)]))

e = ElementLineP1()
a = AssemblerElement(m, e)

A = a.iasm(lambda du,dv: du*dv)
f = a.iasm(lambda v,x: 1.0/(1.0+x**2)*v)

x = direct(A,f,I=np.setdiff1d(np.arange(A.shape[0]),[0]))

exact = lambda x: 1./4.*(np.pi*x-4.0*x*np.arctan(x)+2.0*np.log(1.0+x**2))
y = np.linspace(0,1,30)

m.plot(x)
plt.plot(y,exact(y),'b-')
plt.legend(['finite element','exact'])
m.show()

The following plot is outputted: