Double For cycle plotting

2 次查看（过去 30 天）
Alfonso
Alfonso 2023-2-10
回答： Nehemiae 2023-3-7
Hi, i need help about plotting those double for cycle in this code. When i plot the 2 for cycle, at the interface of the multilayer soil (at Z1 and T1) it seems that matlab doesn't see the 2°layer from Z1 to Z2 and from T1 to T2. The code should represent the trend of the concentration of a solute as a function of depth and time, like the image in the bottom of this question.
those images are old but it should appear something like this. The problem is that i cant connect the plot of the 2 for cycle in the same graph. Please help me <3
function twolayerupwind_2
clc
clear all
%This is the fucking file of the exam exercise
%Explicit Euler Method is used for solving the following
%initial-boundary-value problem:
%Ut+v1*Ux=alpha1Uxx Ut+v2*Ux=alpha2Uxx
%U(z,0)=0 %Cond iniziale
%U(0,t)=g1(t) U(L,t)=g2(t) %Cond. Al bordo
%Problem parameters
%---------------------layer 1----------------------
alpha1 = 20 ; th1=1; Z1 = 100; %T = 10; -
nx = 300; nt1 = 500; v1=40; Dm1=alpha1; C0=50; T1=Z1/v1;
%---------------------Interface--------------------
xh=Z1;
%profondità dell'interfaccia.
%---------------------layer 2----------------------
alpha2 = 2 ; Z2 = 100; v2=5 ; Dm2=alpha2; th2=th1*(v1/v2);
Z=Z1+Z2; nt2 = nt1; T2=Z2/v2;
%Check of problem parameters
if any([alpha1 alpha2 Z1 Z2 T1 T2 Z-xh nx-2 nt1-2 nt2-2]<=0)
error('Check problem parameters')
end
%Stability Layer 1
st = 1; dx = Z/nx;
while 2*alpha1*T1/nt1/dx^2 + T1/nt1/dx*abs(v1) > st
nt1=nt1+1;
end
%Stability Layer 2
while 2*alpha2*T2/nt2/dx^2 + T2/nt2/dx*abs(v2) > st
nt2=nt2+1;
end
%Post-stabilità andiamo allora ad esplicitare i seguenti termini:
dt1 = T1/nt1; r1 = dt1/dx^2; s1 = dt1/dx;
dt2 = T2/nt2; r2 = dt2/dx^2; s2 = dt2/dx;
%Inizialization
nt1= (T1/dt1)+1
nt2= (T2/dt2)+1
h= nx*xh/Z + 1; h = double(uint16(abs(h)));
x=linspace(0,Z,nx+1); t1=linspace(0,T1,nt1+1); t2=linspace(0,T2,nt2+1);
xh = x(h);
u=getPhi(x,nx,Z); U=u; %Fissiamo la nostra funzione concentrazione(u) pari prorpio alla
if v1>=0
p1 = alpha1*r1 + s1*v1; q1 = alpha1*r1;
else
p1 = alpha1*r1; q1 = alpha1*r1 - s1*v1;
end
if v2>=0
p2 = alpha2*r2 + s2*v2; q2 = alpha2*r2;
else
p2 = alpha2*r2; q2 = alpha2*r2 - s2*v2;
end
%Subfunction
%[Per le condizioni iniziali]
function Pisciazz = getPhi(x,nx,Z)
Pisciazz(1:nx+1,1)=0;
%Pisciazz(1:h,1)=0 %Cond iniziale asse z, 1°Strato
%Pisciazz(h+1:nz+1,1)=0 %Cond iniziale asse z, 2°Strato
end
% [Per bordo z=0]
function VarTempo = getVarTempo(t1)
%VarTempo=2
Dt=0.5; C0=50;
if t1<=Dt
VarTempo= C0;
else
VarTempo= 0;
end
end
%Subfunction [Per bordo z=L]
function NoemannNothing = getNoemannNothing(t)
NoemannNothing= u(end-1);
end
%Finite Difference Method [Forward for the time and Central for the space]
for j=2:nt1+1
u(1)=getVarTempo(t1(j));
u(2:h-1)= p1*u(1:h-2) + (1- p1 - q1 )*u(2:h-1) + q1*u(3:h);
%buona u(h)=(th1*Dm1*(-4*u(h-1)+u(h-2))-th2*Dm2*((4*u(h+1)-u(h+2))))/(-3*(th1*Dm1+th2*Dm2));
%u(h)=(th1*Dm1*(-4*u(h-1)+u(h-2))-th2*Dm2*((4*u(h+1)-u(h+2))))/(((v1*th1-v2*th2)*2*dx)-3*(th1*Dm1+th2*Dm2));
u(h)=u(h-1);
figure(1);
plot(u,x,'r*:',U,x,'k');
xlabel('C(z,t1)'); ylabel('z');
title(['t=',num2str(t1(j)),'seconds']);
axis([0 50 0 Z1+Z2]);
pause(0.01);
A(j) = trapz(x,u)
end
for j=nt1+2:nt2+1
u(h)=u(h-1);
u(h+1:end-1)= p2*u(h:end-2) + (1- p2 - q2 )*u(h+1:end-1) + q2*u(h+2:end);
u(end)=getNoemannNothing(t);
figure(1);
plot(u,x,'r*:',U,x,'k');
xlabel('C(z,t2)'); ylabel('z');
title(['t=',num2str(t2(j)),'seconds']);
axis([0 50 0 Z1+Z2]);
pause(0.01);
A(j) = trapz(x,u)
end
legend(['t=',num2str(T1)],'t=0');
figure(2)
plot(0:dt1:T1,A)
end

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

回答（1 个）

Nehemiae
Nehemiae 2023-3-7
Hello,
In order to make the two for loops update the same figure, the “hold on” command can be used in MATLAB. The code written below, and the changes made to ensure that each for loop updates the figure each iteration and clears the unnecessary points, while ensuring that both for loops affect the same figure. Also, I have made changes in the iteration conditions of the second for loop, or else it never executes.
figure;
for j=2:nt1+1
u(1)=getVarTempo(t1(j));
u(2:h-1)= p1*u(1:h-2) + (1- p1 - q1 )*u(2:h-1) + q1*u(3:h);
u(h)=u(h-1);
plot(u,x,'r*:',U,x,'k');
xlabel('C(z,t1)'); ylabel('z');
title(['t=',num2str(t1(j)),'seconds']);
axis([0 50 0 Z1+Z2]);
pause(0.01);
A(j) = trapz(x,u);
end
hold on
for j=nt1+1:2*nt2
u(h)=u(h-1);
u(h+1:end-1)= p2*u(h:end-2) + (1- p2 - q2 )*u(h+1:end-1) + q2*u(h+2:end);
u(end)=getNoemannNothing(t1);
plot(u,x,'r*:',U,x,'k');
xlabel('C(z,t2)'); ylabel('z');
title(['t=',num2str(t2(j-nt1)),'seconds']);
axis([0 50 0 Z1+Z2]);
pause(0.01);
A(j) = trapz(x,u);
hold off
end
The resulting graph is as follows.
The documentation on the hold command (https://www.mathworks.com/help/matlab/ref/hold.html) is useful in understanding the above code.

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by