This is MATLAB code for (ut=uxx ) on {0,1} where

u(x,t)=sin(2pi*x)exp[-(2 pi).^2]*t,

I.C: u(x,t) = sin(2pi*x);

B.C: u(0,t) = u(1,t) =0;

L=1; n=1;

t_final =0.1;

N = 32; %elements

dx =1/N;

dt = 0.01;

x = 0:dx:L;

t=0:dt:t_final;

for i=1:length(t)

y=sin(2*pi*x)*(exp((-(2*pi).^2)*t(i)));

y(1)=0;

y(end)=0;

figure(1)

plot(x,y)

hold on

end

%% Backward Euler for FEM solution

% L = 2;

xx =0:dx:L;

nx = length(xx); %Nodes

M = spdiags([(dx/6)*ones(nx,1), (2*dx/3)*ones(nx,1), (dx/6)*ones(nx,1)],[-1, 0,1],nx,nx) ; %% Mass matrix

A = spdiags([(-1/dx)*ones(nx,1), (2/dx)*ones(nx,1), (-1/dx)*ones(nx,1)],[-1, 0, 1],nx,nx) ;%% Stiffness matrix

delta_t = 0.001;

tspan = 0:delta_t:t_final;

yy = zeros(1,nx);

for nt=1:length(tspan)

yy = ( sin(2*pi*xx/L)*exp((-(2*pi).^2)*tspan(nt)))*(M/(M+delta_t*A));

yy(1)= sin(2*pi*xx(1));

yy(end)=0;

figure(2)

plot(xx,yy)

hold on

end

I used formula for error, error = norm(exact solution- FEM solution)

but discretization error does not converges to zero.

I obtained following errors by changing no of elements;

N Error

4 0.00100000

8 0.00074918

16 0.0007419

32 0.00073473

64 0.00073303

128 0.00073274

256 0.00073274

512 0.00073470

Similar questions and discussions