See discussions, stats, and author profiles for this publication at: https://www.researchgate.
net/publication/277009199
Matlab code to solve heat equation and notes
Research · May 2015
DOI: 10.13140/RG.2.1.5002.6400
CITATIONS READS
0 3,756
1 author:
Sabahat Qasim Khan
Riphah International University
13 PUBLICATIONS 11 CITATIONS
SEE PROFILE
All content following this page was uploaded by Sabahat Qasim Khan on 21 May 2015.
The user has requested enhancement of the downloaded file.
5/21/2015 prog1
clc;
clear all;
close all;
% problem u_t = beat u_xx
% beta 4/pi^2
% exact solution
%
% u(x,t) = exp(-t) * sin(pi/2*x) + exp(-t/4)* sin(pi/4*x)
%
u =@(x,t) exp(-t).*sin(pi/2.*x)+exp(-t/4).*sin(pi/4.*x);
t0 = 0;
tn = 0.08;
x0 = 0;
xn = 4;
dx = 0.2;
dt = 0.04;
x = (x0:dx:xn)';
t = (t0:dt:tn)';
beta = 4/pi^2;
s = beta*dt/dx^2;
nx = (xn-x0)/dx;
nt = (tn-t0)/dt;
% we start indexing from zero
nx_int = nx-1; % number of interior points in spatial dim
nt_int = nt-1; % number of interior points in temporal dim
% construction of tridiagonal matrix
indx1 = ones(1,nx-2);
indx2 = ones(1,nx-1);
a1 = -s*indx1;
a2 = (1+2*s)*indx2;
a3 = -s*indx1;
A = full(gallery('tridiag',a1,a2,a3));
% boundary conditions
% U(0,k) = 0
% U(20,k) = 0
% initial condition
int_cond = @(x) sin(pi/4*x).*(1+2*cos(pi/4*x));
+
% A U(k 1) = U(k) + b
1/5
5/21/2015 prog1
U=zeros(nx+1,nt+1);
U_exact=zeros(nx+1,nt+1);
U(1,:)=0;
U(nx+1,:)=0;
U(:,1) = int_cond(x);
for k=1:nt_int+1
U(2:nx_int+1,k+1) = A\U(2:nx_int+1,k);
end
% exact solution construxtion
for k=1:nt+1
U_exact(1:nx+1,k) =u(x,t(k));
end
Numerical_sol = U
Analytical_sol = U_exact
figure(1)
mesh(t,x,U)
xlabel('Temporal dim')
ylabel('Spatial dim')
zlabel('Numerical sol')
figure(2)
mesh(t,x,U_exact)
xlabel('Temporal dim')
ylabel('Spatial dim')
zlabel('Analytical sol')
figure(3)
mesh(t,x,U_exact-U)
xlabel('Temporal dim')
ylabel('Spatial dim')
zlabel('Error in sol')
Numerical_sol =
0 0 0
0.4655 0.4521 0.4392
0.8968 0.8713 0.8467
1.2630 1.2277 1.1935
1.5388 1.4967 1.4561
1.7071 1.6620 1.6183
1.7601 1.7158 1.6730
1.7000 1.6603 1.6219
fi 2/5
5/21/2015 prog1
1.5388 1.5070 1.4761
1.2967 1.2752 1.2542
1.0000 0.9901 0.9803
0.6787 0.6807 0.6824
0.3633 0.3763 0.3886
0.0820 0.1041 0.1250
-0.1420 -0.1137 -0.0868
-0.2929 -0.2617 -0.2319
-0.3633 -0.3328 -0.3036
-0.3550 -0.3286 -0.3034
-0.2788 -0.2594 -0.2408
-0.1526 -0.1423 -0.1325
-0.0000 0 0
Analytical_sol =
0 0 0
0.4655 0.4518 0.4386
0.8968 0.8707 0.8455
1.2630 1.2268 1.1918
1.5388 1.4957 1.4541
1.7071 1.6609 1.6162
1.7601 1.7147 1.6709
1.7000 1.6594 1.6202
1.5388 1.5063 1.4748
1.2967 1.2748 1.2534
1.0000 0.9900 0.9802
0.6787 0.6810 0.6829
0.3633 0.3769 0.3896
0.0820 0.1048 0.1265
-0.1420 -0.1128 -0.0849
-0.2929 -0.2607 -0.2300
-0.3633 -0.3318 -0.3018
-0.3550 -0.3278 -0.3018
-0.2788 -0.2588 -0.2397
-0.1526 -0.1420 -0.1319
-0.0000 -0.0000 -0.0000
3/5
5/21/2015 prog1
4/5
5/21/2015 prog1
Published with MATLAB® R2014a
5/5
View publication stats
