Fzero 3 variables 1 equation optimisation

I have this equation:

(P.^3)+(0.5.*H)+Q+(1/sin(H*P))=0

I have this function saved in its own func.m file:

function A = func(x,y)
A = (P.^3)+(0.5.*H)+Q+(1/sin(H*P))
end

So it has 3 variables: P,H,Q and I cannot rearrange it in terms of H.

I need to find out for each combination of (P,Q) i.e. for P from 1 to 5 and Q from 1 to 5 the value of H for which A=0. This should give me a matrix of the 25 H values.

I am aware to use "fsolve" but having problems how to do this.

I have looked on Mathworks help which says to do this:

fun = @cos; % function
x0 = [1 2]; % initial interval
x = fzero(fun,x0)

Though I don't understand how to adapt that first line (with the @) to call my function and how to automatically find the interval.

At the moment what I do is choose a value of P and a value of Q and a range for H. Pass these 3 to the function and get back a list of A values. Look at A values to see where it passes through 0 and then set the corresponding values of H either side as the interval in the second bit of code of the fzero function. This gives me a result. Repeat for the other cominations of (P,Q). This is not the best way to do it I know and it time consuming.

How do I code so that I can put in a range of P and a range of Q, and get back a matrix of H values for which that equation equals zero? The intention is then to plot H in terms of P and Q as axes.

Answer 1

[Matt J]

Why not just do

 H=2*(A-Q-P.^3)

[James]

There was an error in my initial post. I have now corrected it. The point it that I cannot just rearrange the equation in terms of H. I know it must equal zero so need to find H value for each combination of P and Q.

[Matt J]

     for P=1:5

      for Q=1:5

       fun= @(H) (P.^3)+(0.5.*H)+Q+(1/sin(H*P));

       H(P,Q)= fzero(fun,H0);

       end

      end

[Matt J]

For such a small problem, this isn't a big deal, but you can optimize computation somewhat as follows

     for P=1:5
      for Q=1:5

       const = (P.^3)+Q;

       fun= @(H) 0.5.*H + (1/sin(H*P)) +const;

       H(P,Q)= fzero(fun,H0);
       end
      end

[James]

I now have this working using the code:

[Hval,fval]=fzero(@(H) func(p,q,H), [-2 -1])

From using func(p,q,H) I see there are two sign changes in the range of H from -5:-1. The first is between -4 and -3, and the second is between -2 and -1.

The solution in the interval -2 to -1 gives:

Hval = -1.8921
fval = -4.4409e-16

Which seems fine, close enough to zero._However..._

The other solution in the interval -4 to -3 gives:

Hval = -3.1416
fval = -5.2664e+14

So the value of the function -5E+14 is clearly not near to zero.

What's wrong with the code?

[Matt J]

I do not understand the change from: fun=@cos to:

It's an Anonymous Function. It's one of several ways MATLAB allows you to define parametrized functions

http://www.mathworks.com/help/optim/ug/passing-extra-parameters.html#brhkghv-9

For which I found the H value for which the equation is zero lies between -4 and -3 so I could then do ... However this is not giving me an output either.

You should be getting an output, or at least an error message. Here is what I get,

      >> fun = @(H) (P.^3)+(0.5.*H)+Q+(1/sin(H*P));
      H0 = [-4 -3];
      H = fzero(fun,H0)

      H =

         -3.1416

Also the original problem of having to manually find out the interval and then define this into fzero.

Since your function has multiple roots, you have to do some analysis to find out approximately where they are and which interval contains the one you're interested in. There is no automated way to handle arbitrary multi-root functions.

In your problem, though, the analysis doesn't look too bad. it is equivalent to finding the intersection of a positive sloping line

 Line(H) = 0.5*H + (P.^3+Q)

and -cosec(H*P). If P and Q are always >=1, then the line will also be greater than 1 for all H>=0. From the graph of the cosecant function, one can then see that the line will have two intersections with -cosec(H*P) in every interval of the form

[2*n-1, 2*n]*pi/P,  n=1,2,3,...

and one intersection will be in the left half of the interval, the other will be in the right. A similar analysis applies for H<0.

[James]

I should have included a graph http://postimg.org/image/qikpocg3h/ plotting it makes more sense of the H values found for function equal to zero. The first solution found is fine but the other, whilst when the function is zero, seems to be when it's shot off to infinity hence the unexpectedly large value.

It seems the way to go is to take wide range of H, spot sign change and use those points as interval to get precise value of H for function at zero. Repeat for other intervals. I had thought there might be a way to do this in code (something like take H if H+1 changed sign then set H and H+1 as interval).

[Matt J]

I should have included a graph

Why? Are we talking about a different function now? As I explained, it is very easy to understand graphically and analytically where the search intervals are going to lie for the function

(P.^3)+(0.5.*H)+Q+(1/sin(H*P))

[James]

I'm just getting muddled now. The main point of this question was help with understanding fzero and that is working now.