How to solve BVP second order ODE using ﬁnite elements with linear shape functions and using weak form formulation?

Solve \(y''(x)-3 y(x) = -x^2\) over \(x=0\ldots 1\) with boundary conditions \(x(0)=0\) and \(x(1)=0\) using piecewise linear trial functions.

Let the trial function be the linear function \(\hat {y}(x)=c_1 x+c_2\). The residual is \(R=\hat {y}''(x)-3\hat {y}(x)+x^2\). Let the test function be \(w(x)\). We need to solve for each element \(i\) the following equation \[ I_i = \int w_i \, R_i\,dx \] Using the weak form, we apply integration by parts to obtain \begin {align*} I_i &= \int _{x_i}^{x_{i+1}} -\frac {dw}{dx}\frac {d\hat {y}}{dx} -3 w(x)\hat {y}(x) + w(x) x^2 \,dx + \left [ w(x) \frac {d\hat {y}}{dx} \right ]_{x_i}^{x_{i+1}} \\ &=0 \tag {96.1}\\ \end {align*}

or \begin {align*} I &= \sum _{i=1}^{i=M}{I_i} \\ &= \sum _{i=1}^{i=M}{ \int _{x_i}^{x_{i+1}} \left ( -\frac {dw}{dx}\frac {d\hat {y}}{dx} -3 w(x)\hat {y}(x) + w(x) x^2 \,dx + \left [ w(x) \frac {d\hat {y}}{dx} \right ]_{x_i}^{x_{i+1}} \right ) }\\ &=0 \end {align*}

Where \(M\) is the number of elements. The above is the \(M\) equations we need to solve for the unknowns \(y_i\) at the internal nodes. In Galerkin method, the test function is \(w_1(x)=H_1(x)\) and \(w_2(x)=H_2(x)\) which comes from writing \(w_i(x) = \frac {d\hat {y}(x)}{dy_i}\)

Rewriting the trial function \(\hat {y}(x)\) in terms of unknown values at nodes results in \[ \hat {y}(x) = H_1(x) y_i + H_2(x) y_{i+1} \] Where \(H_1(x)=\frac {x_{i+1}-x}{h}\) and \(H_2(x)=\frac {x-x_i}{h}\) where \(h\) is the length of each element. Hence (96.1) becomes \begin {align} I_i &= \int _{x_i}^{x_{i+1}} -\begin {Bmatrix}H_1'(x)\\H_2'(x)\end {Bmatrix} \begin {Bmatrix}H_1'(x) H_2'(x)\end {Bmatrix} \begin {Bmatrix}y_i \\y_{i+1}\end {Bmatrix} -3 \begin {Bmatrix}H_1(x)\\H_2(x)\end {Bmatrix} \begin {Bmatrix}H_1(x) H_2(x)\end {Bmatrix} \begin {Bmatrix}y_i \\y_{i+1}\end {Bmatrix} + \begin {Bmatrix}H_1(x)\\H_2(x)\end {Bmatrix} x^2 \,dx + \left [ w(x) \frac {d\hat {y}}{dx} \right ]_{x_i}^{x_{i+1}} \tag {96.2} \end {align}

The terms \(\left [ w(x) \frac {d\hat {y}}{dx} \right ]_{x_i}^{x_{i+1}}\) reduce to \(\left [ w(x) \frac {d\hat {y}}{dx} \right ]_{0}^{1}\) since intermediate values cancel each other leaving the two edge values over the whole domain.

But \(H_1'(x)=\frac {-x}{h}\) and \(H_2'(x)=\frac {x}{h}\). Plugging these into (96.2) and integrating gives the following \begin {align} I_i &= \frac {1}{h} \begin {bmatrix} -1 & 1\\ 1 & -1 \end {bmatrix} \begin {Bmatrix}y_i \\y_{i+1}\end {Bmatrix} -h \begin {bmatrix} 1 & \frac {1}{2} \\ \frac {1}{2} & 1 \end {bmatrix} \begin {Bmatrix}y_i \\y_{i+1}\end {Bmatrix} + \frac {1}{12 h} \begin {bmatrix} 3 x_i^4-4 x_i^3 x_{i+1} + x_{i+1}^4 \\ x_i^4-4 x_i x_{i+1}^3 + 3 x_{i+1}^4 \end {bmatrix} + \left [ w(x) \frac {d\hat {y}}{dx} \right ]_{0}^{1} \tag {96.3} \end {align}

The algebra above was done with Mathematica. In the following code, \(x_1\) is \(x_i\) and \(x_2\) is \(x_{i+1}\)

h1 = (x2 - x)/h; 
h2 = (x - x1)/h; 
int1 = -{{D[h1, x]}, {D[h2, x]}}.{{D[h1, x], D[h2, x]}}; 
int2 = -3 {{h1}, {h2}}.{{h1, h2}}; 
int3 = {{h1}, {h2}} x^2; 
sol1 = Integrate[int1, {x, x1, x2}]; 
sol1 /. {(x1 - x2) -> -h, (x2 - x1) -> h}

{{-(1/h), 1/h}, {1/h, -(1/h)}}

sol2 = Integrate[int2, {x, x1, x2}]; 
sol2 /. {(x1 - x2) -> -h, (x2 - x1) -> h}

{{-h, -(h/2)}, {-(h/2), -h}}

sol3 = Integrate[int3, {x, x1, x2}]; 
Simplify[%]; 
% /. {(x1 - x2) -> -h, (x2 - x1) -> h}

{{(3 x1^4 - 4 x1^3 x2 + x2^4)/(12 h)}, 
 {(x1^4 - 4 x1 x2^3 + 3 x2^4)/(12 h)}}

The above equation \(I_i\) is calculated for each element \(i\), resulting in \(M\) equations where \(M\) is the number of elements. Collecting all these local equations into a global stiﬀness matrix and then adjusting the global stiﬀness matrix for boundary conditions by eliminating corresponding rows and columns, it is then solved for the unknowns \(y_i\) as a system \(Ax=f\) using direct solver. The following code illustrates the method.

Matlab

M   = 9;  %number of elements 
N   = M+1; %number of nodes 
dof = 2; %degree of freedom per element 
kk  = zeros(N); 
f   = zeros(N,1); %force vector 
h   = 1/M; %length of each element 
k_local=1/h*[-1 1;1 -1]-h*[1 1/2;1/2 1]; 
 
for i = 1:M  %loop over all elements 
 x1 = (i-1)*h; 
 x2 = i*h; 
 f_local = -1/(12*h)*[3*x1^4-4*x1^3*x2+... 
       x2^4;x1^4-4*x1*x2^3+3*x2^4]; 
 f(i:i+dof-1) = f(i:i+dof-1)+f_local; 
 kk(i:i+dof-1,i:i+dof-1) = ... 
    kk(i:i+dof-1,i:i+dof-1)+k_local; 
end 
%fix first and last rows since these are B.C. 
kk(1,:)     = 0; 
kk(1,1)     = 1; 
kk(end,:)   = 0; 
kk(end,end) = 1; 
f(1)        = 0; 
f(end)      = 0; 
y           = kk\f; %solve 
 
plot(0:h:1,y,'or') 
hold on; 
plot(0:h:1,y,'-'); 
title(sprintf... 
 ('FEM solution using to ODE using %d elements',M')); 
xlabel('x'); 
ylabel('y(x)'); 
grid

The analytical solution to the ODE is below. We can see it is very close to the FEM solution above with 11 elements. More elements leads to more accurate solution.

Mathematica

Analtical solution

ClearAll[y, x] 
sol= y[x]/.First@DSolve[{y''[x]-3 y[x]==-x^2, 
        y[0] == 0, y[1] == 0}, y[x], x]; 
Plot[sol, {x, 0, 1}, Frame -> True, 
      FrameLabel -> {{"y(x)", None}, 
          {"x", "analytical solution"}}, 
      GridLines -> Automatic, 
      GridLinesStyle -> {{Dashed, LightGray}, 
                         {Dashed, LightGray}} 
    ]

4.14 How to solve BVP second order ODE using ﬁnite elements with linear shape functions and using weak form formulation?