find a zero of a two-variable function

my main.m is:

clear
clc
close all
%
X0 = [9.609 , 32.288]; %initial value close to zero of function
%
f = Trajectory(X0);

The script Trajectory.m is a messy think that returns a value of f. I want to find values of tht two input variables (X0) is the initial vbalues that give f=0

[Dyuman Joshi]

See fsolve (requires Optimization Toolbox)

[Torsten]

f = 0 or f = [0,0] ?

For f = 0, try "fminunc" on f.^2.

[Bruno Luong]

FSOLVE should also work for 0 , no? I believe it calls FMINCON to minimize norm(f(X))^2 in turn.

It doesn't have unique solition for sure in this case

[Sam Chak]

Hi @Douw Gerbrand

Feel free to correct me if my visualization is inaccurate. In general, a typical bivariate function will produce a surface. If you aim to solve the equation f(x, y) = 0, the real solution will yield a line segment, unless the surface is a Light Cone. For instance, solving will lead to the real solution , a straight line.

What outcome do you anticipate?

[X, Y]  = meshgrid(-1:0.05:1);
Z       = X.^3 + Y.^3;
surfc(X, Y, Z, 'FaceAlpha', 0.5)
shading interp
xlabel x, ylabel y, zlabel f(x,y), title('Two-variable function, f(x, y)')

[Douw Gerbrand]

Thank you all for your rich and useful answers to my rather ignorant question. I am a recent refugee to MATLAB from Fortran! I am still getting used to the variety of implicit functions, and their uses and context.

Yes, it is indeed a 2D contouring problem, but the "f" returned by my script "Trajectory" is actualy a positive definite cost function. "f" has only one zero in a wide domain around the intial values for which I have pretty good first estimates.

I will work on the various approaches to see what works best for me. Thanks again.

[Sam Chak]

Hi @Douw Gerbrand

Initially, I assumed that you had a predetermined path along a fixed surface and were seeking to determine a set of 2D coordinates on that static surface geometry, as demonstrated by @John D'Errico. However, based on your clarification, it seems that you are dealing with a motion system, indicated by the keywords "Trajectory" and "ODE."

Consequently, I am now attempting to visualize the three-dimensional shape of your positive definite cost function. The two parameters in the vector X are not coordinates per se, but rather factors that influence the trajectory motion. These factors can determine whether the trajectory is longer or shorter, and faster or slower in reaching a given endpoint. Such variations will directly impact the "shape" of the positive definite cost function.

Interestingly, doesn't this situation resemble the Brachistochrone problem?

[Sam Chak]

Is 'Trajectory' a script or a function? If it is a function, what is the output of the function?

function output = Trajectory()
    output  = ...;
end

[Douw Gerbrand]

Hi Sam:

Oh wow! I have not heard the term Brachistochrone since my second year university course on applied mathematics. Yes, in some senses it is that, but with additional complications. My forces include gravity, but also two aerodunamic forces, parallel to motion and perpendicular to motion. This m,akes the script "Trajecrtory" quite complicated, and non-linear. Anyway, Bruno spotted my error and I now have a working code.

T^hanks to all!

Douw

Answer 1

[Bruno Luong]

If you are not unlucky

Trajectory = @(x) sum(x.^2,2)-34^2-1; % Test function, use use rather own

X0 = [9.609 , 32.288]; %initial value close to zero of function
f1 = @(dx1) Trajectory(X0+[dx1 0]);
dx1 = fzero(f1,0);
X = X0+[dx1 0];
disp(X)

   10.6998   32.2880

Trajectory(X)

ans = 0

fimplicit(@(x,y) Trajectory([x(:) y(:)]).')
hold on
plot(X(1),X(2),'or')

[Douw Gerbrand]

Hi Bruno:

My "Trajectory" is a fairly complicated thing. It involves numerical solution of a pair of non-linear ODE's. I believe it would therefore not be possible to specify it using an "@" as you suggest. Is there a way to use your approach but not using the "@"?

[Bruno Luong]

As I said I gave a simple Trajectory as example, you don't have to use it but plug in your complicated code.on this

X0 = [9.609 , 32.288]; %initial value close to zero of function
f1 = @(dx1) Trajectory(X0+[dx1 0]);
dx1 = fzero(f1,0);
X = X0+[dx1 0];
Trajectory(X)
disp(X)

Answer 2

[John D'Errico]

A two variable function is a surface. (@Sam Chak said this, but did not take the idea to the point of completion.)

But if a surface has any solution, where it becomes zero, then almost always it will have infinitely many solutions. In fact, the tool which draws the zeros of such a function is usually called contour. That is, you want to effectively find a contour plot.

Yes, a contour plot usually is drawn with many contours at different levels, but a contour plot is what you want, all the same. You want to see only ONE contour drawn, at z == 0. Again, a contour plot can be thought of as the set of values (x,y), suze that z(x,y) is a level surface, so z(x,y) is equal to some given constant, here zero.

For an example, I'll use a symbolic function, as there are nice tools we can use. But a function handle will also work, or just an m-file. ANY function of two independent variables will work. Of course, you want it to be a continuous function, as if not, things will get very messy.

syms x y
z = x^3 - y^3 + x^2*y^2 - 3

z = 

again, we cannot simply solve for a zero of that function, as there will be infintiely many such zeros, if any exist. We might do this:

H = fcontour(z);
H.LevelList = 0;
grid on

The curve drawn in cyan is the locus of points that solve the problem z(x,y)==0, so the level set of z, at z==0. As I said, not one solution. The problem is, a tool like fsolve will produce only one solution, it will not understand there are infinitely many such solutions, just finding only one point on one of those curves, and the solution it does find will be a function of the starting value you give it.

And of course, finding the "equation" of those curves here will be a very complicated looking thing, and that is for a rather simple problem to write. We can actually find at least a set of points that will approximate those curves, using the contourc function.

So what is it you want to see? Just a plot of the solutions? Then use tools like fcontour, or contour. I could also have used fimplicit. If your goal is a set of solutions, then we could use contourc. For example...

Xv = linspace(-5,5);
Yv = linspace(-5,5);
[X,Y] = meshgrid(Xv,Yv);
% Note my careful use of the dotted operators here to vectorize the code
% and make it properly evaluated for arrays X and Y.
zfun = @(x,y) x.^3 - y.^3 + x.^2.*y.^2 - 3;
Z = zfun(X,Y);
xy = contourc(Xv,Yv,Z,[0 0])

xy = 2x212

0 -2.3787 -2.3737 -2.3603 -2.3420 -2.3238 -2.3058 -2.2881 -2.2727 -2.2706 -2.2528 -2.2354 -2.2183 -2.2017 -2.1856 -2.1717 -2.1700 -2.1547 -2.1402 -2.1267 -2.1143 -2.1032 -2.0939 -2.0865 -2.0817 -2.0800 -2.0825 -2.0906 -2.1062 -2.1328 64.0000 5.0000 4.9712 4.8990 4.7980 4.6970 4.5960 4.4949 4.4054 4.3939 4.2929 4.1919 4.0909 3.9899 3.8889 3.7978 3.7879 3.6869 3.5859 3.4848 3.3838 3.2828 3.1818 3.0808 2.9798 2.8788 2.7778 2.6768 2.5758 2.4747

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

These are the solutions, split into two curves. Unfortunately contourc puts them into one array. Looking at the first column, we see [0;64]. That tells us the first contour had z==0, and there were 64 points found on that curve. The second branch of solutions had 146 points on it that were found.

So, if you want a set of solutions, then the simple approach is to use contourc. If you want only a plot, then just use one of the contour plotting tools.

[Douw Gerbrand]

Hi John:

Thanks for your useful ideas. Yes, it is a 2D contouring idea, but as explained to Bruno, f is positive definite and has only one zero.

Douw

[Douw Gerbrand]

all I need is the values of the two arguments in "Trajectory.m" that produce a zero value for "f" which is a positive definite cost function (calculated in "Trajectory.m" and returned to my main.m)

I expect to implement some code in main.m that repeatedly calls "Trajectory.m" and finds the values for the two-vector X that results in f=0.000, to within some small tolerance.

[John D'Errico]

Yes. then you are NOT looking to find a zero of the function, but its MINIMUM value. That is simply achieved using one of many tools, perhaps fminsearch, perhaps fminunc, or fmincon, or GA, etc.

[Douw Gerbrand]

Hi John:

Thank you for that. Yes, fminsearch seems a good way to go. I will try that.

Douw

[Douw Gerbrand]

Hi John: I tried this:

clear

clc

close all

%

X0 = [9.609 , 32.288]; %initial value close to zero of function

%

X = fminsearch(Trajectory,X

%

In the command window I see:

Undefined function or variable 'X0'.

Error in Trajectory (line 11)

velocity_init = X0(1);

Error in main (line 7)

X = fminsearch(Trajectory,X0)

what did do wrong?

[Douw Gerbrand]

the line was:

X = fminsearch(Trajectory,X0)

Answer 3

[Bruno Luong]

@Douw Gerbrand If your objective function f is pd quadratic form of 2 variable X, then I believe you should minimize f rather than solve for f = 0.

Also if you use MATLAB ode solver, be careful that it can be an issue of optimizing parameters, see an example here : https://www.mathworks.com/matlabcentral/answers/2105041-using-fmincon-to-solve-objective-function-in-integral-form?s_tid=srchtitle

If you have only 2 parameters to estimate you could use fminsearch.

[Bruno Luong]

Try

X0 = [9.609 , 32.288]; %initial value close to zero of function
X = fminsearch(@(X)Trajectory(X),X0);
% or alternatively a simpler
X = fminsearch(@Trajectory,X0);

Or if it throw an error, please explain us at least the declaration signature of Trajectory and what does it suppose to do precisely with the input/output arguments.

[Douw Gerbrand]

Hi Bruno

You have spotted my error! I was not paying attention to how the variable vector X was being handled in my script "Trajectrory" I fixed that, andd not my code works.

Many thanks.