# 1D Heat Conduction using explicit Finite Difference Method

860 views (last 30 days)
Derek Shaw on 15 Dec 2016
Answered: youssef aider on 12 Feb 2019
Hello I am trying to write a program to plot the temperature distribution in a insulated rod using the explicit Finite Central Difference Method and 1D Heat equation. The rod is heated on one end at 400k and exposed to ambient temperature on the right end at 300k. I am using a time of 1s, 11 grid points and a .002s time step. When I plot it gives me a crazy curve which isn't right. I think I am messing up my initial and boundary conditions. Here is my code.
L=1;
t=1;
k=.001;
n=11;
nt=500;
dx=L/n;
dt=.002;
alpha=k*dt/dx^2;
T0(1)=400;
for j=1:nt
for i=2:n
T1(i)=T0(i)+alpha*(T0(i+1)-2*T0(i)+T0(i-1));
end
T0=T1;
end
plot(x,T1)

KSSV on 15 Dec 2016
First place, it is not giving any curve..there is a error in your code. Please recheck your code once.
Derek Shaw on 15 Dec 2016
I am not sure I understand correctly. In the above I wrote this equation to be iterated With Boundary conditions and Initial Conditions with T0=400k and TL=Ti=300k
I am not sure how to set these boundary conditions in the code. Or if there is a curve I need to derive before doing the iterations.

michio on 15 Dec 2016
Edited: michio on 15 Dec 2016
It seems your initial condition and boundary conditon (x = L) are missing in the code. Try
L=1;
t=1;
k=.001;
n=11;
nt=500;
dx=L/n;
dt=.002;
alpha=k*dt/dx^2;
T0=400*ones(1,n);
T1=300*ones(1,n);
T0(1) = 300;
T0(end) = 300;
for j=1:nt
for i=2:n-1
T1(i)=T0(i)+alpha*(T0(i+1)-2*T0(i)+T0(i-1));
end
T0=T1;
end
plot(T1)

Show 1 older comment
michio on 15 Dec 2016
Glad to know that you figure things out. One additional tip is vectorization instead of for-loop, ie.
for i=2:n-1
T1(i) = T0(i) + alpha*(T0(i+1)-2*T0(i)+T0(i-1));
end
is equivalent to
T1(2:n-1) = T0(2:n-1) + alpha*(T0(3:n)-2*T0(2:n-1)+T0(1:n-2));
The later runs much faster.
Kunpeng Liao on 13 Oct 2017
Are the time step and grid spacing missing? I guess it would be:
if true
T1(2:n-1) = T0(2:n-1) + alpha*dt/dx^2*(T0(3:n)-2*T0(2:n-1)+T0(1:n-2));
end
Torsten on 16 Oct 2017
Look at the definition of "alpha" in the code ...
Best wishes
Torsten.

youssef aider on 12 Feb 2019
here is one, you can just change the boundaries
clear
clc
clf
% domain descritization
alpha = 0.05;
xmin = 0;
xmax = 0.2;
N = 100;
dx = (xmax-xmin)/(N-1);
x = xmin:dx:xmax;
dt = 4.0812E-5;
tmax = 1;
t = 0:dt:tmax;
% problem initialization
phi0 = ones(1,N)*300;
phiL = 230;
phiR = phiL;
% solving the problem
r = alpha*dt/(dx^2) % for stability, must be 0.5 or less
for j = 2:length(t) % for time steps
phi = phi0;
for i = 1:N % for space steps
if i == 1 || i == N
phi(i) = phiL;
else
phi(i) = phi(i)+r*(phi(i+1)-2*phi(i)+phi(i-1));
end
end
phi0 = phi;
plot(x,phi0)
shg
pause(0.05)
end

okncyln on 22 Dec 2017
implicit?