Galerkin FEM solver for 1D elliptic problem

This is a solver for the Dirichlet boundary-value problem

$$ -\phi''(x) + \phi(x) = f(x)\,,\qquad x\in[a,b]\subset \mathbf R. $$

$$ \phi(a) = \phi(b) = 0\,,$$

where $f\in C^\infty([a,b])$ is the given data. We search for numerical solutions $u_h\in V_h$, where $V_h \subset V$ is a subspace of $V=H_0^1((a,b))$. $V_h$ consists of piece-wise polynomials of first and second order, respectively. LAGRANGE polynomials or SPLINES are used as basis functions.

Contents

clear all
close all

scrsz = get(0,'screensize');

a = 0; % left domain boundary
b = 2; % right domain boundary
M = 16; % number of cells (elements)

L = b-a;

syms z % define a symbolic variable z
u_ex(z) = -sin(2*pi/L*z)*exp(z); % fabricated exact solution
%u_ex(z) = (z-a)*(z-b)*(z-(a+b)/2)^4 % fabricated polynomial solution
Du_ex(z) = diff(u_ex(z),1); % derivative
rho(z) = -diff(u_ex(z),2) + u_ex(z); % fabricated right-hand side

x = zeros(1,M+1); % vector of vertices

%switches:
%swi_grid = 'ana';
swi_grid = 'num';

%swi_meth = 'Lag2';
swi_meth = 'spl2'; % splines with p-1 regularity.

swi_rhs = 'exact'; % must use this to compute the right-hand side
%swi_rhs = 'basis';

save_err = true;
swi_showerr = true;

Create a non-uniform grid

if swi_grid=='ana'

    %H(z) = 1/(exp(4*z)+2) % spacing function
    H(z) = sin(2*pi*z) + 1.2; % spacing function
    cst = M/int(1/H,a,b)
    prim(z) = cst*int(1/H);
    prim(z) = prim(z) - prim(a); % set the correct constant in the primitive
    vpa(prim(a))
    vpa(prim(b))

    for j=1:M-1
        vertex = solve(prim(z) - j == 0,z);
        x(j+1) = vpa(vertex);
    end
    x(M+1) = b;

    z = linspace(a,b,100);
    plot(x,zeros(size(x)),'bx',z,vpa(H(z))); % visual check of grid
    xlim([a b])

elseif swi_grid=='num'

    %H = '0.05 + 0.1*(y<0.3) + 0.05*(y>0.5)'
    H = '0.1'
    count = 0;
    yvec = a;
    y = a;
    while y<b
        count = count + 1;
        del = eval(H);
        yvec = [yvec  y + del];
        y = yvec(end);
    end
    M = count
    cst = (b-a)/(y-a);
    xvec=a;
    for i = 1:M
        y = yvec(i);
        del = cst*eval(H);
        xvec = [xvec  xvec(end) + del];
    end
    x = xvec;

    clear y
    y = linspace(a,b,100);
    plot(x,zeros(size(x)),'bx',y,eval(H)); % visual check of grid
    title('Grid spacing')
    xlim([a b])

end
H =

0.1


M =

    20

Choose the approximation space

if swi_meth=='Lag1'

    Lagrange1

elseif swi_meth=='Lag2'

    Lagrange2

elseif swi_meth=='spl1'

    spline1

elseif swi_meth=='spl2'

    spline2_ref

end

Diagnostics

nod_ab = [a nodes b]; % vector of nodes (vertices and midpoints)
z = linspace(a,b,200); % set the symbolic variable
figure('position',[200 200 scrsz(3)*0.7 scrsz(3)*0.4])
subplot(1,2,1)
plot(z,eval(u_ex),'r',nod_ab,u,'ob')
legend('exact','numerics')
set(gca,'fontsize',16)
subplot(1,2,2)
plot(z,eval(Du_ex),'r',nod_ab,Du_plot,'ob')
legend('exact','numerics')
set(gca,'fontsize',16)

dxvec = x(2:end) - x(1:end-1);
dx = max(dxvec) % grid parameter
q25 = 	3/4*x(1:end-1) + 1/4*x(2:end);
mids = 	1/2*x(1:end-1) + 1/2*x(2:end);
q75 = 	1/4*x(1:end-1) + 3/4*x(2:end);

% L2 error:
L2 = 0;
for i=1:M

    z = x(i);
    e_000 = eval(u_ex) - u_000(i);
    z = q25(i);
    e_025 = eval(u_ex) - u_025(i);
    z = mids(i);
    e_050 = eval(u_ex) - u_050(i);
    z = q75(i);
    e_075 = eval(u_ex) - u_075(i);
    z = x(i+1);
    e_100 = eval(u_ex) - u_100(i);

    %L2 = L2 + dxvec(i)/2*( abs(e_000)^2 + abs(e_100)^2 ); % Trapezoid rule (order 2).
    %L2 = L2 + dxvec(i)/6*( abs(e_000)^2 + 4*abs(e_050)^2 + abs(e_100)^2 ); % Simpson rule (order 4).
    L2 = L2 + dxvec(i)/90*( 7*abs(e_000)^2 + 32*abs(e_025)^2 + 12*abs(e_050)^2 + 32*abs(e_075)^2 + 7*abs(e_100)^2 ); % Boole rule (order 6).

end
L2 = sqrt(L2)

% H1 error:
H1 = 0;
for i=1:M

    z = x(i);
    e_000 = eval(Du_ex) - Du_000(i);
    z = q25(i);
    e_025 = eval(Du_ex) - Du_025(i);
    z = mids(i);
    e_050 = eval(Du_ex) - Du_050(i);
    z = q75(i);
    e_075 = eval(Du_ex) - Du_075(i);
    z = x(i+1);
    e_100 = eval(Du_ex) - Du_100(i);

    %H1 = H1 + dxvec(i)/2*( abs(e_000)^2 + abs(e_100)^2 ); % Trapezoid rule (order 2).
    %H1 = H1 + dxvec(i)/6*( abs(e_000)^2 + 4*abs(e_050)^2 + abs(e_100)^2 ); % Simpson rule (order 4).
    H1 = H1 + dxvec(i)/90*( 7*abs(e_000)^2 + 32*abs(e_025)^2 + 12*abs(e_050)^2 + 32*abs(e_075)^2 + 7*abs(e_100)^2 ); % Boole rule (order 6).

end
H1 = sqrt(L2^2 + H1)

if save_err
    save(['errors/err_' num2str(M) '_' swi_meth '_' swi_rhs],'L2','H1','dx'); % save errors to file
end

if swi_showerr

    MMvec = [5 11 20 40 81 161 320];
    for i=1:7
        err = load(['errors/err_' num2str(MMvec(i)) '_' swi_meth '_' swi_rhs]);
        hvec(i) = err.dx; % vector of grid parameters h
        L2vec(i) = err.L2; % vector of L2 norms
        H1vec(i) = err.H1; % vector of H1 norms
    end

    figure('position',[300 300 scrsz(3)*0.7 scrsz(3)*0.4])
    if swi_meth=='Lag1' | swi_meth=='spl1'

        loglog(hvec,L2vec,'ob',hvec,H1vec,'xg',hvec(3:end),H1vec(3)*hvec(3:end)/hvec(3),'-k', ...
            hvec(3:end),L2vec(3)*(hvec(3:end)/hvec(3)).^2,'--k','linewidth',2,'markersize',14)
        legend('L^2 error','H^1 error','h^1','h^2','location','northwest')

    elseif swi_meth=='Lag2' | swi_meth=='spl2'

        loglog(hvec,L2vec,'ob',hvec,H1vec,'xg',hvec(3:end),H1vec(3)*(hvec(3:end)/hvec(3)).^2,'-k', ...
            hvec(3:end),L2vec(3)*(hvec(3:end)/hvec(3)).^3,'--k','linewidth',2,'markersize',14)
        legend('L^2 error','H^1 error','h^2','h^3','location','northwest')

    end
    xlabel('h')
    ylabel('errors')
    title(swi_meth)
    set(gca,'fontsize',16)

end
dx =

    0.1000


L2 =

   7.5015e-04


H1 =

    0.0416