Why if condition does not store my variables?

4 次查看（过去 30 天）
Meva
Meva 2016-12-3
评论： Meva 2016-12-4
Hello,
Please forgive me my post is long. But I debugged and I came up with no solution. I want to store my u, p and gap arrays at time t=2, t=3, t=4, and t=5, then plot all. To do this, I put if commands associated with the times above as can be seen from the code below. At first, Matlab gave me the plots. However,now matlab says undefined variable u2 which is the first stored variable.
clear all;close all;clc;
Tsim=4;
ii=101;
dt=10^(-5);
tstart = 1;
dx=0.01;
Nt=(Tsim/dt)+1;
% Nt = 1;
ii_vect = 1:ii;
x = (ii_vect-1)*dx;
%----------------------------------------------------------------
b = 0.5; %|
kappa = 1; %|
gbar = .1; %|
M = 1; %|
g = gbar/M; %|
Gamma = 0.2; %|
F = kappa.*x.*(1-x); %|
h = kappa*b*(1-b); % h0 = F(b); %|
theta = kappa*(1-2*b); % htheta0 = F'(b); %|
%----------------------------------------------------------------
pi=4.*atan(1);
hD=0;
thetaD=0;
gap0 = 0.5;
cb=gap0;
ub=zeros(1,ii); u=zeros(1,ii); p=zeros(1,ii); q=zeros(1,ii); gap=zeros(1,ii);
gap=gap0+(-F+h+theta.*(x-b));
ub=(cb-hD.*(x-b)-thetaD/2.*(x-b).^2)./gap;
dc=-10^(-5);
%______________ TIME ITERATION____________%
for nt=1:Nt
flag=0;
mmm=1;
m=1;
c=cb;
t=tstart+(nt-1)*dt;
while (flag==0)
u=(c-hD.*(x-b)-thetaD/2.*(x-b).^2)./gap;
p(1)=0.5*(1-u(1)^2);
q(1)=0;
for i=2:ii
q(i)=q(i-1)-dx*(u(i-1)-ub(i-1))/dt;
p(i)=0.5*(1-u(i)^2)+q(i);
end
st(m)=p(ii)-0;
m=m+1;
if (m==2)
c=c+dc;
end
if (m==3)
c=(c*st(1)-(c-dc)*st(2))/(st(1)-st(2));
end
if (m==4)
mmm=mmm+1;
if (mmm==2)
m=1;
else
flag=1;
end
end
end
sumint1=zeros(1,ii); sumint2=zeros(1,ii);
sumint1=0.5*(p(1)); sumint2=0.5*(p(1))*(-b);
for k=2:ii-1
xx=(k-1)*dx;
sumint1=sumint1+(p(k));
sumint2=sumint2+(p(k))*(xx-b);
end
hDDOT=(sumint1*dx - M*g)/M;
hD=hD+dt*hDDOT;
h=h+dt*hD;
thetaDDOT=sumint2*dx/(Gamma*M);
thetaD=thetaD+dt*thetaDDOT;
theta=theta+dt*thetaD;
hL=h-theta*b;
hR=h+theta*(1-b);
gap=gap0+(-F+h+theta.*(x-b));
%____THESE IF CONDITIONS ARE NOT RECOGNISED !___%
if t==2 %%(nt==(2-tstart)/dt +1)
u2=u;
p2=p;
gap2=gap;
end
if t==3 %%(nt==(3-tstart)/dt +1)
u3=u;
p3=p;
gap3=gap;
end
if t==4 %%(nt==(4-tstart)/dt +1)
u4=u;
p4=p;
gap4=gap;
end
if t==5 %%(nt==(5-tstart)/dt +1)
u5=u;
p5=p;
gap5=gap;
end
end
ub(:)=u(:);
cb=c;
%_______PLOTTING_______%
figure(1)
subplot(2,2,1)
plot(x,u2,'k-')
hold on
plot(x,u3,'k:')
hold on
plot(x,u4,'k--')
hold on
plot(x,u5,'k-.')
hold on
legend('t=2', 't=3', 't=4','t=5')
subplot(2,2,2)
plot(x,p2,'k-')
hold on
plot(x,p3,'k:')
hold on
plot(x,p4,'k--')
hold on
legend('t=2', 't=3', 't=4','t=5')
subplot(2,2,[3,4])
plot(x,gap2,'k-')
hold on
plot(x,gap3,'k:')
hold on
plot(x,gap4,'k--')
hold on
plot(x,gap5,'k-.')
legend('t=2', 't=3', 't=4','t=5')

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

采纳的回答

the cyclist
the cyclist 2016-12-3
I have not gone through your code in detail, but I have an educated guess.
When you do your tests such as
if t==2
t is a floating-point variable which it looks to me might not be exactly equal to 2 (or whatever), because you are calculating it via a for loop that will not lead to exact integers. [It might be equal to 2.000000000001, for example.]
I expect you need to either calculate t in a way that ensure it is exactly equal to the integer you expect, or make your comparison in a way the includes a tolerance.
  7 个评论
显示 5更早的评论
Meva
Meva 2016-12-3
Very many thanks @cyclist ! So, would you please mind giving me a recommendation for time iteration then?
t
is the floating point now. So how can I set t iteration and how can I declare t so that it would not be a floating point ? Or is there any alternative way to do this:
Tsim=4;
dt=10^(-5);
tstart = 1;
Nt=(Tsim/dt)+1;
for nt=1:Nt
t=tstart+(nt-1)*dt
...
end
Thanks so much. It will solve my problem and anyone having trouble with floating points.
Meva
Meva 2016-12-4
Hello again @cyclist,
Having look at this on the web, the best solution I have come across is
Nt=400001;
instead of the Nt line in above. Just that. I was wondering that if I have any better way?

请先登录，再进行评论。

更多回答（0 个）

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by