PiecewiseFunctionApproximation

Download this page as a Matlab LiveScript

Table of contents

Approximating a function using piecewise elements
Matlab implementation of a piecewise linear approximation
\(L_2\) error measure
Error convergence study

Approximating a function using piecewise elements

So we want to approximate the function

\[u\left(x\right)=2x\;\sin \left(2\pi x\right)+3\;\textrm{for}\;x\in \left\lbrack 0,1\right\rbrack\]

using piecewise elements of order 1 or 2.

This means finding all nodal values \(u_i\) to the approximation (the blue circles in the figure above). This problem can be formulated using Galerkin as such: Find \(u_i\) such that

\[\int_{0}^{1}r(x,u_{i})\varphi_{i}dx=0\] \[=\int_0^1 \left(U\left(x,u_i\right) - u\left(x\right)\right)\varphi_j \;d\;x\;=0\]

with \(U\left(x,u_i\right)=\sum_i \varphi_i \left(x\right){u}_i\)

Now, this is the same old Galerkin expression as we had in the previous section. In order to apply a piecewise approximation, we shall further expand the expressions to formulate a matrix system!

Expanding the expression in the integral we get

\[\int_0^1 \left(\sum_i \varphi_i \left(x\right)\varphi_j \left(x\right){\;u}_i -u\left(x\right)\varphi_j \left(x\right)\right)\;d\;x\;=0\] \[=\int_0^1 \sum_i \varphi_i \left(x\right)\varphi_j \left(x\right){\;u}_i \;d\;x-\int_0^1 u\left(x\right)\varphi_j \left(x\right)\;d\;x=0\]

move \(u_i \;\)out

\[\sum_i \int_0^1 \varphi_i \left(x\right)\varphi_j \left(x\right)\;d\;x\;\;{\;u}_i =\int_0^1 u\left(x\right)\varphi_j \left(x\right)\;d\;x\]

which is equivalent to the matrix system

\(M_{\textrm{ij} } u_i =f_i\) or \(\mathit{\mathbf{M} }\;\mathit{\mathbf{u} }=\mathit{\mathbf{f} }\)

Solving this system gets us all nodal values for \(u_i\).

Piecewise approximation - The mesh

To describe any mesh we need a coordinate matrix (vector in 1D)

\[ \mathbf{x}=\begin{bmatrix}x_{1}\\ x_{2}\\ x_{3}\\ x_{4}\\ x_{5} \end{bmatrix} \]

and a nodal connectivity matrix

\[ \mathbf{T}=\begin{bmatrix}1 & 2\\ 2 & 3\\ 3 & 4\\ 4 & 5 \end{bmatrix} \]

4 linear elements will lead to 5 equations to be solved for \(u_i ,i=1,2,\ldotp \ldotp \ldotp ,5\)

\[ \mathbf{u}=\begin{bmatrix}u_{1}\\ u_{2}\\ u_{3}\\ u_{4}\\ u_{5} \end{bmatrix} \]

An element \(e\) between nodes \(x_i\) and \(x_j\) contributes only to equations \(i\;\)and \(j\). We can use the \(e\)'th row of \(\mathbf T\) to see which nodes are associated with the element \(e\).

Example: Element mass matrix for element 3 between node 3 and 4:

\[\begin{array}{l} \;\;\;\;\begin{array}{ccccc} \;1 & \;2 & \;\;3 & \;\;\;4 & \;\;5 \end{array}\\ \begin{array}{c} 1\\ 2\\ 3\\ 4\\ 5 \end{array}\underset{\mathit{\mathbf{M}}}{\underbrace{\left\lbrack \begin{array}{ccccc} & & & & \\ & & & & \\ & & M_{33} & M_{34} & \\ & & M_{43} & M_{44} & \\ & & & & \end{array}\right\rbrack } } \end{array}\]

\[ \def\arraystretch{2.5} \begin{array}{l} M_{33} =\int_{x_3 }^{x_4 } \varphi_3 \varphi_3 dx\\ M_{34} =\int_{x_3 }^{x_4 } \varphi_3 \varphi_4 dx\\ M_{43} =\int_{x_3 }^{x_4 } \varphi_4 \varphi_3 dx\\ M_{44} =\int_{x_3 }^{x_4 } \varphi_4 \varphi_4 dx\\ \end{array} \]

Local Numbering

We can use local node numbering to simplify things:

Now, for each element \(e\) we define \(x_1\) and \(x_2\), compute \(\varphi\) and then compute \(\mathbf{M}^e\) and \(\mathbf{f}^e\)

Element stiffness matrix

\[\mathbf{M}^e =\int_{x_1 }^{x_2 } \boldsymbol\varphi^T \boldsymbol\varphi dx\]

\[ \def\arraystretch{2.5} \mathbf{M}^{e}=\begin{bmatrix}\int_{x_{1}}^{x_{2}}\varphi_{1}\varphi_{1} \;dx & \int_{x_{1}}^{x_{2}}\varphi_{1}\varphi_{2} \;dx\\ \int_{x_{1}}^{x_{2}}\varphi_{2}\varphi_{1} \;dx & \int_{x_{1}}^{x_{2}}\varphi_{2}\varphi_{2} \;dx \end{bmatrix} \]

and element load vector

\[\mathbf{f}^e =\int_{x_1 }^{x_2 } \varphi^T u(x) \;dx\]

\[ \def\arraystretch{2.0} \mathbf{f}^e =\left\lbrack \begin{array}{c} \int_{x_1 }^{x_2 } u(x) \varphi_1 \;dx\\ \int_{x_1 }^{x_2 } u(x) \varphi_2 \;dx \end{array}\right\rbrack \]

Utilizing a piecewise approximation leads to the element mass matrix and element load vector:

\[M_e =\int_{x_1 }^{x_2 } {\varphi }^T \varphi \;d\;x\]

and

\[f_e =\int_{x_1 }^{x_2 } {\varphi }^T u\left(x\right)\;d\;x\]

which need to be assembled into the global mass matrix and load vectors.

Assembly process

For each element \(e\), we add \(\mathbf{M}^e\) and \(\mathbf{f}^e\) to the correct place in \(\mathbf{M}\) and \(\mathbf{f}\).

Element 1: M([1,2],[1,2])

Element 2: M([2,3],[2,3])

Element 3: M([3,4],[3,4])

Element 4: M([4,5],[4,5])

We need to add each element contribution to the global mass matrix.

\[\begin{array}{l} M_{\textrm{ii}}^{\textrm{new}} = M_{\textrm{ii}}^{\textrm{old}} + M_{11}^e \\ M_{\textrm{ij}}^{\textrm{new}} = M_{\textrm{ij}}^{\textrm{old}} + M_{12}^e \\ \vdots \end{array}\]

\[\begin{array}{l} f_i^{\textrm{new}} =f_i^{\textrm{old}} +f_1^e \\ f_j^{\textrm{new}} =f_j^{\textrm{old}} +f_2^e \end{array}\]

or on matrix form

\[ M(\text{row},\text{col}) \]

where row denotes the array containing the row indices and col the column indices. Thus we could write

\[ M([i,j],[i,j]) \]

to mean a sub-matrix of \(M\) with rows \(i\) and \(j\) and columns \(i\) and \(j\).

The corresponding Matlab (or Julia) syntax is

M([i,j],[i,j]) = M([i,j],[i,j]) + Me

All nodes that have several elements will get contributions from each element. We are "accumulating" mass to these nodes from the elements.

Contribution from element 1

\[{\mathit{\mathbf{M}}\left(\left\lbrack 1,2\right\rbrack ,\left\lbrack 1,2\right\rbrack \right)}_{2\times 2} ={\mathit{\mathbf{M}}\left(\left\lbrack 1,2\right\rbrack ,\left\lbrack 1,2\right\rbrack \right)}_{2\times 2} +M_{2\times 2}^{e_1 }\]

\[\begin{array}{l} \;\;\;\;\begin{array}{ccccc} \;\;\;\;\;1 & \;\;\;2 & 3 & 4 & 5 \end{array}\\ \begin{array}{c} 1\\ 2\\ 3\\ 4\\ 5 \end{array}\left\lbrack \begin{array}{ccccc} M_{11}^1 & M_{12}^1 & & & \\ M_{21}^1 & M_{22}^1 & & & \\ & & & & \\ & & & & \\ & & & & \end{array}\right\rbrack \end{array}\]

Contribution from element 2

\[{\mathit{\mathbf{M}}\left(\left\lbrack 2,3\right\rbrack ,\left\lbrack 2,3\right\rbrack \right)}_{2\times 2} ={\mathit{\mathbf{M}}\left(\left\lbrack 2,3\right\rbrack ,\left\lbrack 2,3\right\rbrack \right)}_{2\times 2} +M_{\;}^{e_2 }\]

\[\begin{array}{l} \;\;\;\;\begin{array}{ccccc} \;\;\;\;\;1 & \;\;\;\;\;\;\;\;2 & \;\;\;\;\;\;\;\;3 & 4 & 5 \end{array}\\ \begin{array}{c} 1\\ 2\\ 3\\ 4\\ 5 \end{array}\left\lbrack \begin{array}{ccccc} M_{11}^1 & M_{12}^1 & & & \\ M_{21}^1 & M_{22}^1 +M_{11}^2 & M_{12}^2 & & \\ & M_{21}^2 & M_{22}^2 & & \\ & & & & \\ & & & & \end{array}\right\rbrack \end{array}\]

etc

For four linear elements, this leads to the equation system

\[ \left(\begin{array}{c} 0.083\,u_1 +0.042\,u_2 =0.39\\ 0.042\,u_1 +0.17\,u_2 +0.042\,u_3 =0.85\\ 0.042\,u_2 +0.17\,u_3 +0.042\,u_4 =0.72\\ 0.042\,u_3 +0.17\,u_4 +0.042\,u_5 =0.45\\ 0.042\,u_4 +0.083\,u_5 =0.27 \end{array}\right) \]

Solving the linear system

We end up with the linear system of equations

\[ \mathbf{Mu=f} \]

which is solved by using a linear equations solver, in Matlab we do this by

u = M\f

We will do this in code now using linear element and symbolic integration.

Matlab implementation of a piecewise linear approximation

clear

syms x
ue(x) = 2*x*sin(2*pi*x)+3;

N = 4;
nnod = N+1;
xnod = linspace(0,1,N+1);
nodes = [(1:nnod-1)', (2:nnod)'];

Phi = @(x,x1,x2) [(x2-x)/(x2-x1), (x-x1)/(x2-x1)];

M = zeros(nnod,nnod); F = zeros(nnod,1);
for iel = 1:N
    inod = nodes(iel,:);
    xc = xnod(inod);
    x1 = xc(1); x2 = xc(2);
    phi = @(x)Phi(x,x1, x2);

    Me = int(phi(x).'*phi(x), x,x1,x2);
    fe = int(phi(x).'*ue(x), x,x1,x2);
    
    M(inod,inod) = M(inod,inod) + Me;
    F(inod) = F(inod) + fe;
end

u = M\F;

figure
fplot(ue,[0,1],'r','LineWidth',2); hold on
plot(xnod,u,'-bo')

\(L_2\) error measure

Now, we will compute the average error between the exact function and our Galerkin approximation, we shall use the \(L_2\) norm as a measure of error. This basicaly means squaring the difference between the functions over the entire domain and taking the square root.

\[L_2 \textrm{error}:=\sqrt{\int_0^1 {\left(U-u\right)}^2 \;d\;x}\]

So we need to compute the difference between the function for each element and adding it all up.

\[L_2 \textrm{error}:=\sqrt{\int_0^1 {\left(U-u\right)}^2 \;d\;x}=\sqrt{\sum_e^N \int_e^{\;} {\left(U_e -u\right)}^2 \;d\;x}\]

where

\[U_e \left(x\right)=\varphi \left(x\right){\mathit{\mathbf{u} } }_e \overset{\textrm{lin}\;\textrm{ele} }{=} \left\lbrack \begin{array}{cc} \varphi_1 & \varphi_1 \end{array}\right\rbrack \left\lbrack \begin{array}{c} u_1 \\ u_2 \end{array}\right\rbrack\]

L2err = 0;
for iel = 1:N
    inod = nodes(iel,:);
    x1 = xnod(inod(1)); x2 = xnod(inod(2));
    phi = @(x)Phi(x,x1, x2);
    U = phi(x)*u(inod);
    L2err = L2err + double(int((U-ue)^2, x,x1,x2));
end
L2err = sqrt(L2err)

L2err = 0.1106

Error convergence study

Now we can start changing \(N\) to see how the error goes down. If we plot \(N\) versus the error we get a convergence plot.

Linear approximation

clear
syms x
ue(x) = 2*x*sin(2*pi*x)+3;

NN = [2,4,8,16,32];
L2Err = [];
Time = [];

figure
fplot(ue,[0,1],'r','LineWidth',2); hold on
hp = plot(NaN,NaN,'-b.');

for N = NN
    tic
    nnod = N+1;
    xnod = linspace(0,1,N+1);
    nodes = [(1:nnod-1)', (2:nnod)'];
    
    Phi = @(x,x1,x2) [(x2-x)/(x2-x1), (x-x1)/(x2-x1)];
    
    M = zeros(nnod,nnod); F = zeros(nnod,1);
    for iel = 1:N
        inod = nodes(iel,:);
        xc = xnod(inod);
        x1 = xc(1); x2 = xc(2);
        phi = @(x)Phi(x,x1, x2);
    
        Me = vpaintegral(phi(x).'*phi(x), x,x1,x2);
        fe = vpaintegral(phi(x).'*ue(x), x,x1,x2);
        
        M(inod,inod) = M(inod,inod) + Me;
        F(inod) = F(inod) + fe;
    end
    
    u = M\F;
    t1 = toc;
    Time(end+1) = t1;

    L2err = 0;
    for iel = 1:N
        inod = nodes(iel,:);
        x1 = xnod(inod(1)); x2 = xnod(inod(2));
        phi = @(x)Phi(x,x1, x2);
        U = phi(x)*u(inod);
        L2err = L2err + double(vpaintegral((U-ue)^2, x,x1,x2));
    end
    L2err = sqrt(L2err);

    L2Err(end+1) = L2err;

    set(hp,'XData',xnod, 'YData',u);
    title(sprintf("Approximation using %d linear elements",N))
    pause(0.1)
    drawnow
end

Quadratic approximation

syms x1 x2 x3 x
A = [1 x1 x1^2
     1 x2 x2^2
     1 x3 x3^2];
C = A\eye(3);
Phi = simplify([1 x x^2]*C);
Phi = matlabFunction(Phi);

Time2 = []; L2Err2 = [];

figure;
for N = NN
    tic
    nnod = N*2+1;
    xnod = linspace(0,1,nnod);
    nodes = [(1:2:nnod-2)', (2:2:nnod-1)', (3:2:nnod)'];
    neq = nnod*1;
    
    M = zeros(nnod,nnod); F = zeros(nnod,1);
    for iel = 1:N
        inod = nodes(iel,:);
        xc = xnod(inod);
        x1 = xc(1); x2 = xc(2); x3 = xc(3);
        phi = @(x)Phi(x,x1, x2, x3);
    
        Me = vpaintegral(phi(x).'*phi(x), x,x1,x3);
        fe = vpaintegral(phi(x).'*ue(x), x,x1,x3);
    
        M(inod,inod) = M(inod,inod) + Me;
        F(inod) = F(inod) + fe;
    end
    
    u = M\F;
    t1 = toc;
    Time2(end+1) = t1;
    
    L2err = 0;
    for iel = 1:N
        inod = nodes(iel,:);
        xc = xnod(inod);
        x1 = xc(1); x2 = xc(2); x3 = xc(3);
        phi = @(x)Phi(x,x1, x2, x3);
        U = phi(x)*u(inod);
        L2err = L2err + double(vpaintegral((U-ue)^2, x,x1,x3));
    end
    L2err = sqrt(L2err);
    L2Err2(end+1) = L2err;
    
    
    clf
    fplot(ue,[0,1],'r','LineWidth',1); 
    hold on
    
    for iel = 1:N
        inod = nodes(iel,:);
        xc = xnod(inod);
        x1 = xc(1); x2 = xc(2); x3 = xc(3);
        phi = @(x)Phi(x,x1, x2, x3);
        U = phi(x)*u(inod);
    
        fplot(U, [x1,x3])
        plot(xc,u(inod),'.b');
    end
    title(sprintf("Approximation using %d quadratic elements",N))
%     if N == 2
%         exportgraphics(gca,"PiecewiseQuadratic.gif","Append",false)
%     end
%     for i = 1:10
%         exportgraphics(gca,"PiecewiseQuadratic.gif","Append",true)
%     end
end

Error convergence graphs

figure; hold on
plot(NN,L2Err,'-ko','DisplayName','Linear')
plot(NN,L2Err2,'-ks','DisplayName','Quadratic')
set(gca,'xscale','log','yscale','log')
xlabel("N"); ylabel("L_2 error")
xticks(NN)
xticklabels({'2','4','8','16','32'})
grid on
legend show

Note how the quadratic approximation not only has a lower initial error but also a steeper slope, meaning it converges faster.

figure; hold on
plot(NN,Time,'-ko','DisplayName','Linear')
plot(NN,Time2,'-ks','DisplayName','Quadratic')
set(gca,'xscale','log','yscale','log')
xlabel("N"); ylabel("Time (s)")
xticks(NN)
xticklabels({'2','4','8','16','32'})
yticks([0.5,1,2,4,8,16])
grid on
legend("show",location="northwest")