How can I fix my code of the heat equation?

4 visualizzazioni (ultimi 30 giorni)
Jesus
Jesus il 9 Nov 2023
Commentato: William Rose il 13 Nov 2023
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

Risposte (1)

William Rose
William Rose il 9 Nov 2023
Modificato: William Rose il 10 Nov 2023
[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.
  9 Commenti
Jesus
Jesus il 11 Nov 2023
Thank you both for your help.

Accedi per commentare.

Categorie

Scopri di più su Dates and Time in Help Center e File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by