Help with a Script using nested for and fzero function

Hi guys,

I need to implement a script that will find the zeros of a function and after that plot the value of alpha(x axis) vs beta(y axis).

the function is given by:

(cos(phi)*sin(beta-acos(cos_phi))-a)+(a-cos(phi)*sin(alpha-acos(cos_phi))*exp^(-(beta-alpha)/tang(acos(cos(phi))))

The program that i wrote till now is:

for a=0:0.2:1
    for cos_phi=0:0.2:1
        for alpha=0:0.1:pi
            beta_calc= @(beta,alpha,cos_phi,a) (cos(phi)*sin(beta-acos(cos_phi))-a)+(a-cos(phi)*sin(alpha-acos(cos_phi))*exp^(-(beta-alpha)/tang(acos(cos(phi)))))
            teste=fzero(@(beta) beta_calc(beta), [pi/2 2*pi])
            plot(alpha,beta)
        end
    end
end

But obviously is not running and im stuccoed :(

Can someone help me?

I am new on matlab so this is a kind of high level for me...

Thanks in advance!

[Walter Roberson]

What is tang() ?

Answer 1

[Walter Roberson]

Code attached.

I had to assume that tang() means tan()

The search range for fzero has to be constrained so that beta >= alpha: otherwise the sign of the term in the exp() can change, leading to the problem that the endpoints of the fzero search are not of opposite sign.

You will notice that all of the graphs are exactly the same. If you work through the formula you posted, you will find that it is of the form

(something - a) + (a - something_else)

which is equivalent to (something - a + a - something_else), and so the a's cancel out.

[Jucimar Carpe]

Hi Walter, thank you for your attention!

I didn't understand why all the graphs you have plotted are equal because i expected they looked different. The solution as you proposed seems logical but the real graph i need to plot is like the attached file. Could you take a look at it?

[Walter Roberson]

I would suggest that you made a mistake in the formula. Perhaps you could post an image of the formula?

[Jucimar Carpe]

Hi Walter!

To reduce the quantity os variables to iterate, we can substitute the "phi" parameter by "acos(cos_phi)", so the expression could be followed by:

((cos_phi*sin(beta-acos(cos_phi))-a)+(a-cos_phi*sin(alpha-acos(cos_phi))*exp(-(beta-alpha)/tan(acos(cos_phi)))))

[Jucimar Carpe]

The way you implemented the scrip it generates a "figure" for each new "a" parameter, the correct graph should alocate the "a"s at the same graph.

[Walter Roberson]

Your expression uses acos(cos(phi)) which might not be the same as phi, which is a mis-optimization on your part. You need to work with phi. You can take cos(phi) at the time you assign phi, the way I show. This is a single execution of cos() for each phi, and that is more efficient (and more correct) than the multiple acos(cos_phi) that you have.

I plotted on different figures because it is impossible to tell the graphs apart for different a values. Look at what you have:

((cos_phi*sin(beta-acos(cos_phi))-a)+(a-cos_phi*sin(alpha-acos(cos_phi))*exp(-(beta-alpha)/tan(acos(cos_phi)))))

so that is

((cos_phi*sin(beta-SOMETHING)-a)+(a-cos_phi*sin(alpha-phi)*exp(-(beta-alpha)/tan(phi))))

which is

((cos_phi*SOMETHING-a)+(a-cos_phi*SOMETHINGELSE*exp(-SOMETHINGOTHER/ANDANOTHER)))

which is

((SOMETHING-a)+(a-SOMETHINGELSE))

and that rearranges to

SOMETHING - a + a - SOMETHINGELSE

which is SOMETHING - SOMETHINGELSE + (a - a)

and so a has no effect on the calculation.

This is why you need to recheck the calculation. I recommended that you post an image of the formula so that we can look over it to see where you might have made a mistake in your (). For example perhaps (cos_phi*sin(beta-acos(cos_phi))-a) should instead be (cos_phi*sin(beta-acos(cos_phi)-a)) which moves the -a part into the sin() call.

[Jucimar Carpe]

Walter, see the attached figure, this is the original formula. We can't take out the "a" because it's multiplying by an exp function.

[Jucimar Carpe]

Now the picture...

[Walter Roberson]

As expected, you had a bracket in the wrong position.

I have attached a revised version that adheres to the given graph as closely as feasible, including the tick labels and grid, and a good approximation of the colors.

I did not, however, replicate the groups of plots. I was able to determine that the rule for any given value of a was that the plot starts at alpha = asind(a) and runs until 90 + asind(a) (so, the two places where sin(alpha)=a), but I did not manage to figure out what the rule was for the height, so I was also not able to text() appropriate labels into place.

Parts of the graph are mostly empty. There are places where the fzero() is not able to find a solution because the endpoints are the same sign.

This happens in particular for a = 0.8, for which the function starts slightly negative at beta = pi/2, rises to slightly above 0 and falls again, ending notably below zero at beta = 2*pi. Even if one were to adjust the zero-crossing algorithm to be able to find a solution, one would have to ask which of the two solutions would be appropriate.

[Jucimar Carpe]

Hi walter, thank you very much for you support!

I was trying to understand the script that you wrote and the graph. I was not able to understand why the graph i sent to you is different the graph the graph your program does. I mean, only graph all the alphas starts plotting at 90º and up, on your graph the alpha starts to plot at 90º and goes under.

Could help to understand that?

[Walter Roberson]

Your sample graph shows plots of a range of alpha_1 values from 0 to 180. I use that entire range as alpha values. However, for alpha values greater than 90, fzero is not able to find solutions for some combinations, and the other solutions are essentially on top of each other.

[Jucimar Carpe]

I made a few modifications on your code to match with what i desire. Example: i suppressed the num_a e etc and apply directly into the index of each parameter.

I notice 1 error on the part:

for cos_phi_idx = 1 : length(cos_phi_vals) cos_phi = cos_ phi_vals(cos_phi_idx);

you entered just with phi_vals on black marked above.

on the cos_phi_vals i started with a minimum value possible 0.0001

With this new code I've attached it seems closely to what i need, i just don't know test what a have to do to make graphs go to the right (up to 90º). :(

[Jucimar Carpe]

fsolve could be able to find another zeros that fzero isn't finding?

[Walter Roberson]

The correction to cos_phi_vals and the change to 0.0001 was the only real change in your version.

The num_* is there for efficiency: there is no point in re-calculating all those things every iteration of the loops, especially with nested loops.

Using rad2degree is good for readability but it does have more overhead than *180/pi.

Using .*exp instead of *exp does not matter because you are working with scalars at that point.

Using (tan(phi)) instead of tan(phi) makes the code more difficult to read: although it is not wrong, when the reader sees that you have emphasized the grouping like that, they have to stop and think about all the various circumstances under which the extra () level could make a difference, all the reasons why you might feel it is important to emphasize the grouping, only to realize that there are no reasons...

In the example plot, it is not clear to me whether, for any given a value, only asin(a) to 90+asin(a) are to be iterated over, or if instead the whole 0 to 180 range is to be iterated over but then the result is to be mapped (linearly?) into the asin(a) to 90+asin(a) range for presentation purposes.

[Jucimar Carpe]

HI walter, your are 100% correct abut the changes that I've made.

About the script, all the values of alpha (from 0 ˜ 180) need to be iterated. Ive noted that if you change the interval of: 

try this_beta = fzero(f, [max(alpha,pi/2), 2*pi])

to:

try this_beta = fzero(f, [pi -pi])

the graph starts to look a little bit with the graph i sent to you and i dont understand why this happens. For me its not logical that all plots starts from 90º (exactly same point).

If look on the graph i sent, all the plots starts at different points of alpha and beta given a different `a´ , so maybe this could be caused because of some wrong indexing.

If you notice, the first plot `a=0´ should start at alpha=180º, but in fact it is starting at 90º.

I don't know what to do anymore...

thanks in advance....

[Walter Roberson]

Okay, we are getting much closer. See attached.

I replaced the direct calls to fzero with calls to a routine intended to search for multiple zeros. The routine divides the given range up into the given number of subdivisions, and looks for all transitions in sign under the assumption that the user specified enough subdivisions that all of the transitions will be present. It finds extracts the transition intervals and does fzero() on each of the intervals to get the more exact location. The routine returns all of the detected transitions.

With the aid of this routine, I was able to find that each combination has between one and three zeros. Some of the zeros are not crossings, being points where the curve just touches zero from below and then goes negative again; I found that I was able to detect those by adding in the special points that are simple fractions of pi (e.g., 1/8 * pi). Or at least I was able to detect them well enough that I was able to find at least one zero for each combination.

I decided to plot all three zeros. In the graph that is produced, the largest of all of the zeros is associated with the ^ marker, including if that is the only zero. In cases where there are more than one zero, the smallest of the zeros is associated with the v marker. In cases where there are three zeros, the middle zero is associated with the + marker.

You will notice that a number of points are drawn below 90, including most of the lower and middle markers. In the plot you see from running the code, none of the upper markers are below 90. That is, however, only because I tweaked the code. In the situation where cos(phi) is sufficiently close to 1, the markers are below 90; in particular for cos(phi) = 1 exactly, the markers are at asind(a), which is never more than 90.

For alpha above roughly 130, there are lower beta greater than 90 for most a values. For a = 0.8 there are also middle beta greater than 90 -- so if the cut-off for drawing is to be [90 360] then you will still get some values plotted.

You might notice that the solid lines (upper beta) appear to vanish at sufficiently high alpha. Those lines are not absent: they are merely all the same positions, so the line drawn last is the one whose color is visible.

[Jucimar Carpe]

ok, i will try to run into.

Thank you very for your kindly support!!!!