Incorrect Square Wave Plot

Question

I'm embarrassed to ask this as I believe I may be missing something obvious, but I just can't see where I am going wrong. As part of a larger program I am investigating the application of discretisation methods to approximate the convection equation on a square wave. However, I have noticed that in some cases the bounds of my square wave are being incorrectly applied. It should have the initial conditions of 1 when 10.25<=X<=10.5 and 0 everywhere else. This is an example of the issue:

L=20;                   %domain length
dx=0.01;                %space step
nx=(L/dx);              %number of steps in space
x=0:dx:(nx-1)*dx;       %steps along x
sqr=zeros(1,nx);        %pre-allocate array space

for j=1:nx
    if ((10.25<=x(j))&&(x(j)<=10.5))  
        sqr(j)=1;
    else
        sqr(j)=0;
end
d=plot(x,sqr,'r','LineWidth',2);    axis([10.1 10.6 0 1.1]);
drawnow;

In this case, the wave displays incorrectly, with x=10.5 taking a value of 0 and not 1, like so:

https://dl.dropboxusercontent.com/u/8037738/project/square_wave.png

The weird thing is that if I change the domain length to a different value it sometimes displays correctly. This is when the domain length is set to 30 and it displays correctly:

https://dl.dropboxusercontent.com/u/8037738/project/square_wave_correct.png

I really don't understand because my x array is always discretised in 0.01 intervals, so it will never 'miss' 10.5 when it loops through. I hope I have explained this problem adequately enough and if it is a stupid mistake on my part I apologise in advance.

Answer 1

You forgot the end when you finish with the if

L=20;                   %domain length
dx=0.01;                %space step
nx=(L/dx);              %number of steps in space
x=0:dx:(nx-1)*dx;       %steps along x
sqr=zeros(1,nx);        %pre-allocate array space

eps = 1.e-8;    

for j=1:nx
    if ((10.25-eps<=x(j))&&(x(j)<=10.5+eps))  
        sqr(j)=1;
    else
        sqr(j)=0;
    end
end
d=plot(x,sqr,'r','LineWidth',2);    axis([10.1 10.6 0 1.1]);
drawnow;

I also recomend you to use a small tolerance when using comparisons with float values, like the x(j)<=10.5

if you se the difference at that point, matlab tell you

x(1051)-10.5

ans =

     1.776356839400250e-15

so because of the round-off the value of x(j) was lightly larger that 10.5, not equal as you were expecting

Answer 2

Your issue is that the x values do not fall exactly on the space step. For example, if you look at the following:

j=1051;
format long;
x(j)

ans =

10.500000000000002

The answer from sebas fixes this issue by adding the tolerance to each of the bounds. You could also try rounding x(j) to the nearest 100th before performing the logic in the if statement:

roundn(x(j), -2)