A convergence prcedure that doesn't work

1 次查看（过去 30 天）
vaggelis vaggelakis
Hello!!
I use the Crank-Nicolson numerical method to calculate temperatures on a metalic plate(aluminium) 10cmx30cm.
I assume that the left vertical side is always in 0 degrees Celsius and the right vertical side in 100 degrees.Both horizontal sides are insulated. So we have a heat diffusion problem. Crank-Nicolson numerical method, says in few words, we divide the plate in intervals using points, each point 's temperature is calculated by using the temperatures of its adjacent points in x,y directions.
T(i,j) = 1/2*a*dt*[E(i-1,j)-2E(i,j)+E(i+1,j)+ T(i-1,j)-2T(i,j)+T(i+1,j)]/dx^2 +
+ 1/2*a*dt*[E(i,j-1)-2E(i,j)+E(i,j+1)+ T(i,j-1)-2T(i,j)+T(i,j+1)]/dy^2 +
+E(i,j)
E is the temperature for the previous time point(we have also divided time in dt intervals) For t=0 E is known(20 degrees)
For the uknown values of T we use an arbitary value(lets say 1) and by repeating the calculation for all the points' temperatures T values must converge on the right ones.
I wrote a very simple algorithm that calculates T(i,j) for the first second(dt=1) when the plate is divided by 2cm intervals (dx=dy=2cm this means i have 6 points on y and 16 on x direction) All i care about is the convergence procedure. It works fine!!
THE PROBLEM IS, it doesn't work if i just use dx=dy=1cm which means a 11X31 matrix Here is my code that works
clear all close all clc
E=zeros(6,16); T=ones(6,16);
dx=2; dy=2; dt=1; a=0.8418; for i=1:6 for j=2:15 E(i,j)=20; end end
for i=1:6 E(i,1)=0; E(i,16)=100; T(i,1)=0; T(i,16)=100; end
for m=1:50 for j=2:15 T(1,j)=(4/3)*T(2,j)-(1/3)*T(3,j); T(6,j)=(4/3)*T(5,j)-(1/3)*T(4,j); end for i=2:5 for j=2:15 T(i,j)=(1/2)*dt*a*(E(i-1,j)-2*E(i,j)+E(i+1,j)+T(i-1,j)-2*T(i,j)+T(i+1,j))/dx^2+(1/2)*dt*a*(E(i,j-1)-2*E(i,j)+E(i,j+1)+T(i,j-1)-2*T(i,j)+T(i,j+1))/dy^2 + E(i,j);
end
end
end
for j=2:15
T(1,j)=(4/3)*T(2,j)-(1/3)*T(3,j);
T(6,j)=(4/3)*T(5,j)-(1/3)*T(4,j);
end
and here is the same code that doesn't work
clear all close all clc
E=zeros(11,31); T=ones(11,31);
dx=1; dy=1; dt=1; a=0.8418; for i=1:11 for j=2:30 E(i,j)=20; end end
for i=1:11 E(i,1)=0; E(i,31)=100; T(i,1)=0; T(i,31)=100; end
for m=1:50 for j=2:30 T(1,j)=(4/3)*T(2,j)-(1/3)*T(3,j); T(11,j)=(4/3)*T(10,j)-(1/3)*T(9,j); end for i=2:10 for j=2:30 T(i,j)=(1/2)*dt*a*(E(i-1,j)-2*E(i,j)+E(i+1,j)+T(i-1,j)-2*T(i,j)+T(i+1,j))/dx^2 + (1/2)*dt*a*(E(i,j-1)-2*E(i,j)+E(i,j+1)+T(i,j-1)-2*T(i,j)+T(i,j+1))/dy^2 + E(i,j);
end
end
end
for j=2:30
T(1,j)=(4/3)*T(2,j)-(1/3)*T(3,j);
T(11,j)=(4/3)*T(10,j)-(1/3)*T(9,j);
end

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by