Fzero with multivariable + fixed parameters (can not find how to make it)

Defining function:

# Check the angle parabola and straight line crossing
function [th]=finda(k,t,th0);
  #th=0;
  th=sec(t+th0)*tan(t+th0)-k/sin(t);
endfunction

Then command line variables definition and using fzero():

octave:27> th0=atan(sqrt((sqrt(101)-1)/2))
th0 = 1.1314
octave:28> k=6
k = 6
octave:29> fzero(@finda(k,t,th0), [1])
error: fun(0.99): subscripts must be either integers 1 to (2^63)-1 or logicals
error: called from
    fzero at line 206 column 10

It should find value of ‘t’ and k and th0 should be considered as fixed variables for the fzero().
Something looks wrong about this.

Below would work:

function [th]=finda(t);
  th0=atan(sqrt((sqrt(101)-1)/2));
  k=6;
  th=sec(t+th0)*tan(t+th0)-k/sin(t);
endfunction
octave:36> fzero(@finda, [1])
ans = 0.7781
octave:37> fzero(@finda, [2])
ans = 3.9196

But th0 and k should not be defined INSIDE the function, but as argument.

Then if t is found, let assume t=0.7781... and r=k/sin(t)=8.5483....
It will be converted to cartesian:

octave:45> x=r*cos(t)
x = 6.088784323982802
octave:46> y=r*sin(t)
y = 5.999999999999999

however, the y has lost something and k=6 has become to 5.999999...
So the point is not exactly on the line.
it might have not a big impact though if the original value from k=6 is used for the y, but…

Eventually, if code is ready it should find all the crossing points of the parabola and the bounding box (that is at (-2,2) to (8,6)):

The coordinate of the parabola (either x or y) should be exactly integer for the coordinate that is on the bounding box line. But it seems not easily handled.
For example x=8, y=some non integer at about ~0.6 calculated by the finda.

Looking into the help of fzero it clearly says:
Find a zero of a univariate function.

If you want to use a multivariate equivalent you should use fsolve.

Only the ‘t’ is to be solved (single variable), and the ‘k’ and ‘th0’ should be able to enter
there as ‘parameters’ not to be solved. The function should accept three variables,
where only one of then is to be solved.
I might be using wrong command. I’m not sure what the manual say about this situation.
How can the ‘k’ and ‘th0’ be given as ‘fixed’ argument for function, and only ‘t’ for solving?
How can the optimization (or solving) variables be separated from other parameters to be given to the function?
Below sample case would work also:

global th0=atan(sqrt((sqrt(101)-1)/2));
global k=6;
function [th]=finda(t);
  global th0;
  global k;
  th=sec(t+th0)*tan(t+th0)-k/sin(t);
endfunction
octave:4> fzero(@finda, [1])
ans = 0.7781

But I don’t want to use global to give the ‘th0’ and ‘k’ into the function.
Also using fzero would be more convenient for the speed:

octave:6> tic; for i=1:1001, fzero(@finda, [1]); end; toc
Elapsed time is 2.82726 seconds.
octave:7> tic; for i=1:1001, fsolve(@finda, [1]); end; toc
Elapsed time is 4.62057 seconds.

But entering the extra parameters looks difficult.
There must be some other way than the global, but not seen here.
Somehow I should be able to prevent the solving/optimization of some parameters/arguments, while the others are to be solved/optimized.

octave:8> fsolve(@finda(k,t,th0), [1])
error: fcn(1+1.49012e-08): subscripts must be either integers 1 to (2^63)-1 or logicals
error: called from
    __fdjac__ at line 52 column 17
    fsolve at line 291 column 12
octave:9>

fsolve also can not work it out what parameters are to be optimized (it should be ‘t’ only).
Of course I can just use the global variables if that is only way, but I would not like to declare
these variables as global, because they are just for this function only.
Exact code is below (and with 2 functions finda & findb for y & x correspondingly):
plot_parabola.m (4.4 KB)
It would work with this ‘global’, but… is that the right way to make it?
Also if I want to get all the solutions [1]…[6] at once, is it possible?

octave:49> fsolve(@findb, [1])
ans = 0.7408
octave:50> fsolve(@findb, [2])
ans = 3.2320
octave:51> fsolve(@findb, [3])
ans = 3.2320
octave:52> fsolve(@findb, [4])
ans = 3.8824
octave:53> fsolve(@findb, [5])
ans = 6.3736
octave:54> fsolve(@findb, [6])
ans = 6.3736
octave:55> fsolve(@findb, [1:6])
error: operator *: nonconformant arguments (op1 is 1x6, op2 is 1x6)
error: called from
    findb at line 4 column 7
    fsolve at line 242 column 8

All the solutions from t=0 to t=2*pi should be generated at once.

Also the parameter yp(1) is larger than ‘k=6’:

octave:63> xp(1)
ans = 6.0888
octave:64> yp(1)
ans = 6.0000
octave:65> yp(1)-6
ans = 1.6875e-14

So it makes some issue with the plot edge.

# To avoid bounding box edge violation
yp(yp>6)=6;

And similarly at the right edge it is not exactly at the edge:

octave:7> xp(end)
ans = 8.0000
octave:8> xp(end)-8
ans = -1.1902e-13

Maybe I just modify this manually: x(end)=8
If I run below code:
for i=0:5, fzero(@findb,[i]), end
it also adds the ‘singularities’ to the solution… I’m guessing.
Using:

thmin=fzero(@finda,[1],optimset("TolX",eps));
thmax=fzero(@findb,[3],optimset("TolX",eps));

But still error exists:

octave:30> xp(end)-8
ans = -1.1902e-13
octave:31> yp(1)-6
ans = 1.6875e-14

Remaining may need be fixed manually: xp(end)=8; yp(1)=6;
It could not get any better with this option. Maybe impossible.

I couldn’t get the ‘gradient’ out:
[thmin,fval,info,g,h]=fzero(@finda,[1],optimset("TolX",2e-16*eps,"MaxIter",500,"GradObj","o

Also I think the question for how to use the fzero is also:
How can I determine direction of optimization?
For example if I start from angle t=-1.1314 and want to find
first point below t=-1.1314 or if I want find first point above t=-1.1314
Can I determine start point and direction?
Or if I want to find all the solution between -pi to +pi, then how…
Below seems to work…

k=6;
thmin=fminbnd(@finda,-pi/2-th0,-th0);
k=8;
thmax=fminbnd(@findb,-th0,pi/2-th0);

The fixed parameters (‘k’, ‘th0’ still must be entered through the ‘global’, because the system can not allow entering them through the function call (ans as a ‘fixed’ parameters).

If I want to find all the intersection point of the box and the parabola… looks difficult to be done.

Why don’t you use xlim and ylim?

Cause I want the curve end exactly at the edge (also not before edge).
But more I think it as a mathematical challenge and to see if I can get these
tools running for what I try.
For example if there is 2 solutions in the given interval (-th0…pi/2-th0):

octave:95> t1=fminbnd(@findy,-th0,pi/2-th0)
t1 = -0.050127
octave:96> t1=fminbnd(@findy,-th0,pi/4-th0)
t1 = -0.8859

I should somehow split the interval to get the both results out.
But splitting the interval is very risky, cause it can not be pre-estimated easily how to split it.
Comparing to command "roots": That can give all the roots (or zeroes).
while fzero can not give all the results (and also can not accept the solution interval as input parameter.).
Considering this parabola and a straight line, they can have 0, 1 or 2 intersection points.
But the fminbnd, can find only 0 or 1 intersection points (cases with 2 intersection points seem to be rejected).
Maybe completely different solution is required then…
Below looks fine (2 roots/results are produced):

octave:8> roots([1 0 -1])
ans =
  -1
   1

But as the parabola is in the polar format and/or tilted, that is not easily made with roots-command.

The usual way to do this is simply to wrap the base function in an anonymous function which sets the constant parameters correctly and leaves only one free variable. As pointed out, fzero only works for univariate functions. Below is an example Octave session showing that approach.

octave:4> f1 = @(k,th0,t) sec(t+th0)*tan(t+th0)-k/sin(t);
octave:5> th0=atan(sqrt((sqrt(101)-1)/2))
th0 = 1.1314
octave:6> k = 6
k = 6
octave:7> f2 = @(t) f1(k, th0, t)
f2 =

@(t) f1 (k, th0, t)

octave:8> fsolve (f2, 1)
ans = 0.7781
octave:9> fsolve (f2, 2)
ans = 0.7781
octave:10> fsolve (f2, 3)
ans = 0.7781
octave:11> fsolve (f2, 4)
ans = 3.9196

It looks hard to understand/see how it works…
For me it looks like fsolve would not take the ‘th0’ and ‘k’ as fixed argument this case.
It might work, but that is hard to see and if k or th0 is fixed number this would not work.
like in case: f1(-2,1.5,t).

fsolve requires a function with only one variable. Hence, the declaration of f2 is f2 = @(t) ..., and there is only one variable present. fsolve never “sees” k and th0 because they are not in the f2 function declaration. This is the standard solution recommended by Matlab for this use case, but you already found a solution with global variables which also works. Choose whichever one makes the most sense for you.

This is maybe good, but I can’t see why and it looks very hard to read this kind of code.
Also I think I’ve gone ‘wrong’ in this type of development, as the fminbnd, fsolve, fzero all
seem to be quite slow to operate.

octave:9> tic; for i=1:1001, roots([1 0 -1]);  end; toc;
Elapsed time is 0.18701 seconds.

This looks more than 10x faster than fzero, etc…
Maybe just convert the parabola to cartesian functions and try again.
Then it looks hard to be calculated as y=f(x) for horizontal line intersections
and with x=f(y) for vertical line sections.
Roots[] might not be needed either as it is just 2nd degree polynomial.

In case of rotation in cartesian, it need ‘parametric form’ and then apply the rotation matrix
to get it done:
https://www.algebra.com/algebra/homework/Quadratic-relations-and-conic-sections/Quadratic-relations-and-conic-sections.faq.question.715758.html
The roots then are possibly easier to be found with this form.

In case if the fzero operated with the fixed parameters directly and find all the ‘roots’ at certain interval that could be preferred, like:
fzero(@findx,xmin,xmax)
and result is all the roots (or [] if there is no root at all), not just one.
I guess it doesn’t want to work, because solution is to use parametric plot instead:
plot_parabola.m (6.1 KB)
The density of he plot is maybe difficult to adjust to be correct… It might work though…
I couldn’t find any test case when it fails.
Numerical errors could be failing it at the corner regions.

I guess it’s a matter of perspective. I think using two anonymous functions is easier than writing separate m-files and using global variables.

I don’t think so. The original equation is non-linear

sec(t+th0)*tan(t+th0)-k/sin(t)

and so requires a solver adapted for that. If you change the problem then you can change your solving technique. If you want to find the zero-crossings of a polynomial then roots which just needs to do Gaussian Elimination to get an answer is going to be faster.

I tried to get it working better, but even working, the line looks suspicious:
plot_parabola.m (5.8 KB)
ones(length(t0),1)*real(bb(1));
There might be better/shorter way to make it…
The code optimization could be harder to be made.
Also below looks similar like problem:

  tv=linspace(-a^2/h,a^2/h,91)(1:90);
  th=linspace(b^2/h,-b^2/h,91)(1:90);

But to drop out the last point seem not too easy.

The plotting might work now:
plot_ellipse.m (12.9 KB)
plot_parabola.m (5.8 KB)
but coding doesn’t look too clean…