# MATLAB Finite Difference Method using Matlab

#### mpm

I am curious to know if anyone has a program that will solve for 2-D Transient finite difference.

I have a project in a heat transfer class and I am supposed to use Matlab to solve for this. However, when I took the class to learn Matlab, the professor was terrible and didnt teach much at all.

Any help would be appreciated.

Also, can I possibly use like a separation of variables process in Matlab and that count as finite difference?

I know there is always more than one way to do things in Matlab and still get the same thing but I dont know that those two things can qualify as the same thing.

Related MATLAB, Maple, Mathematica, LaTeX, Etc News on Phys.org

#### Clausius2

Gold Member
mpm said:
I am curious to know if anyone has a program that will solve for 2-D Transient finite difference.

I have a project in a heat transfer class and I am supposed to use Matlab to solve for this. However, when I took the class to learn Matlab, the professor was terrible and didnt teach much at all.

Any help would be appreciated.

Also, can I possibly use like a separation of variables process in Matlab and that count as finite difference?

I know there is always more than one way to do things in Matlab and still get the same thing but I dont know that those two things can qualify as the same thing.
I think you should describe how is your problem and how you plan to solve it, and maybe we could advice you.

#### mpm

I guess a more specific problem would help. I dont want to type it all out so here is a link.

http://www.me.ttu.edu/files/pantoya/S05 - Matlab Finite Diff DP2 asgnmnt.pdf

The professor gave us a simliar example and he wrote some code for it. I think I might be able to manipulate it into fitting my problem but Im not exactly sure because the conditions are different on each side of the block.

Here is a link to his code.

http://www.me.ttu.edu/files/pantoya/ftcs_block_cavity.m

Also here is a link to his problem he solved.

http://www.me.ttu.edu/files/pantoya/Block with cavity example.doc

Can maybe someone tell me what part needs to be changed to alter the conditions I need? If I know which part needs to be changed maybe I can figure it out. I think I know which parts to remove for the hole in the block of his problem.

Any help is appreciated.

Holy cow that's some messy code. Do yourself a favor and represent the entire surface with a single vector. Then your lapacian becomes a tridiagonal matrix, and the boundary conditions are much easier to enforce that way (you don't need a million if elseif statements).

But if you don't know how to do that, it looks like you just pretty much use the code that was given, but you need to alter the boundary conditions to match your problem. The problem you are given has much simpler conditions than the one that the code is given for, so it shouldn't be that hard.

Then it looks likes he's using a Forward Euler scheme for the time integration, and Forward Euler often blows up when it's not supposed to. There's a parameter involving delta T and the coefficients of the equation that controls the stability of Forward Euler. You want to place this value in the unit disc on the left half of the real axis to make Forward Euler stable.

#### Clausius2

Gold Member
Holy cow that's some messy code.
The fact is that professor has shown not to have any idea about elegance at programming. That problem can be solved in 60 lines or so.

Last edited:

I don't see the need for matrix inversion or Gauss-Siedel because he's already solving explicitly for the values at the new time. The heat equation explicitly in time is $$u^{n+1} = u^n + Au^n$$, and he's given u(x,y,0). I was just suggesting that by making A a matrix and u a vector he could reduce all the stepping around the surface to one matrix-vector product.

Now if he were solving Poisson's equation, $$Au=f$$, then he would need matrix inversion or Gauss-Siedel. But even then, I would recommend Krylov subspace methods over both and still go with matrix-vector form.

#### Clausius2

Gold Member
I don't see the need for matrix inversion or Gauss-Siedel because he's already solving explicitly for the values at the new time.
What inversion are you talking about???

Take a look at my post, Have I ever mentioned inverting anywhere?? :tongue2: :rofl:

#### mpm

Ok. I have edited my professor's code and taken out the cavity. Or at least I think I did. Im sure I need to change some boundaries but I am not quite sure what I need to change yet. Since you say it is messy, maybe there is an easier way for me to do this.

I noticed you mention Laplace and triadiagonal. I know of a code that does this I think but I am not sure how I need to manipulate it.

I apologize for having to click the download button. This was just on the Matlab file sharing thing and I think this may be what you speak of.

If this is what you speak of, maybe you can help me manipulate it into working.

I appreciate everyones help here. I feel a little resentful in the fact that I was jipped of learning Matlab because my professor was angry about getting a sophomore level class when he usually teaches senior level classes.

Thanks to all.

Last edited by a moderator:

I'm afraid I have only created more confusion here. Forget everything I said about matrices. You can do it the way your professor did and it will work fine. Why don't you show us what you have changed his code to and we can tell you if you're on the right track.

#### mpm

Here is my code that I have. Sorry I cut and pasted it on here, but I dont have a server to put the file on. Im getting it to run, but there are problems somewhere because my plot doesnt come out right.

Any help with my errors would be appreciated. I probably need to change a coordinate with an i or j somewhere and maybe some other things.

Thanks.

mpm

clear
clc

%define material properties
k = 0.401; %conduction coefficient units = W/cm*K
Cp = 385; %specific heat units = J/kg*K
rho = 0.008933; %density units = kg/cm^3

%define geometry
Lx = 50; Ly = 20;

dx = 1; %differential x length of cells
Nx = Lx/dx; %number of x nodes
dy = 1; %differential y length of cells
Ny = Ly/dy; %Number of Y nodes

Cx = round(0.4/dx); Cy = round(0.25/dy);

tf = 60; %units of seconds
dt = 10; %seconds

Box_x = 0; %define lower left corner of box
Box_y = 0;

%Initialize mesh
T(:,:,1)=25*ones(Nx,Ny);
T(1,:,1)=100;

time(1) = 0;
time(2) = dt;
t=1;
while time(t+1)<=tf %index time step increments
for i=1:Nx; %index in x direction
for j=1:Ny; %index in y direction
if i==1; %identify left boundary
T(i,j,t+1)=100;
elseif i>=2 & i<=(Cx-1);
if j==1; %lower boundary
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/(dx) + dx*(T(i,j+1,t)-T(i,j,t))/dy + dy*(T(i+1,j,t)-T(i,j,t))/dx)+T(i,j,t);
elseif j==Ny; %upper boundary
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
else %all other vertical intermediate nodes on left boundary
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dx*(T(i,j+1,t)-T(i,j,t))/dy + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
end
T(i,1,t+1)=T(i,2,t+1);
if j==1; %lower boundary
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/(dx) + dx*(T(i,j+1,t)-T(i,j,t))/dy + dy*(T(i+1,j,t)-T(i,j,t))/dx)+T(i,j,t);
elseif j==Ny; %top boundary
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
end
if j==1; %bottom
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/(dx) + dx*(T(i,j+1,t)-T(i,j,t))/dy + dy*(T(i+1,j,t)-T(i,j,t))/dx)+T(i,j,t);
elseif j==Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
end
elseif i==Nx;
if j==1;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/(dx) + dx*(T(i,j+1,t)-T(i,j,t))/dy + -10)+T(i,j,t);
elseif j>1&j<Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx - 10 + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
elseif j==Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + -10 + dx*(T(i,j-1,t)-T(i,j,t))/dy)+T(i,j,t);
end
end
end
end
t = t+1;
time(t+1) = time(t)+dt;
end

x(1) = Box_x;
y(1) = Box_y;
for i=2:Nx;
x(i) = x(i-1)+dx;
end
for j=2:Ny;
y(j) = y(j-1)+dy;
end

figure(1); clf;
Tfinal = T(:,:,end);
mesh(y,x,Tfinal)
title(strcat('Solution at time = ',num2str(time(end))));
xlabel('y')
ylabel('x')

true_time = time(1:end-1);
clear time;
time = true_time;

figure(2); clf;
center_y = round(Ly/2);
for e = 1:length(time)
T1(e) = T(2,center_y,e)
T4(e) = T(end,center_y,e)
end
plot(time,T1,'k',time,T4,'b')
title('Temp history of specific center line nodes')
xlabel('time')
ylabel('T ^oC')

Ok, first of all it looks like Cx and Cy were coordinates for the cavity in the other example, so those shouldn't show up anywhere in your code. You erased most of them but still have a crucial elseif depending on one. In fact, if you just change that Cx to Nx, I think you will find that you get a very nice looking solution. Unfortunately it's not the solution to the problem you want to be solving.

Second, when you cut out one elseif statement about i==Cx you also need to kill the stuff below it. As it is now your first elseif statement contains repeated equations for the top and bottom boundaries.

Most importantly though, rather than just changing the geometry, you also need to change the actual boundary conditions. Your professer's problem held the left side at a constant temperature and had a flux leaving the right side. Your problem appears to be heating the left side and holding the right side constant. So your staments regarding the right side, i==Nx, should just set T = 100. Now since you are heating the left side of the box, i == 1, you want to be increasing the temperature at that boundary by some amount at each time step. You can figure out exactly what to put by thinking about what 2000 W means.

Hopefully this points you in the right direction.

#### mpm

I just want to make sure I get a couple of things straight.

In the Original code, CLx and Cly were used for the cavity dimensions. Are you sure Cx and Cy are for the cavity as well?

If that is the case, do I need to remove any line of code that involves Cx and Cy or do I need to change that to Nx and Ny.

As for the elseif statements, do I need to remove all of that stuff involved with elseif and Cx? Or do I need to remove all elseif statements? Just want to be clear on this.

Once I get to this part, I may can handle the rest about the boundary conditions. I will just have to see. I just want to make sure all of my geometry coding is correct before I mess with boundary conditions.

One last question, if I mess dx and dt, should I keep those relatively small (ie.e < 1cm) or larger (>1cm)? It seems like there might be a stability problem if I go too small. Is that the problem or is it just the coding causing the problem?

I think Clx and Cly are the size of the cavity, but Cx and Cy are where it starts or stops on the domain. I admit that I didn't look at it long enough to figure out what they were referring to, but since your domain is just a big rectangle, the only boundaries you should need are Nx and Ny.

You are working your way from bottom to top and then left to right (assuming (1,1) is the bottom left corner), so I think the general outline should be something like this.

if x == left wall,
if y = bottom, do something - bottom left corner​
elseif y = middle, do something - left side​
elseif y = top, do something - top left corner​
end​
elseif x is in the middle
if y == bottom, do something - bottom​
elseif y = middle, do something -- This will handle most of the points​
elseif y = top, do something - top​
end​
elseif x is on the right boundary
if y = bottom do something - bottom right corner​
elseif y = middle do something - right side​
elseif y = top do something - top right corner​
end​
end

Those statements should only depend on Nx, Ny, 1, and 1.

As for timesteps, this is just a wild guess from my limited knowledge of semi-discrete Fourier analysis, but the fact that you have a Laplacian makes me say set $$\Delta T < \Delta X^2$$ to keep it stable. You'll figure it out when you do the stability analysis.

Last edited:

#### mpm

Thanks for your time and help on this. Hopefully I can figure it out from what you have told me. At least I have a direction to go towards. Hopefully I am not annoying you with my questions.

I have until Thursday to finish this so hopefully I can figure it out.

#### mpm

%ftcs_block_cavity.m
%
%forward in time centered in space
%based on heat transfer energy balance at all nodes
%
%two dimensional finite difference solution to a rectangular block with
%a rectangular crack or cavity

clear
clc

%define material properties
k = 0.401; %conduction coefficient units = W/cm*K
Cp = 385; %specific heat units = J/kg*K
rho = 0.008933; %density units = kg/cm^3

%define geometry
Lx = 50; Ly = 20;

dx = 1; %differential x length of cells
Nx = Lx/dx; %number of x nodes
dy = 1; %differential y length of cells
Ny = Ly/dy; %Number of Y nodes

tf = 60; %units of seconds
dt = 1; %seconds

Box_x = 0; %define lower left corner of box
Box_y = 0;

%Initialize mesh
T(:,:,1)=25*ones(Nx,Ny);
T(1,:,1)=100;

time(1) = 0;
time(2) = dt;
t=1;
while time(t+1)<=tf %index time step increments
for i=1:Nx; %index in x direction
for j=1:Ny; %index in y direction
if i==1; %identify left boundary
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/(dx) + dx*(T(i,j+1,t)-T(i,j,t))/dy + 2000)+T(i,j,t);
elseif j>1&j<Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + 2000 + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
elseif j==Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + 2000 + dx*(T(i,j-1,t)-T(i,j,t))/dy)+T(i,j,t);
end

if j==1; %lower boundary
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/(dx) + dx*(T(i,j+1,t)-T(i,j,t))/dy + dy*(T(i+1,j,t)-T(i,j,t))/dx)+T(i,j,t);
elseif j==Ny; %upper boundary
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
else %all other vertical intermediate nodes on left boundary
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dx*(T(i,j+1,t)-T(i,j,t))/dy + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
end

T(i,1,t+1)=T(i,2,t+1);
if j==1; %lower boundary
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/(dx) + dx*(T(i,j+1,t)-T(i,j,t))/dy + dy*(T(i+1,j,t)-T(i,j,t))/dx)+T(i,j,t);
elseif j==Ny; %top boundary
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
end

if j==1; %bottom
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/(dx) + dx*(T(i,j+1,t)-T(i,j,t))/dy + dy*(T(i+1,j,t)-T(i,j,t))/dx)+T(i,j,t);
elseif j==Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
end

if j==1;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/(dx) + dx*(T(i,j+1,t)-T(i,j,t))/dy + dy*(T(i+1,j,t)-T(i,j,t))/dx)+T(i,j,t);
elseif j==Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
end

if j==1;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/(dx) + dx*(T(i,j+1,t)-T(i,j,t))/dy + dy*(T(i+1,j,t)-T(i,j,t))/dx)+T(i,j,t);
elseif j>1&j<Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/(dx) + dx*(T(i,j+1,t)-T(i,j,t))/dy + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/dy)+T(i,j,t);
elseif j==Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
end

if j==1;
elseif i==Nx;
T(i,j,t+1)=100;
end
end
end
end
t = t+1;
time(t+1) = time(t)+dt;
end

x(1) = Box_x;
y(1) = Box_y;
for i=2:Nx;
x(i) = x(i-1)+dx;
end
for j=2:Ny;
y(j) = y(j-1)+dy;
end

figure(1); clf;
Tfinal = T(:,:,end);
mesh(y,x,Tfinal)
title(strcat('Solution at time = ',num2str(time(end))));
xlabel('y')
ylabel('x')

true_time = time(1:end-1);
clear time;
time = true_time;

figure(2); clf;
center_y = round(Ly/2);
for e = 1:length(time)
T1(e) = T(2,center_y,e);
T4(e) = T(end,center_y,e);
end
plot(time,T1,'k',time,T4,'b')
title('Temp history of specific center line nodes')
xlabel('time')
ylabel('T ^oC')

#### mpm

Here is my latest code

while time(t+1)<=tf %index time step increments
for i=1:Nx; %index in x direction
for j=1:Ny; %index in y direction

if i==1; %identify left boundary
if j==1;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/(dx) + dx*(T(i,j+1,t)-T(i,j,t))/dy + 2000)+T(i,j,t);
elseif j>1&j<Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + 2000 + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
elseif j==Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + 2000 + dx*(T(i,j-1,t)-T(i,j,t))/dy)+T(i,j,t);

elseif i>1 & i<Nx; %identify middle boundary
if j==1;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
elseif j>1&j<Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
elseif j==Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);

elseif i ==Nx; %identify right boundary
if j==1;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
elseif j>1 & j<Ny;
T(i,j,t+1) = [100];
elseif j==Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);

end
end
end
end
end

end
t = t+1;
time(t+1) = time(t)+dt;
end

x(1) = Box_x;
y(1) = Box_y;
for i=2:Nx;
x(i) = x(i-1)+dx;
end
for j=2:Ny;
y(j) = y(j-1)+dy;
end

figure(1); clf;
Tfinal = T(:,:,end);
mesh(y,x,Tfinal)
title(strcat('Solution at time = ',num2str(time(end))));
xlabel('y')
ylabel('x')

true_time = time(1:end-1);
clear time;
time = true_time;

figure(2); clf;
center_y = round(Ly/2);
for e = 1:length(time)
T1(e) = T(2,center_y,e);
T4(e) = T(end,center_y,e);
end

#### mpm

Here is my last post of code.

%ftcs_block_cavity.m
%
%forward in time centered in space
%based on heat transfer energy balance at all nodes
%
%two dimensional finite difference solution to a rectangular block with
%a rectangular crack or cavity

clear
clc

%define material properties
k = 0.401; %conduction coefficient units = W/cm*K
Cp = 385; %specific heat units = J/kg*K
rho = 0.008933; %density units = kg/cm^3

%define geometry
Lx = 50; Ly = 20;

dx = 1; %differential x length of cells
Nx = Lx/dx; %number of x nodes
dy = 1; %differential y length of cells
Ny = Ly/dy; %Number of Y nodes

tf = 60; %units of seconds
dt = 1; %seconds

Box_x = 0; %define lower left corner of box
Box_y = 0;

%Initialize mesh
T(:,:,1)=25*ones(Nx,Ny);
T(1,:,1)=100;

time(1) = 0;
time(2) = dt;
t=1;
while time(t+1)<=tf %index time step increments
for i=1:Nx; %index in x direction
for j=1:Ny; %index in y direction

if i==1; end %identify left boundary
if j==1;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/(dx) + dx*(T(i,j+1,t)-T(i,j,t))/dy + 2000)+T(i,j,t);
elseif j>1&j<Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + 2000 + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
elseif j==Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + 2000 + dx*(T(i,j-1,t)-T(i,j,t))/dy)+T(i,j,t);
end
end

if i>1 & i<Nx; end %identify middle boundary
if j==1;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
elseif j>1&j<Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
elseif j==Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
end
end

if i==Nx; end %identify right boundary
if j==1;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
elseif j>1 & j<Ny;
T(i,j,t+1) = [100];
elseif j==Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
end
end
t = t+1;
time(t+1) = time(t)+dt;
end

x(1) = Box_x;
y(1) = Box_y;
for i=2:Nx;
x(i) = x(i-1)+dx;
end
for j=2:Ny;
y(j) = y(j-1)+dy;
end

figure(1); clf;
Tfinal = T(:,:,end);
mesh(y,x,Tfinal)
title(strcat('Solution at time = ',num2str(time(end))));
xlabel('y')
ylabel('x')

true_time = time(1:end-1);
clear time;
time = true_time;

figure(2); clf;
center_y = round(Ly/2);
for e = 1:length(time)
T1(e) = T(2,center_y,e);
T4(e) = T(end,center_y,e);
end
plot(time,T1,'k',time,T4,'b')
title('Temp history of specific center line nodes')
xlabel('time')
ylabel('T ^oC')

#### mpm

LeBrad, Here is my coding as of 3/23/05 1:06 AM Central Time.

I took your adivce and copied and pasted my professors T equations. I manipulated them into what I thought was right. 5 of the 9 equations work I think. They dont give me problems at least when I run them.

The other 4 I guess I didnt manipulate in the write way. They are marked as comments. if you have any suggestions feel free to give them.

Here is my code.

Thanks for everything.

%ftcs_block_cavity.m
%
%forward in time centered in space
%based on heat transfer energy balance at all nodes
%
%two dimensional finite difference solution to a rectangular block with
%a rectangular crack or cavity

clear
clc

%define material properties
k = 0.401; %conduction coefficient units = W/cm*K
Cp = 385; %specific heat units = J/kg*K
rho = 0.008933; %density units = kg/cm^3

%define geometry
Lx = 50; Ly = 20;

dx = 1; %differential x length of cells
Nx = Lx/dx; %number of x nodes
dy = 1; %differential y length of cells
Ny = Ly/dy; %Number of Y nodes

tf = 60; %units of seconds
dt = 1; %seconds

Box_x = 0; %define lower left corner of box
Box_y = 0;

%Initialize mesh
T(:,:,1)=25*ones(Nx,Ny);
T(1,:,1)=100;

time(1) = 0;
time(2) = dt;
t=1;
%while time(t+1)<=tf %index time step increments
for i=1:Nx; %index in x direction
for j=1:Ny; %index in y direction

if i==1; end %identify left boundary
if j==1; %bottom left corner
% T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i+1,j,t)-T(i,j,t) + T(i,j+1,t)-T(i,j,t) + 2000 +T(i,j,t);
elseif j>1&j<Ny; %left side
% T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i+1,j,t)-T(i,j,t) + 2000 + T(i,j-1,t)-T(i,j,t)+T(i,j,t);
elseif j==Ny; %top left corner
% T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i+1,j,t)-T(i,j,t) + 2000 + T(i,j-1,t)-T(i,j,t)+T(i,j,t);
end
end

if i>1 & i<Nx; end %identify middle boundary
if j==1; %bottom
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i-1,j,t)-T(i,j,t) + T(i,j+1,t)-T(i,j,t) + T(i+1,j,t)-T(i,j,t) + T(i,j,t);
elseif j>1&j<Ny; %middle
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i+1,j,t) + T(i-1,j,t) + T(m,j+1,t) + T(i,j-1,t) + T(i,j,t);
elseif j==Ny; %top
% T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i-1,j,t)-T(i,j,t) + T(i+1,j,t)-T(i,j,t) + T(i,j-1,t)-T(i,j,t) + T(i,j,t);
end
end

if i==Nx; end %identify right boundary
if j==1; %bottom right corner
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i-1,j,t)-T(i,j,t) + T(i,j+1,t)-T(i,j,t)+T(i,j,t);
elseif j>1 & j<Ny; % right side
T(i,j,t+1) = [100];
elseif j==Ny; %top right corner
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i-1,j,t)-T(i,j,t) + T(i,j-1,t)-T(i,j,t)+T(i,j,t);
end
end
t = t+1;
time(t+1) = time(t)+dt;
end

x(1) = Box_x;
y(1) = Box_y;
for i=2:Nx;
x(i) = x(i-1)+dx;
end
for j=2:Ny;
y(j) = y(j-1)+dy;
end
time = [0:dt:tf];
figure(1); clf;
Tfinal = T(:,:,end);
mesh(y,x,Tfinal)
title(strcat('Solution at time = ',num2str(time(end))));
xlabel('y')
ylabel('x')

true_time = time(1:end-1);
clear time;
time = true_time;

figure(2); clf;
center_y = round(Ly/2);
for e = 1:length(time)
T1(e) = T(2,center_y,e);
T4(e) = T(end,center_y,e);
end
plot(time,T1,'k',time,T4,'b')
title('Temp history of specific center line nodes')
xlabel('time')
ylabel('T ^oC')

You need to delete the END statments that appear like this
if i ==1; end;
You have 3 of those that shouldn't be there.

I fixed your time counting, and I changed the line near the top that says
T(Nx,:,1) = 100; That sets the right side to 100 to start. Before it was setting the left side.

You need another set of parentheses in your equations so your coefficients multiply all of the differences.

This is basically the heat equation
$$a*\frac{\partial T}{\partial t} = \nabla^2 T= \frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2}$$

It is saying that the change in temperature with respect to time depends on second spatial derivatives of the temperature. Since we don't know the function T, we can't take its derivatives. But we can approximate them as differences. The derivative of a line just tells you its slope. And the slope is just the rise over the run. So we can approximate a first derivative by saying $$\frac{\partial T}{\partial x} = \frac{T_{i+1} - T_i}{\Delta x}$$ where T(i+1) and T(i) are the temperatures at two nodes adjacent in the x direction.

To get a second derivative, you just take the derivative of two first derivative points, $$\frac{\partial^2 T}{\partialx x^2} = \frac{\frac{T_{i+1} - T_i}{\Delta x} - \frac{T_{i}-T_{i-1}}{\Delta x}}{\Delta x}$$

To get the formulas for the derivatives with respect to y, just replace the i's with j's.

So when you put it all together, your equation for the temperature at a node T(i,j) will be a function of $$T_{i+1,j}, T_{i,j}, T_{i-1,j}, T_{i,j+1}, T_{i,j-1}$$

This means that the demperature at a point depends on the temperatures of the 4 points around it at the previous step. Now, the heat equation is only definded in the middle of your rectange and ALONG the edges, but not normal to the edges. So any point in the middle of the rectangle will have a formula involving those 5 T() values.

When you're on the left edge, you can't use T(i-1,j) because there is no node there. You can still use j+1 and j-1 though, which means you're still using the second derivative along the boundary. So if you keep those and ditch the T(i-1) term, you're left with T(i,j), T(i,j+1), T(i,j-1), T(i+1,j). This means you have a second derivative along the boundary and a first derivative normal to the boundary, which is what the heat equation and your boundary conditions say.

So, when putting the T equations in the right spot, look at where you are on the boundary and figure out which derivatives you need to worry about there and then that will tell you which indexed terms need to be in that particular equation. You can also think about how a point on the top of the rectangle cannot access T(i,j+1), so that term shouldn't be there. That's another good way to check if you have the right equation in the right spot.

The full equations you will need are mostly the same as the ones from your professor's code. So use this to find which ones go where and understand why.

Here is what you had with a few small things fixed. (If you poste code again, put it between the code tags to save space and preserve indentations.)

Code:
%ftcs_block_cavity.m
%
%forward in time centered in space
%based on heat transfer energy balance at all nodes
%
%two dimensional finite difference solution to a rectangular block with
%a rectangular crack or cavity

clear all
clc

%define material properties
k = 0.401; %conduction coefficient units = W/cm*K
Cp = 385; %specific heat units = J/kg*K
rho = 0.008933; %density units = kg/cm^3

%define geometry
Lx = 50; Ly = 20;

dx = 1; %differential x length of cells
Nx = Lx/dx; %number of x nodes
dy = 1; %differential y length of cells
Ny = Ly/dy; %Number of Y nodes

tf = 60; %units of seconds
dt = 1; %seconds

Box_x = 0; %define lower left corner of box
Box_y = 0;

%Initialize mesh
T(:,:,1)=25*ones(Nx,Ny); % Set T initial
T(Nx,:,1)=100; % Set right side to initial value

time = [dt:dt:tf]; % make a vector of time values
for t = 1:length(time) % loop through the time values
% time(count) = count*dt;
for i=1:Nx %index in x direction
for j=1:Ny %index in y direction

if i==1 %identify left boundary
if j==1 %bottom left corner
% T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i+1,j,t)-T(i,j,t) + T(i,j+1,t)-T(i,j,t) + 2000 +T(i,j,t);
elseif j>1&j<Ny %left side
% T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i+1,j,t)-T(i,j,t) + 2000 + T(i,j-1,t)-T(i,j,t)+T(i,j,t);
elseif j==Ny %top left corner
% T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i+1,j,t)-T(i,j,t) + 2000 + T(i,j-1,t)-T(i,j,t)+T(i,j,t);
end
end

if i>1 & i<Nx; %identify middle boundary
if j==1; %bottom
%   T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i-1,j,t)-T(i,j,t) + T(i,j+1,t)-T(i,j,t) + T(i+1,j,t)-T(i,j,t) + T(i,j,t);
elseif j>1&j<Ny %middle
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i+1,j,t) + T(i-1,j,t) + T(i,j+1,t) + T(i,j-1,t) + T(i,j,t);
elseif j==Ny %top
% T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i-1,j,t)-T(i,j,t) + T(i+1,j,t)-T(i,j,t) + T(i,j-1,t)-T(i,j,t) + T(i,j,t);
end
end

if i==Nx %identify right boundary
if j==1 %bottom right corner
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i-1,j,t)-T(i,j,t) + T(i,j+1,t)-T(i,j,t)+T(i,j,t);
elseif j>1 & j<Ny % right side
T(i,j,t+1) = [100];
elseif j==Ny %top right corner
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i-1,j,t)-T(i,j,t) + T(i,j-1,t)-T(i,j,t)+T(i,j,t);
end
end % end right boundary if
end % end y loop
end % end x loop
end % end time loop

x(1) = Box_x;
y(1) = Box_y;
for i=2:Nx;
x(i) = x(i-1)+dx;
end
for j=2:Ny;
y(j) = y(j-1)+dy;
end
time = [0:dt:tf];
figure(1); clf;
Tfinal = T(:,:,end);
mesh(y,x,Tfinal)
title(strcat('Solution at time = ',num2str(time(end))));
xlabel('y')
ylabel('x')

true_time = time(1:end-1);
clear time;
time = true_time;

figure(2); clf;
center_y = round(Ly/2);
for e = 1:length(time)
T1(e) = T(2,center_y,e);
T4(e) = T(end,center_y,e);
end
plot(time,T1,'k',time,T4,'b')
title('Temp history of specific center line nodes')
xlabel('time')
ylabel('T ^oC')

#### mpm

My code runs now, I just need to mess with my boundary conditions a little more. I am getting a line similar to what I need but I dont think its the right one.

I will try and look at your comments more in depth and see if I can figure it out.

Thanks again, you have been a great help.

#### mpm

Oh, one more thing. I noticed whenever I change my dt and my dx and dy, my program starts running nuts. Do I need to keep one or all of these in a certain range?

For my assignment I have to change these but Im not sure if there is a boundary I need to remember.

Im guessing this hasto do with the stability criterion.

Is this correct?

#### mpm

I think I have figured out all of my code. My only concern now is the stability criteria. Im not sure what values for dt I need to use. Also Im not sure if dx effects this as well.

I believe the solution is supposed to converge but its not right now. I havea feeling this has somethign to do with the stability criteria.

Any thoughts?

Is there any method that you were told to use for finding the stability criteria? I will show you what I remember of one way, but there may be an easier way to do it.

Ok, suppose you're working with the equation $$\frac{\partial u}{\partial t} = C* \frac{\partial^2 u}{\partial x^2}$$

Assume a solution of the form $$\hat{U}(t) = e^{ikx}$$

Now if you plug it into the discrete equation

$$\frac{u^{n+1}-u^n}{\Delta t} = \frac{C}{\Delta x^2} (u^n_{i+1} - 2u^n_i + u^n_{i-1})$$

you will get

$$\frac{\hat{U}^{n+1}-\hat{U}^n}{\Delta t} e^{ikx} = \frac{C}{\Delta x^2} \hat{U}^n ( e^{ik(x+\Delta x)} - 2e^{ikx} + e^{ik(x-\Delta x)})$$
Then pull out exp(ikx) on the right and move some stuff from the left to the right to get

$$\hat{U}^{n+1} e^{ikx}= \hat{U}^n e^{ikx}+ \frac{C \Delta t}{\Delta x^2}\hat{U}^n e^{ikx} (e^{ik\Delta x} -2 + e^{-ik\Delta x})$$

Now cancel the exponential and pull out U to get

$$\hat{U}^{n+1} = \alpha \hat{U}^n$$

where alpha is all the other constant terms gathered up, something like (1+C*delt/delx^2*stuff).

Now look at that equation. It says that the U for the next time step is just U from the previous time step multiplied by a constant. So if you want your method to be stable (meaning as time gets large the U values don't go to infinity), you want alpha to be less than one. So alpha has a dependence on delt/delx^2 in it, so if you write it out and set it less than one, you will find a value that the ratio must be less than. So based on a deltaX, you can pick a deltat to keep it stable.

If you play with it for a while you should be able to find alpha. I don't want to take all the fun out of it for you. A hint is that it involes a trig function (sin^2 I think).

#### JerryWang

May i ask a question here????
If I have a 3-D wave equation, I need to solve it by finite difference method....
How can I do???
Thanks

### Physics Forums Values

We Value Quality
• Topics based on mainstream science
• Proper English grammar and spelling
We Value Civility
• Positive and compassionate attitudes
• Patience while debating
We Value Productivity
• Disciplined to remain on-topic
• Recognition of own weaknesses
• Solo and co-op problem solving