Plotting Only 4 Time steps

조회 수: 9 (최근 30일)
Fares
Fares 2022년 11월 13일
댓글: Fares 2022년 11월 14일
Hello,
I really failed to make my plot draw only the results of 4 time spots out of the whole time steps, I tried to use hold on but it drew every single iteration. I need just to show 4 screen shots of the evoloving temperature disutribution (begining, middle, middle 2, final), can anyone help me pls.
Thank you
  댓글 수: 1
the cyclist
the cyclist 2022년 11월 13일
Here is the code, for those who do not care to download it ...
%Seconf Deravitive solution
set(0,'DefaultAxesFontName','Times New Roman')
set(0,'DefaultAxesFontSize',17)
clear *
L=0.025; % copper wire length
N = 100; % you decide, # of discretization points
% the more you increase the N number, the more it is close to the BC
NN = N + 2; % 2 represent the ghost cells.
dx = L/N;
x = linspace (-dx/2,L+dx/2,NN) % define temp. solution dx/2 3dx/2 etc..
T = 293 + sin(2*pi*x/L);
in = 2:1:(NN-1); % Internal cell
rhs = in+1 % Right hand side
lhs = in-1 %Left hand side
d2Tdx2 = (T(rhs)-2*T(in)+T(lhs))/(dx^2); %1st derivatve of the heat equation shown in the slide page 37
Ta = -sin(2*pi*x/L)*(2*pi/L)^2; %analytical
T = 293 +sin(2*pi*x/L);
%Ta = cos(2*pi*x/L)*(2*pi/L); % analytical derivative (maybe)
figure(2),clf
plot(x(in),d2Tdx2,'o-' ), hold on% we have put n to make the vectors match by excluding the ghost cells.
plot(x(in),Ta(in),'r-')
legend('Numerical d2Tdx2','Analytical dTdx')
%%%%%%%%%%%%%%%Start the HW%%%%%%%%%%
%solve dTdt = alpha*d2Tdx2
%using Euler method and central differnce
% for space derivative
alpha = 10.1; %m2/s
T = 300*ones(1,102)% Initalize
Tmin = 310; Tmax=300;
T(1) = 2*Tmin - T(2) %0.5*(T(1)+T(2)) = Tmin
T(NN) = 2*Tmax - T(NN-1);
dt = 0.000000001; %timestep in second
simutime = 1; %seconds
simusteps = round(simutime/dt);
for(t=1:simusteps)
d2Tdx2 = (T(rhs)-2*T(in)+T(lhs))/(dx^2);
dT = alpha*dt*d2Tdx2; %get dT in each cell
T(in) = T(in) + dT;
%set BC
%If you want to isolate the left side
% you can add T(1) = T(2)
% and comment the line below.
T(1) = 2*Tmin - T(2) %0.5*(T(1)+T(2)) = Tmin
T(NN) = 2*Tmax - T(NN-1);
if(mod(t,10))
figure(3);
plot(x(in),T(in),'r-','LineWidth',1), hold on
h=xlabel('Length (m)');
h=ylabel('Temperature (K)');
title('Part A');
end
end

댓글을 달려면 로그인하십시오.

답변 (1개)

the cyclist
the cyclist 2022년 11월 13일
I expect that this if statement is not doing what you expect ...
if(mod(t,10))
The lines following that will be executed whenever mod(t,10) is non-zero, which is most iterations.
  댓글 수: 3
the cyclist
the cyclist 2022년 11월 14일
I made a few changes to your code. The most important is the change to that mod() function, so now it only plots every 5,000 iterations.
Note that your code runs for a billion iterations.
set(0,'DefaultAxesFontName','Times New Roman')
set(0,'DefaultAxesFontSize',17)
clear *
L=0.025; % copper wire length
N = 100; % you decide, # of discretization points
% the more you increase the N number, the more it is close to the BC
NN = N + 2; % 2 represent the ghost cells.
dx = L/N;
x = linspace (-dx/2,L+dx/2,NN); % define temp. solution dx/2 3dx/2 etc..
T = 293 + sin(2*pi*x/L);
in = 2:1:(NN-1); % Internal cell
rhs = in+1; % Right hand side
lhs = in-1; %Left hand side
d2Tdx2 = (T(rhs)-2*T(in)+T(lhs))/(dx^2); %1st derivatve of the heat equation shown in the slide page 37
Ta = -sin(2*pi*x/L)*(2*pi/L)^2; %analytical
T = 293 +sin(2*pi*x/L);
%Ta = cos(2*pi*x/L)*(2*pi/L); % analytical derivative (maybe)
figure(2),clf
plot(x(in),d2Tdx2,'o-' ), hold on% we have put n to make the vectors match by excluding the ghost cells.
plot(x(in),Ta(in),'r-')
legend('Numerical d2Tdx2','Analytical dTdx')
%%%%%%%%%%%%%%%Start the HW%%%%%%%%%%
%solve dTdt = alpha*d2Tdx2
%using Euler method and central differnce
% for space derivative
alpha = 10.1; %m2/s
T = 300*ones(1,102);% Initalize
Tmin = 310; Tmax=300;
T(1) = 2*Tmin - T(2); %0.5*(T(1)+T(2)) = Tmin
T(NN) = 2*Tmax - T(NN-1);
dt = 1.e-9; %timestep in second
simutime = 1; %seconds
simusteps = round(simutime/dt);
for t=1:simusteps
d2Tdx2 = (T(rhs)-2*T(in)+T(lhs))/(dx^2);
dT = alpha*dt*d2Tdx2; %get dT in each cell
T(in) = T(in) + dT;
%set BC
%If you want to isolate the left side
% you can add T(1) = T(2)
% and comment the line below.
T(1) = 2*Tmin - T(2); %0.5*(T(1)+T(2)) = Tmin
T(NN) = 2*Tmax - T(NN-1);
figure(3)
hold on
if mod(t,5e3)==0
plot(x(in),T(in),'-r','LineWidth',1)
xlabel('Length (m)');
ylabel('Temperature (K)');
title('Part A');
end
end
Fares
Fares 2022년 11월 14일
That helps a lot.
Thank you very much for your help, it is most appreciated.

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 MuPAD에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by