How can I fix my code of the heat equation?

I need to observe the behavior of the graph from the heat equation; however when I run my code I just get an orange graph that does that does not show any heat behavior, it shows the time movement in the title but not the heat movement in the graph. Any help will be greatly apreciated.This is what I have maybe someone can point me in the direction of my mistake. Thank you.

clc; close all;
N=80;
a=0;
b=1;
c=0;
d=1;
x=50;
y=50;
x=linspace(a,b,x+1);
y=linspace(c,d,y+1);
h=(b-a)/N;
[x,y]=meshgrid(x,y);
dt=0.001;
nt=5000;
k=0.01;
eps=k*dt/(h^2);
pos=[0.2,0.2];
f= exp(-((x-pos(1)).^2+(y-pos(2)).^2)./0.01);
uold=f;
unew=0*f;

for n=1:nt
    for i=2:x
        for j=2:y
            w=f(i-1,j)-2*f(i,j)+ f(i+1,j);
            v=f(i,j-1)-2*f(i,j)+f(i,j+1);
            r=w+v;
            unew(i,j)= (f(i,j)+(eps))*(r)
        
        end
    end
 unew(1,:)=0;
 unew(:,1)=0;
 unew(end,:)=0;
 unew(:,end)=0;

 uold =unew;
 pcolor(x,y,unew)
 shading interp
 colormap hot
 title(['Temperature distribution at time t= ', num2str(n*dt)])
 drawnow
end

Answer 1

[William Rose]

@Jesus,

[Edit: Correct "x(i)=pos(1) & y(j)=pos(2)" to "x(i)==pos(1) & y(j)==pos(2)".]

[Edit: Fixed the update formula.]

Am I correct to infer from your code that there is a heat source at location pos?

Am I correct that u(i,j) is the teperature at position x(i),y(j)?

What is the temperature at the source? I do not see that in the code.

Am I correct that temperature is zero along the edges?

I do not udenrstand why you use the exponential function in your solution. I realize that the exponential is the analytic solution. IT seems you are also partly using the finite difference method for the second derivative, but I don;t understand how you are combining these approaches.

I would use a simple forward finite difference method with a small time step and see what happens. I realize this is not as elegant as Crank-Nicholson or other higher order approaches, but it is simple and should work, as long as the time step is small enough compared to the spatial grid.

My loops would be something like this:

uold=zeros(nx,ny); unew=uold;
for n=1:nt
    for i=2:nx-1
        for j=2:ny-1
            if x(i)==pos(1) & y(j)==pos(2)
                unew=Tsource;
            else
                unew(i,j)=uold(i,j)+(uold(i-1,j)+uold(i+1,j)+uold(i,j-1)+uold(i,j+1)-4*uold(i,j))*k*dt/h^2;
            end
        end
    end
    uold=unew;
    %plot unew, and wait a bit before continuing
end

The loop above does not update the edge elements, so they will stay at u=0 (their initial state) throughout.

Good luck.

[Jesus]

Thank you.

[William Rose]

@Jesus, you are welcome. Please note that I had a mistake in my formula above and it is fixed now. I ran the code using the ideas above, and it worked, with k=0.01, nx=101 (which means h=0.01), and dt=0.001. The quantity k*dt/h^2 should be <0.5 for stability of the solution method I proposed. With the values I used, k*dt/h^2=0.10, so it satisfies the stability criterion. The plot looked stable after about 1-2 seconds. I ran the simulation to 5 seconds, to be sure it had reached a steady state. The final plot is below.

[Jesus]

Thank you. Your graphis more advance that what mine need to be. I only need a 2D graph yours is a 3D. I did some modifiying to the code and it is kind of working, but I am not sure if that is all it was supose to do. The graph show similar to the one I am trying to get, but I observe t change from 0 until 1 and I did not see any movement of the orange cicle(heat), so I am not sure the code is working correctly. Plus I have to use a nt big enough to heat the whole plate, and I tried some huge numbers and nothing happened, so I am not sure it is not possible or I have a mistake in the code. I was just looking at the data from my variable unew and I think I might have something backwards beacuse the values start at zero and increace slowly but after a while the values go back to zero and I do not think that suppose to happend, I attached a screen shot of what I am taking about.

clc; close all;
N=80;
a=0;
b=1;
c=0;
d=1;
xs=50;
ys=50;
x=linspace(a,b,xs+1);
y=linspace(c,d,ys+1);
h=(b-a)/N;
[x,y]=meshgrid(x,y);
dt=0.001;
nt=1000;
k=0.01;
eps=k*dt/(h^2);
pos=[0.2,0.2];
f= exp(-((x-pos(1)).^2+(y-pos(2)).^2)./0.01);
uold=f;
unew=f*0;

for n=1:nt
    for i=2:xs
        for j=2:ys
            w=f(i-1,j)-2*f(i,j)+ f(i+1,j);
            v=f(i,j-1)-2*f(i,j)+ f(i,j+1);
            r=w+v;
            unew(i,j)= f(i,j)+eps*(r);
            
        end
    end

 unew(1,:)=0;
 unew(end,:)=0;
 unew(:,1)=0;
 unew(:,end)=0;

 uold =unew;


 pcolor(x,y,unew)
 shading interp
 colormap hot
 title(['Temperature distribution at time t= ', num2str(n*dt)])
 drawnow
end

[Walter Roberson]

for n=1:nt
    for i=2:xs
        for j=2:ys
            w=f(i-1,j)-2*f(i,j)+ f(i+1,j);
            v=f(i,j-1)-2*f(i,j)+ f(i,j+1);
            r=w+v;
            unew(i,j)= f(i,j)+eps*(r);

You do not change f anywhere in the nested loops, so every iteration of n you are calculating the same thing.

[Walter Roberson]

clc; close all;
N=80;
a=0;
b=1;
c=0;
d=1;
xs=50;
ys=50;
x=linspace(a,b,xs+1);
y=linspace(c,d,ys+1);
h=(b-a)/N;
[x,y]=meshgrid(x,y);
dt=0.001;
nt=1000;
k=0.01;
eps=k*dt/(h^2);
pos=[0.2,0.2];
f= exp(-((x-pos(1)).^2+(y-pos(2)).^2)./0.01);
uold=f;
unew=f*0;

for n=1:nt
    for i=2:xs
        for j=2:ys
            w=uold(i-1,j)-2*uold(i,j)+ uold(i+1,j);
            v=uold(i,j-1)-2*uold(i,j)+ uold(i,j+1);
            r=w+v;
            unew(i,j)= uold(i,j)+eps*(r);
        end
    end

 unew(1,:)=0;
 unew(end,:)=0;
 unew(:,1)=0;
 unew(:,end)=0;

 uold =unew;


 pcolor(x,y,unew)
 shading interp
 colormap hot
 title(['Temperature distribution at time t= ', num2str(n*dt)])
 drawnow
end

[Jesus]

Thank you very much Walter Roberson, now it is working correctly that solve the issue that I was having, I really appreciate it. Thank you for you help.

[William Rose]

@Jesus, you have been helped by the person, @Walter Roberson, who knows more about Matlab, and who has helped more people, than anybody.

[Jesus]

Thank you both for your help.

[William Rose]

@Jesus, you are welcome.