below are the codes to solve heat transfer using implicit and explicit method but my implicit method is showing huge error, what is wrong on the implicit method

clc

clear

close all

nx = 10;

nt = 50;

a = 0;

b = 1;

t0 = 0;

tf = 0.2;

dx = (b-a)/(nx-1);

dt = (tf-t0)/(nt-1);

x = a:dx:b;

t = t0:dt:tf;

s = dt/dx^2

% Mesh Figure

figure()

mesh(x, t, zeros(nt, nx), 'EdgeColor', 'k', 'FaceAlpha', 0.5);

xlabel('X');

ylabel('t');

zlabel('Solution');

title('Mesh');

%% Analytical solution

UA = zeros (nx,nt);

for j = 1:nt

for i = 1:nx

UA(i,j) = sin(pi*x(i))*exp(-pi^2*t(j));

end

end

figure()

contourf (UA,200,'linecolor','non')

xlabel('X')

ylabel('t')

title('Analytical Solution')

colormap(jet(256))

colorbar

caxis([0,1])

%% Numerical Solution (Explicit Scheme)

UN = zeros (nx,nt);

% initial condition

UN(:,1) = sin(pi*x);

for j = 1:nt-1

for i = 2:nx-1

UN(i,j+1) = s*UN(i-1,j)+(1-2*s)*UN(i,j)+s*UN(i+1,j);

end

end

figure()

contourf (UN,200,'linecolor','non')

xlabel('X')

ylabel('t')

title('Numerical Solution (Explicit Solution)')

colormap(jet(256))

colorbar

caxis([0,1])

%% Numerical Solution (Implicit Scheme)

UT = zeros (nx,nt);

% initial condition

UT(:,1) = sin(pi*x);

for j = 1:nt-1

for i = 2:nx-1

UT(i-1,j) = -s*UT(i-1,j)+(1+2*s)*UT(i,j)-s*UT(i+1,j);

end

end

figure()

contourf (UT,200,'linecolor','non')

xlabel('X')

ylabel('t')

title('Numerical Solution (Implicit Solution)')

colormap(jet(256))

colorbar

caxis([0,1])

%% Explicit Error

E = abs(UA-UN);

figure()

contourf (E,200,'linecolor','non')

xlabel('X')

ylabel('t')

title('Explicit Error')

colormap(jet(256))

colorbar

%% Implicit Error

E = abs(UA-UT);

figure()

contourf (E,200,'linecolor','non')

xlabel('X')

ylabel('t')

title('Implicit Error')

colormap(jet(256))

colorbar

Answer 1

[Torsten]

Implicit scheme means:

(ut(i,j+1)-ut(i,j))/dt = (ut(i-1,j+1)-2*ut(i,j+1)+ut(i+1,j+1))/dx^2 

or

ut(i,j+1)/dt - (ut(i-1,j+1)-2*ut(i,j+1)+ut(i+1,j+1))/dx^2 = ut(i,j)/dt

or

-ut(i-1,j+1)/dx^2 + (2/dx^2+1/dt)*ut(i,j+1) - ut(i+1,j+1)/dx^2 = ut(i,j)/dt  (2 <= i <= nx-1)

This is a tridiagonal linear system of equations to be solved for the vector ut(:,j+1) in each time step.