Lagrange basis of order p=2

Contents

Shape functions

phi0 = conv([2 -1],[1 -1]);
phi1 = conv([4 0],[-1 1]);
phi2 = conv([2 -1],[1 0]);

Dphi0 = polyder(phi0);
Dphi1 = polyder(phi1);
Dphi2 = polyder(phi2);

xiax = 0:0.01:1;
figure
plot(xiax,polyval(phi0,xiax),'b',xiax,polyval(phi1,xiax),'g',xiax,polyval(phi2,xiax),'r')

Integrals on the reference element

pr = polyint(conv(phi0,phi0));
M00 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(phi0,phi1));
M01 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(phi0,phi2));
M02 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(phi1,phi1));
M11 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(phi1,phi2));
M12 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(phi2,phi2));
M22 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(Dphi0,Dphi0));
S00 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(Dphi0,Dphi1));
S01 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(Dphi0,Dphi2));
S02 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(Dphi1,Dphi1));
S11 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(Dphi1,Dphi2));
S12 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(Dphi2,Dphi2));
S22 = polyval(pr,1) - polyval(pr,0);

N = M * 3 - (M+1); % dimension of the FEM space (= number of basis functions)

C = zeros(1,N); % fill the vector of cell volumes with zeros
C(1:2:end) = x(2:end) - x(1:end-1); % insert the cell volumes
% Attention: the first element of C is C_0 !! --> C(i) = C_{i-1}
% Only the odd positions are filled.
x_large = zeros(1,N); % vector of vertices with zeros in between.
x_large(1:2:end) = x(1:end-1); % the last vertex is omitted.

nodes = zeros(1,N); % node vector
nodes(1:2:end) = (x(2:end)+x(1:end-1))/2;
nodes(2:2:end-1) = x(2:end-1);

Sparse matrix assembly

nz = 1; % counter of non-zero elements (+1)

i = 1; % boundary conditions on the left:
ivec(nz) = i; % vector of row indices
jvec(nz) = i; % vector of column indices
Svec(nz) = 1/C(i)*S11; % vector of stiffness matrix entries
Mvec(nz) = C(i)*M11; % vector of mass matrix entries
nz = nz+1;
ivec(nz) = i;
jvec(nz) = i+1;
Svec(nz) = 1/C(i)*S12;
Mvec(nz) = C(i)*M12;
nz = nz+1;

i = 2; % boundary conditions on the left:
ivec(nz) = i;
jvec(nz) = i;
Svec(nz) = 1/C(i-1)*S22 + 1/C(i+1)*S00;
Mvec(nz) = C(i-1)*M22 + C(i+1)*M00;
nz = nz+1;
ivec(nz) = i;
jvec(nz) = i+1;
Svec(nz) = 1/C(i+1)*S01;
Mvec(nz) = C(i+1)*M01;
nz = nz+1;
ivec(nz) = i;
jvec(nz) = i-1;
Svec(nz) = 1/C(i-1)*S12;
Mvec(nz) = C(i-1)*M12;
nz = nz+1;
ivec(nz) = i;
jvec(nz) = i+2;
Svec(nz) = 1/C(i+1)*S02;
Mvec(nz) = C(i+1)*M02;
nz = nz+1;

% i odd:
for i=3:2:N-2

    ivec(nz) = i;
    jvec(nz) = i;
    Svec(nz) = 1/C(i)*S11;
    Mvec(nz) = C(i)*M11;
    nz = nz+1;
    ivec(nz) = i;
    jvec(nz) = i+1;
    Svec(nz) = 1/C(i)*S12;
    Mvec(nz) = C(i)*M12;
    nz = nz+1;
    ivec(nz) = i;
    jvec(nz) = i-1;
    Svec(nz) = 1/C(i)*S01;
    Mvec(nz) = C(i)*M01;
    nz = nz+1;

end

% i even:
for i=4:2:N-3

    ivec(nz) = i;
    jvec(nz) = i;
    Svec(nz) = 1/C(i-1)*S22 + 1/C(i+1)*S00;
    Mvec(nz) = C(i-1)*M22 + C(i+1)*M00;
    nz = nz+1;
    ivec(nz) = i;
    jvec(nz) = i+1;
    Svec(nz) = 1/C(i+1)*S01;
    Mvec(nz) = C(i+1)*M01;
    nz = nz+1;
    ivec(nz) = i;
    jvec(nz) = i-1;
    Svec(nz) = 1/C(i-1)*S12;
    Mvec(nz) = C(i-1)*M12;
    nz = nz+1;
    ivec(nz) = i;
    jvec(nz) = i+2;
    Svec(nz) = 1/C(i+1)*S02;
    Mvec(nz) = C(i+1)*M02;
    nz = nz+1;
    ivec(nz) = i;
    jvec(nz) = i-2;
    Svec(nz) = 1/C(i-1)*S02;
    Mvec(nz) = C(i-1)*M02;
    nz = nz+1;

end

i = N-1; % boundary conditions on the right:
ivec(nz) = i;
jvec(nz) = i;
Svec(nz) = 1/C(i-1)*S22 + 1/C(i+1)*S00;
Mvec(nz) = C(i-1)*M22 + C(i+1)*M00;
nz = nz+1;
ivec(nz) = i;
jvec(nz) = i+1;
Svec(nz) = 1/C(i+1)*S01;
Mvec(nz) = C(i+1)*M01;
nz = nz+1;
ivec(nz) = i;
jvec(nz) = i-1;
Svec(nz) = 1/C(i-1)*S12;
Mvec(nz) = C(i-1)*M12;
nz = nz+1;
ivec(nz) = i;
jvec(nz) = i-2;
Svec(nz) = 1/C(i-1)*S02;
Mvec(nz) = C(i-1)*M02;
nz = nz+1;

i = N; % boundary conditions on the right:
ivec(nz) = i;
jvec(nz) = i;
Svec(nz) = 1/C(i)*S11;
Mvec(nz) = C(i)*M11;
nz = nz+1;
ivec(nz) = i;
jvec(nz) = i-1;
Svec(nz) = 1/C(i)*S01;
Mvec(nz) = C(i)*M01;
nz = nz+1;

A = sparse(ivec,jvec,Svec + Mvec); % system matrix
figure
spy(A)
title('A')
Mass = sparse(ivec,jvec,Mvec); % mass matrix
figure
spy(Mass)
title('M')

Right-hand side vector

if swi_rhs=='exact'

    syms xi

    for i = 1:2:N

        clear z igrand han int
        z = x_large(i) + xi*C(i); % is the variable transformation R_{i-1}
        igrand(xi) = eval(rho)*poly2sym(phi1,xi); % the integrand in xi-space
        han = matlabFunction(igrand(xi));

        int = C(i)*integral(han,0,1); % numerical integration

        rhovec(i) = int;

    end

    for i = 2:2:N-1

        clear z igrand han int
        z = x_large(i-1) + xi*C(i-1); % is the variable transformation R_{i-2}
        igrand(xi) = eval(rho)*poly2sym(phi2,xi); % the integrand in xi-space
        han = matlabFunction(igrand(xi));
        int = C(i-1)*integral(han,0,1); % numerical integration

        rhovec(i) = int;

        clear z igrand han int
        z = x_large(i+1) + xi*C(i+1); % is the variable transformation R_{i}
        igrand(xi) = eval(rho)*poly2sym(phi0,xi); % the integrand in xi-space
        han = matlabFunction(igrand(xi));
        int = C(i+1)*integral(han,0,1); % numerical integration

        rhovec(i) = rhovec(i) + int;

    end

elseif swi_rhs=='basis'

    z = nodes; % set the symbolic variable
    rhovec = (Mass * eval(rho)')';

end

Solve with direct solver

u = A \ rhovec';
u = [0 u' 0]; % adding boundary conditions is possible in the Lagrange basis).

Get the values of $u_h$ and $u_h'$ at the vertices and at the midpoints:

u_000 = u(1:2:end-2);
u_025 = 3/8*u(1:2:end-2) + 3/4*u(2:2:end-1) - 1/8*u(3:2:end);
u_050 = u(2:2:end-1);
u_075 = -1/8*u(1:2:end-2) + 3/4*u(2:2:end-1) + 3/8*u(3:2:end);
u_100 = u(3:2:end);

Du_000 = 1./C(1:2:end).*( -3*u(1:2:end-2) + 4*u(2:2:end-1) - u(3:2:end) );
Du_025 = 1./C(1:2:end).*( -2*u(1:2:end-2) + 2*u(2:2:end-1) );
Du_050 = 1./C(1:2:end).*( -u(1:2:end-2) + u(3:2:end) );
Du_075 = 1./C(1:2:end).*( -2*u(2:2:end-1) + 2*u(3:2:end) );
Du_100 = 1./C(1:2:end).*( u(1:2:end-2) - 4*u(2:2:end-1) + 3*u(3:2:end) );

u_plot = u;

Du_plot = zeros(1,N+2);
Du_plot(1) = Du_000(1);
for i=1:M-1
    Du_plot(2*i) = Du_050(i);
    Du_plot(2*i+1) = Du_000(i+1);
end
Du_plot(2*M) = Du_050(M);
Du_plot(2*M+1) = Du_100(M);