Solutions
July 15, 2025
Problem 1
Given:
−(a(x)u0 )0 + q(x)u = f (x) x(0, 1), u(0) = u(1) = 0 (1)
3
where, a(x) = x + 1, q(x) = 1, and f (x) = x − 5x − 1 (2)
1. The formula for the exact solution is u(x) = x2 − x. This is obtained by
assuming u(x) = ax2 + bx + c and substituting in equation 1.
2. Equation (1) is expanded as follows:
au00 + a0 u0 − qu + f (x) = 0 (3)
We assume the solution to be of the linear form:
φ = Ni φ i + Nj φ j (4)
where φe is the elemental approximation of u(x) and φi , φj are the values
at nodes i, j [Link] shape functions Ni and Nj are given as:
xj − x x − xi
Ni = , Nj = (5)
L L
The residue of the function,
R(φe ) = aφ00 + aφ0 − qφ + f (6)
(φ = φ(x) etc.) Galerkin assumes that the integral of residual of the
function weighted over the shape function is zero. i.e,
Z xj
Ni (aφ00 + aφ0 − qφ + f ) dx = 0 (7)
xi
Substituting the value of φ from (4) and simpifying, we get
Z xj Z xj
x
a(Ni0 )2 + qNi2 φi + aNi0 Nj0 + qNi Nj φj dx = Ni f dx+[aNi φ0 ]xji
xi xi
(8)
1
[Link]
Figure 1: Cable driven robot configuration
2
[Link]
Figure 2: Problem 1
Similarly, for the second shape function, from the equation:
Z xj
Nj (aφ00 + aφ0 − qφ + f ) dx = 0 (9)
xi
we get:
Z xj Z xj
0 0 x
aNi Nj + qNi Nj φi + a(Nj0 )2 + q(Nj )2 φj dx = Nj f dx+[aNj φ0 ]xji
xi xi
(10)
3. In matrix form, equations (8,10) are:
R xj R xj 0 0 R xj x
a(N 0 )2 + q(Ni )2 dx Ni f dx + [aNi φ0 ]xji
aN N + qNi Nj dx φi
Rxxij 0 i 0 R xxij i 0 2j = xi
x x (11)
Nj f dx + [aNj φ0 ]xji
aNi Nj + qNi Nj dx xi a(Nj ) + q(Nj )2 dx φj
R j
xi xi
The above elemental equation is in the form.
[K] [φ] = [F + I] (12)
The above equation is assembled to the global array and solved to get the
nodal values of u(x).
4. The above equation (11) is assembled and solved using MATLAB. The
integrations are performed using a 2-point quadrature. Corresponding
files are: ‘Problem1a.m, localstiffnessmatrix1a.m, forcematrix1a.m’
5. L-2 error and order of convergence is calculated using MATLAB (same
code): Nodes | L-2 error | Order of convergence
----------------------------------
3
10 9.32635e-05 /
20 2.38711e-05 1.96605
40 6.04043e-06 1.98254
80 1.51942e-06 1.99113
6. For piecewise quadratic element, the equations (4) and (5) will be changed
as:
φ = Ni φi + Nj φj + Nk φk (13)
where
2 −4 2
Ni = (x − xj ) (x − xk ) , Nj = 2 (x − xi ) (x − xk ) , Nk = 2 (x − xi ) (x − xj )
L2 L L
(14)
Proceeding the same way as in the case of linear elements, we get the
matrix equation:
R xk R xk 0 0 R xk
a(N 0 )2 + q(Ni )2 dx [aN 0 N 0 + qNi Nk ] dx
aN N + qNi Nj dx
Rxxik 0 i 0 R xxik i 0 2j R xxki i0 k0 φi
2
aNi Nj + qNi Nj dx a(N j ) + q(N j ) dx aN j N k + qN j N k dx φj
Rxxik x
R xik 0 0 x
R xik
[aNi0 Nk0 + qNi Nk ] dx a(Nk0 )2 + q(Nk )2 dx φk
xi xi
aNj Nk + qNj Nk dx xi
R xk x
N f dx + [aNi φ0 ]xki
R xxik i x
= R xi Nj f dx + [aNj φ0 ]xki (15)
xk x
xi
Nk f dx + [aNk φ0 ]xki
7. The solution to problem using quadratic elements is carried out using the
MATLAB codes ‘Problem1b.m, localstiffnessmatrix1b.m, forcematrix1b.m’
L-2 error and order of convergence is calculated in the code:
Nodes | L-2 error | Order of convergence
----------------------------------
10 3.14852e-15 /
20 1.69003e-14 2.42431
40 9.43906e-14 2.48159
80 4.58948e-13 2.28161
8. From the results, it is clear that the solution converges faster for higher
order elements than the linear element.
Problem 2
Given:
S q
W 00 (x) − W (x) = , 0 < x < l, W (0) = W (l) = 0 (16)
D 2D
where q = 200, S = 100, D = (8.8)107 , l = 50
Note: The formula given in the question for comparison seems to be incor-
rect. The exact formula is:
√S √S √S √S
e− D x −1 + e D x −e D l + e D x q
W (x) = √S (17)
2 1 + e Dl S
4
[Link]
Figure 3: Exact solution
1. The exact solution is plotted using the MATLAB file exactsolution.m
2. The approximate solution to equation (16) is found out using the FEM
formulation as described above.
• For piecewise linear function, the equation (11) is used by setting a
S q
=1 and substituting q = D , f = − 2D . Corresponding MATLAB
files are ’Problem2a.m, localstiffnessmatrix2a.m, forcematrix2a.m’.
• For piecewise quadratic function, the equation (15) is used by setting
S q
a =1 and substituting q = D , f = − 2D . Corresponding MATLAB
files are ’Problem2b.m, localstiffnessmatrix2b.m, forcematrix2b.m’.
3. The results are plotted in figures below for linear and quadratic elements
4. L2- error for linear elements:
Nodes | L-2 error | Order of convergence
----------------------------------
10 5.85028e-10 /
20 1.49705e-10 1.96638
40 3.78801e-11 1.98261
L2- error for quadratic elements:
Nodes | L-2 error | Order of convergence
----------------------------------
10 1.44965e-12 /
20 9.17055e-14 3.98255
40 5.81562e-15 3.979
5
[Link]
Figure 4: Comparison between results
(1.96 + 0.84)2 × 1.092 × 2
n= = 51
0.62
(1.96 + 0.84)2 × 1.092 × 2
n= = 51.4
0.62