solution of poisson problem with finite element method in 1D

Question

I am going to solve $$-u''=-1-10\delta (x-0.2)$$ $$u(0)=u(1)=0$$ by finite element method in 1D. The exact solution that I got is $$u(x) = \left\{ {\begin{array}{*{20}{c}} {\frac{{{x^2}}}{2} - \frac{{17}}{2}x~~~~~~~~~~if~~~x \in [0,0.2]} \\ {\frac{{{x^2}}}{2} + \frac{3}{2}x - 2~~~~~if ~~~~x \in [0.2,1]} \\ \end{array}} \right.$$ and I put $f(x)=-1-10*dirac(x-0.2)$ in my matlab code. The problem is the solution of the finite element problem is exactly the same as of the solution of equation $-u''=-1$.

Answer 1

The simplest way to implement Dirac loading for linear elements is to simply put a node at $x=0.2$ and add $-10$ to your load vector at the respective slot.

Here is an example code using Python, NumPy and spfem, see last 8 lines or so. You can easily do the same thing in Matlab.

from spfem.mesh import MeshLine                                                          
from spfem.assembly import AssemblerElement                                              
from spfem.element import ElementLineP1                                                  
from spfem.utils import direct                                                           
import numpy as np                                                                       
import matplotlib.pyplot as plt                                                          

# build mesh                                                                             
p=np.arange(0,1.2,0.2)                                                                   
t=np.arange(0,len(p)-1)                                                                  
t=np.vstack((t,t+1))                                                                     
m=MeshLine(np.array([p]),t)                                                              
m.refine(4)                                                                              

# assemble stiffness matrix and load vector                                             
a=AssemblerElement(m,ElementLineP1())                                                    
A=a.iasm(lambda du,dv: du*dv)                                                            
f=a.iasm(lambda v:-1*v)                                                                  

# boundary dofs                                                                          
B=np.array([0,len(p)-1])                                                                 
# interior dofs                                                                          
I=np.setdiff1d(np.arange(m.p.shape[1]),B)                                                

# add -10 to load vector element corresponding to x=0.2                                                               
f0=0*f                                                                                   
f0[1]=-10                                                                                

# solve system                                                                           
x=direct(A,f+f0,I=I)                                                                     

m.plot(x)                                                                                
m.show()                                                                                 

The solution looks like this:

The rationale for this is the following. The load vector corresponds to the linear form

$$-\int_0^1 v \,dx - \int_0^1 10 \delta(x-0.2) v \,dx = -\int_0^1 v\,dx - 10 v(0.2).$$

If you put $i$'th node at $x=0.2$ then the $i$'th element of the load vector is

$$f_i = -\int_0^1 \phi_i \,dx - 10.$$

Thus, the contribution from the Dirac loading is just $-10$.