Only the first if-statement block executes,

Hi there!

I wrote an if-statement, followed by an if-else statement, followed by another if-else statement, followed by the end keyword.

The goal is to use different formulas for different values of alpha, something like the below.

However, it seems that only the first if-statement block executes -- for all possible values of alpha, even when alpha already exceeds that interval.

Where is my mistake? Thanks in advance.

    alpha = linspace(0,3*pi/2,1000)
    
    if 0 <= alpha <= pi/2
    vx = ...;
    vy = ...;
    
    elseif pi/2 < alpha <= pi
    vx = ...; 
    vy = ...;
    
    elseif pi < alpha <= 3*pi/2
    vx = ...;
    vy = ...;

    end

Answer 1

[Walter Roberson]

if 0 <= alpha <= pi/2

MATLAB interprets that as

if ((0 <= alpha) <= pi/2)

the first part, 0 <= alpha, produces a logical value, 0 or 1.

The second part compares that 0 or 1 to pi/2 and finds that 0 or 1 is always less than pi/2 so the overall test always succeeds.

MATLAB is not defined as chaining operations. In practice, chaining operations like you show is permitted in Symbolic Mathematics contexts, especially involving the piecewise() operator.

It is safest to always write the expanded version,

if 0 <= alpha && alpha <= pi/2

[Noob]

Hi Walter!

I figured this out. Thanks for your help!

[Steven Lord]

Also note that Code Analyzer will warn you about this construct in sufficiently recent release of MATLAB. If you look at the attached myfun.m:

dbtype myfun.m

1     function y = myfun(x)
2     if 0 <= x <= 2
3     	y = 1;
4     else
5     	y = 2;
6     end
7     end

you can get the list of Code Analyzer messages using the codeIssues function.

issues = codeIssues("myfun.m");

myfun.m has only one Code Analyzer message. You'll have to scroll to the right to see the whole thing.

disp(issues.Issues.Description)

Expressions like a <= b <= c are interpreted as (a <= b) <= c. Typically, to test a <= b <= c mathematically, if all arguments are numeric scalars, use (a <= b) && (b <= c), otherwise use (a <= b) & (b <= c).

[Noob]

Hi Steve! I've seen that warning message, but basically ignored it, because other errors were ocurring: It turns out I can't use the && symbol for vectors, since I have alpha = linspace ( ... ), so I used the single & to fix this issue. Then, I applied Walter's answer, which I've known about (Matlab gives the warning message) for a few days but ignored, since I had other errors to fix. Now, I'll all set. Thanks!

[Walter Roberson]

When alpha is a vector you need to be careful about

if 0 <= alpha & alpha <= pi/2

This is considered equivalent of

if all(0 <= alpha & alpha <= pi/2)

and only succeeds of all of the alpha are in range.

If you have if/elseif testing vectors, then typically you should be switching to logical indexing:

mask = 0 <= alpha & alpha <= pi/2;
vx(mask) = some computation in alpha(mask)

mask = pi/2 < alpha & alpha <= pi;
vx(mask) = some computation in alpha(mask)

[Noob]

Hi Walter,

I'm actually finding that I can only run code for one interval of alpha at a time.

So alpha = linspace(0, pi/2, 100), run code, plot data, hold on, change to alpha = linspace(pi/2 + 1e-10, pi,100), run code, plot data, hold on, etc. works just fine.

However, if I try to run code for the entire alpha = (0, 2*pi, 1001) interval, then the code doesn't run, and I get an error message.

So, I guess I had thought that the role of the else-if statements were to help me not have to manually switch the alpha interval and run my code multiple times.

But it turns out this is not the case.

What do you think the issue is?

In addition to the above attempts, I also tried using a for loop, and a while loop, but neither seems better than just manually changing the alpha interval myself, and running code and plotting code again, which, to me, is a lot faster. But I am guessing there's a way to do this without having to manually switch the intervals.

Thanks!

[Walter Roberson]

You should use logical indexing like I show -- no loop needed.

[Noob]

Hi Walter,

What is mask and vx(mask)?

Thanks,

[Noob]

Let me try this now, Walter. Thanks for your time!

[Walter Roberson]

Pre-allocating:

    alpha = linspace(0,3*pi/2,1000);
    vx = nan(size(alpha));
    vy = nan(size(alpha));
    

[Walter Roberson]

mask is a temporary variable. You can use

mask = 0 <= alpha & alpha <= pi/2;
vx(mask) = some computation involving alpha(mask)
vy(mask) = some computation involving alpha(mask)

or you can instead

vx(0 <= alpha & alpha <= pi/2) = some computation involving alpha(0 <= alpha & alpha <= pi/2)
vy(0 <= alpha & alpha <= pi/2) = some computation involving alpha(0 <= alpha & alpha <= pi/2)

repeating the logical computation each time. It is more efficient to do the logical computation only once and store the result in a tempory variable, but it is up to you whether you do that or not.

[Noob]

Hi Walter,

This worked beautifully. But I must admit: I don't totally understand why.

Here's my attempt:

mask will bring back 1s and 0s. For the 1s, the computation will execute for the corresponding values of alpha.. For all the alpha's with corresponding values of 0, the computation won't execute.

That is all there is to it.

And, Matlab is happy, because logical indexing is fast, and Matlab doesn't have much more evaluating to do.

Does that sound about right?

Thanks so much, this is a super neat method.

[Torsten]

Here are some examples for this method (called "logical indexing"):

https://blogs.mathworks.com/steve/2017/03/07/logical-indexing-2/?s_tid=blogs_rc_1

[Noob]

Thanks again Walter, and thanks Torsten!

This method is super neat.

Have a great night!

[Noob]

Hi Walter,

Do you purposely want to re-use the variable mask, to have it overwritten by the next set of computations for the next interval of alpha values?

Why not use mask1 for one set of alpha, mask2 for another set of alpha, etc.?

Thanks!

[Noob]

Hi Walter,

How come re-using mask doesn't cause the calculations to be overwritten and lost?

Instead, the calculations are all kept, and plotted nicely.

I would guess to use mask1, mask2, etc. but that seems not necessary.

Thanks!

[Steven Lord]

Can you dynamically create variables with numbered names like x1, x2, x3, etc.? Yes.

Should you do this? The general consensus is no. That Discussions post explains why this is generally discouraged and offers several alternative approaches.

In the case of your code, you don't need to keep all the masks around after you're finished using them for indexing. Use one, throw it away, pick up a new mask, throw it away, etc. In that case reusing the same variable name isn't a problem. Assigning new values to that name (without using indices) throws away the old values.

[Walter Roberson]

How come re-using mask doesn't cause the calculations to be overwritten and lost?

You start out with vx and vy all nan.

Then you step case by case, selecting locations in vx and vy according to alpha, and writing to those selected locations.

If somehow you had multiple logical masks matching a single location, then you would overwrite the results at that location. For example if you accidentally coded

mask = -pi < alpha & alpha <= 3*pi/2;

then you would be overwriting previous calculations. But as long as you do not make mistakes in your logical masks, your masks should be disjoint, and nothing should be overwritten.

[Noob]

Hi Walter / Steven,

Can I use something like vx(mask) = ... if vx is a vector quantity and not a scalar quantity?

It seems to not work for vectors.

Thanks!

[Walter Roberson]

Certainly.

alpha = linspace(0,3*pi/2,1000);
mask = 0 <= alpha & alpha < pi;
vx = nan(size(mask));
vx(mask) = alpha(mask).^2;
mask = pi <= alpha & alpha < 3*pi/2;
vx(mask) = -alpha(mask);

plot(alpha, vx)

[Noob]

Hi Walter,

If vx is some vector quantity, say, a 3x1 column vector for each alpha, then it doesn't work.

Matlab gives the error:

Unable to perform assignment because the left and right sides have a different number of elements.

[Walter Roberson]

alpha = linspace(0,3*pi/2,1000);
mask = 0 <= alpha & alpha < pi;
vx = nan(3,length(mask));
vx(:,mask) = [alpha(mask).^2;
              alpha(mask).^3 - alpha(mask).^2;
              -alpha(mask)];
plot(alpha, vx)

[Noob]

Hi Walter and Steve,

All looks ok to me now. Thanks again for your help!

[Noob]

Hi Walter,

Could the above logical indexing technique cause issues for ode45?

ode45 solved my equations tonight, but it returned all NaN values.

I suspect there's an issue with using the logical indexing, like the pre-allocating, and stuff.

Not sure ...

Matlab doesn't throw any more error messages, so I did some decent debugging of my code throughout today.

The solutions are just all NaN ...

[Noob]

Hi Walter,

It indeed looks like the logical indexing technique is making ode45 return all NaN solutions. When I use a simpler function, without needing logical indexing, the code runs fine, and I get ode solutions for all time. What do you think the issue could be?

Here's a sample of my logical indexing code, for instance:

mask = 0 <= alpha & alpha <= pi/2;

               vrelB             =  nan(3, length(mask));  
               vrelBxp           =  nan(size(mask));
               vrelByp           =  nan(size(mask));
               vrelB(:,mask)     =  vrel(:,mask) + (w/2)*omega*jp*ones(1,length(alpha(mask)));          
               vrelBxp(:,mask)   =  dot(vrelB(:,mask), ip*ones(1,length(alpha(mask))));
               vrelByp(:,mask)   =  dot(vrelB(:, mask), jp*ones(1,length(alpha(mask))));    

And I write this sort of code until alpha = 2*pi.

I made beautiful figures with this technique.

However, now I want to solve equations, using ode45, and the above technique is now causing issues. I wonder if solutions are just being overwritten or something ...

[Noob]

Hi Walter,

The above function, using logical indexing, is in a function file.

Then I use anonymous functions to pass the above function into ode45.

Matlab doesn't throw any errors, but all the solutions returned are NaN.

[Walter Roberson]

The mathematics used for all of the ode*() variable-step solvers requires that the implementation is continuous to the second derivative (the functions for internal Jacobians.) If your implementation is not continuous to the second derivative, then if you are lucky ode45() returns an error about not being able to resolve at a particular time. If you are not lucky then ode45() does not notice and silently produces wrong answers.

If you are using logical indexing in your ode45 function implementation, then chances are high that the implementation is not continuous to the second derivative.

Generally speaking, you need to use event functions to detect discontinuities (up to the second derivative) and signal termination; then you restart ode45() at the terminal state of the previous ode45() call.

In the case where the discontinuities are predictable completely by time, instead of using event functions, you can instead control the tspan parameter.

[Noob]

Hi Walter,

Would switching to if / elseif statements be a better implementation for ode45?

I tried, but the code doesn't even run; I get an unrecognized variable error message from Matlab.

From the figures I made (using logical indexing), I do notice some corner points at alpha = pi/2 and 3pi/2; this is likely where the function is not differentiable. If my suspicions are true, how could I circumvent this issue of non-differentiability at these two points?

Thanks!

[Walter Roberson]

Switching to if/elseif will not help. The implementation is either twice differentiable or it is not.

There are two fundamental strategies to deal with this situation:

1. Carefully construct the implementing code so that it is twice differentiable at each boundary condition. This might involve using bridge equations near the boundary conditions to match the second derivative
2. Use event functions to detect the discontinuities and terminate ode45() at the boundary conditions, and then resume ode45() from where you left off

When there are boundaries at multiples of then it is common for the equations to discontinuous on the first or second derivative; the cases that are continuous tend not to need conditional logic to express.

[Noob]

Hi Walter!

Ok, I guess I have more math modeling work to do now! Let's see what happens. I am guessing I need to learn how to make a sine wave with kinks / corners on it ... smoothed out, so that it is differentiable everywhere. Or, the issue could be a combination of errors. But the non-smooth behavior of my function is very clear, since I plotted many figures first and foremost, before even writing a dynamics code to give to ode45. So, I now need to approximate my non-smooth function with one that is smooth. I'll read up on this, and google around to learn of what techniques are typically used. If you have any advice in this regard, please let me know! I wonder if I can get rid of using conditional logic altogether; at the moment, it seems I am stuck with it.

Thanks so much!

[Walter Roberson]

Note that the issue of boundary conditions does not typically itself lead to NaN values being returned. NaN values being returned typically is either due to a mistake in the implementation (such as forgetting the equality test of a boundary condition), or else due to dynamics that "run away" to +/- infinity.

[Noob]

What's an equality test?

Do you think my using 0 <= alpha is problematic, and that I should instead use 0 < alpha & alpha == 0?

Also, yes, I found a glaringly obvious error in my equations!

In my equations, I was dividing a vector by another vector; something like 4i + 5j divided by 8i + 10j is nonsensical.

So, I'll work to fix up this issue now.

[Walter Roberson]

Suppose that you had

mask = 0 <= alpha & alpha < pi/2;
%some code here
mask = pi/2 < alpha & alpha < pi;
%some code here

then the code would not account for the case of alpha == pi/2 . You need to be careful at the boundaries of your conditions.