889 views (last 30 days)

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)

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)

http://jp.mathworks.com/matlabcentral/fileexchange/59916-simple-heat-equation-solver could make a good example for you.

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

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

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!
## 2 Comments

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/316950-1d-heat-conduction-using-explicit-finite-difference-method#comment_413805

⋮## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/316950-1d-heat-conduction-using-explicit-finite-difference-method#comment_413805

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/316950-1d-heat-conduction-using-explicit-finite-difference-method#comment_413833

⋮## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/316950-1d-heat-conduction-using-explicit-finite-difference-method#comment_413833

Sign in to comment.