Why if condition does not store my variables?

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')

Answer 1

[the cyclist]

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.

[Meva]

Thanks so much for your answer @the cyclist. However, before posting as

if t==2

I posted

if nt==(2-tstart)/dt + 1

which gives t=2. It still did not work unfortunately!

Please do you have any idea?

[the cyclist]

After a closer look, I'm pretty convinced I'm right. Executing the code you posted, it does not get into this conditional:

      if  t==4   %%(nt==(4-tstart)/dt +1)  
          u4=u;
          p4=p;
          gap4=gap;
      end

but when I substitute this code instead:

      if  abs(t-4)<1.e-8   %%(nt==(4-tstart)/dt +1)  
          u4=u;
          p4=p;
          gap4=gap;
      end

it does get into it. When it does, the value of t-4 is -4.440892098500626e-16.

Note that if I simply display the value of t, even in "long" format, it will display as

4.000000000000000

but that

sprintf('%22.18f',t)

will display

ans =

    3.999999999999999556

This is an expected level of floating point imprecision.

[Meva]

Thank you very very much @cyclist. When you say

1.e-8

is that 10^(-8) ?

When I decrease time step

dt

to

dt = 10^(-9)

I must use something smaller than that , e.g.

1.e-10

Am I right? Thanks so much.

[Meva]

Also, there is another thing that bothers me a lot. The results for gap and u are very small now, comparing to the previous results which were obtained using if conditions such that

if (t==5)

So, I was wondering why the current results and the previous results are so different. If previously Matlab found the results using

    if (t==5)

, is not it the results for t=5 ? Or what were previous results representing?

[the cyclist]

Yes, 1.e-8 represents 10^(-8).

My best guess regarding your point about the gaps is that if you set the tolerance too large, then there could be many points which will trigger the condition, some of which are too far away from the values of t you want. Only the last one gets saved, not the one that is closest to the actual value of t.

[Meva]

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]

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?