Derivation of Finite Element Equations Using Variational Approach

Derivation of Finite Element Equations Using Variational Approach

Yu Yu Paing

Department of Mathematics Polytechnic University (Maubin) Maubin, Myanmar

Abstract – This paper presents the derivation of the discrete finite element system for the Poisson equation. Based on the variational formulation, the domain is partitioned into sub-intervals to define piecewise linear basis functions. The coefficients of the stiffness matrix and the load vector are computed through integration over these intervals, leading to the construction of a linear system of equations. This process demonstrates how a continuous boundary value problem is transformed into a solvable matrix form.

Keywords – Poisson equation; variational formulation; piecewisse linear basis functions

1. INTRODUCTION

   This document describes the systematic derivation of finite element equations from the weak formulation of boundary value problems using variational approach. By defining piecewise linear basis functions, we construct the coefficient stiffness matrix and the load vector required to solve linear system . Furthermore, we demonstrate the numerical implementation of these equations using Matlab to obtain the approximation solution . The Matlab code, numerical results, and corresponding plots are provided to compare the finite element approximation with the exact solution, thereby validating the computational efficiency and accuracy of the method.

2. STIFFNESS MATRIX FOR ELEMENT BY ELEMENT

   We consider the equation

   (1)

   which is a linear system of equations with equations in unknowns, . In matrix form the linear system (1)

   can be written as

   problem, a uniform Cartesian mesh

   where , defining the intervals . By assembling element by element:

   The idea is to break up the integration element by element, so that for any integrable function , we have

   The stiffness matrix can be written

   Note that each interval has only two nonzero basis functions. This leads to

   where is the matrix with elements with

   are -vectors.

   The matrix is called the stiffness matrix and the load vector, with terminology from early applications of FEM in structural mechanics. In the finite element method for the numerical solution of elliptic partial differential equations, the stiffness matrix is a matrix that represents the system of linear equations that must be solved in order to ascertain an approximate solution to the differential equation. For model

   The nonzero contribution from a particular element is

   Thus we finally assemble the tridiagonal matrix

   the two by two local stiffness matrix. Similarly, the local vector is

   and the global load vector can also be assembled element by element:

3. Derivation of finite element equations

   Consider the boundary value problem for the Poisson equation

   in (0, 1),

   The weak formulation is

   for every .

   In the element , there are only two nonzero basis functions centered at and respectively:

   where and are defined only on one particular element. We can concentrate on the corresponding contribution to the stiffness matrix and load vector from the two nonzero basis functions.

   It is easy to verify that

   Therefore, the local stiffness matrix is

   Set and = (0, 1), to construct the finite element approximation of this problem, we subdivide (0,1) into N subintervals with length by the points

   where . The finite element approximation is

   for every

   where ={ is continuous piecewise linear,

   We approximate by , made up of piecewise linear functions. Thus if then must be linear for all and since , must be continuous,

   i.e., for all .

   Further as and . A typical element of will have a graph as shown in Figure 1 below. Since the function is linear in each interval , it is completely determined by its values at and in this interval. Thus the degrees of freedom are located at the nodes

   (since and have fixed values). We

   define the basis as follows:

   and the stiffness matrix is assembled as follows:

   Here, is the continuous piecewise linear function that takes the value 1 at node point and the value 0 at other node points (see Figure 2). Since , it can be written as a linear combination of the basis functions:

   where the coefficients are the unknowns to be determined. Then we derive a linear system of equations for the

   coefficients by substituting the approximate solution for the exact solution in the weak form

   i.e.,

   Next we choose the test function as successively, to get the system of linear equations

   or in the matrix-vector form:

4. EXAMPLE FOR FINITE ELEMENT SOLUTIONS USING MATLAB CODES

Consider the problem

Let I = (1, 2) be divided into a uniform mesh , below we

calculate the finite element approximation for . Firstly, we find the variational formulation to the problem by multiplying (2) by a test function such that

and integrate over the interval (1, 2) we get

for every .

The problem is now to find such that

for every .

Using the linear basis functions for , can be written as a linear combination

Thus, since .

where

Since and , we

obtain

The matrix is symmetric (i.e., ) and positive definite (i.e., ). We can observe that if . Thus the matrix is tridiagonal,

i.e., only the elements in the main diagonal and the two adjoining diagonals may be different from zero.

We also have

and

Thus we get the linear system:

The above matrix symmetric, positive definite and tridiagonal. It can be shown to be monotone, i.e., its inverse has all entries positive. Thus if then as for all , the right hand side vector will have all components positive and so

> 0 for all .

The following Figure 3 is the approximation with depicted with the exact solution . The Matlab code, numerical solution, and the corresponding result plot (see Figure 3) are presented below.

Matlab Code for Example

% Define problem parameters a=1; b=2; %Interval (1, 2) N=5; h=1/N; %Mesh

%Generate Nodes x=linspace(a, b, N+1); disp( Nodes(xvalues): )

%Display values disp(x)

disp([ Mesh h= ,num2str(h)]); A=zeros(N-1, N-1);

for i=1:N-1 A(i, i)=2/h; if i>1

A(i, i-1)=A(i, i-1)-1/h;

A(i-1, i)=A(i-1, i)-1/h;

end end

disp( Stiffness Matrix A: ); disp(A);

%Initialize load vector F=zeros(N-1, 1);

%Compute load vector for i=1:N

xL=x(i); xR=x(i+1); xm=(xL+xR)/2; f_mid=6*xm;

if i>1

F(i-1)=F(i-1)+f_mid*h/2; end

if i<=N-1 F(i)=F(i)+f_mid*h/2; end

end

%Display load vector disp( Load Vector F: ); disp(F);

%Solve for FEM approximation U=A\F;

%Display FEM solution

disp( FEM Approximation U: ) disp(U)

%plot x1=x y1=[0;U;0]

%Exact solution a=1:0.01:2;

u_exact=-a.^3+7*a-6; z=u_exact;

plot(x1,y1, — ,a,z, – )

title( Solid line: Exact solution, Dotted line: FEM

0.7680

Figure 1

Figure 2

Figure 3

ACKNOWLEDGMENT

Approximation solution ) xlabel( x )

ylabel( u_{h}(x) & u(x) ) Result

>> FEM1

Nodes(xvalues):

1.0000 1.2000 1.4000 1.6000 1.8000 2.000

Mesh h=0.2 Stiffness Matrix A:

Load Vector F: 1.4400

1.6800

1.9200

2.1600

FEM Approximation U: 0.6720

1.0560

1.1040

I would like to thank Professor Dr Aung Kyaw, Head of the Department of Mathematics, University of Yangon for giving the pportunity to do this research and Dr Aung Zaw Myint, Department of Mathematics, University of Mandalay for his valuable advice and guidance throughout this study.

REFERENCES

1. J. Kiusalaas, Numerical methods in engineering with MATLAB, Pennsylvania State University, 2005.
2. R. Pretap, Getting started with MATLAB: A quick introduction for science and engineers, Oxford Univeristy Press, 2013.
3. Z. Li, Z. Qiao and T. Tang, Numerical solution of differential equations, Cambridge University Press, Cambridge, 2018.