How to select an interval in fsolve ?

I am using fsolve to solve a set of two non linear equations in two variables. But I have no idea which interval I should select. Is there any way or procedure to guess the interval ? please help. Thanks.

Answer 1

[Star Strider]

It depends on the functions. If there is a global minimum (one unique solution), where you start will likely not be important. If you doubt that there is a global minimum, and if there could be several regional minima, use the Global Optimization Toolbox patternsearch function. The other option is to use a genetic algorithm. Even if you do not have the Global Optimization Toolbox, genetic algorithms are easy to write in MATLAB.

[Chetna Sharma]

@Thank you Star. I tried solving it using both patternsearch and ga in optimization toolbox, but it gives me an error stating - '----------------------------- Optimization running. Error running optimization. Your objective function must return a scalar value.' Can you please guide how can I remove this and get my global minimum. Thanks Chetna.

[Star Strider]

My pleasure.

I will have to know more about the problem you want to solve, what your data are, and I will need to see your objective function.

[Chetna Sharma]

Hello Star, I tried to find out the interval by using fmincon in Global optimization toolbox.The funtion doesn't return to near zero value.Instead it gives me ...e+18 value. And also to run this you need to put a start value. But how do I know, which start value I should put to get the expected result.

My code is below, I am expecting z(2) and z(1) to be somewhat like to near 1.3 and 0.7 respectively.

function F = test5(z) a = 5*10-6; lambda = 632*10^-9; n1 = 1.455; n2 = 1.45; k0 = 2*pi/lambda; k1 = 2*pi*n1/lambda; k2 = 2*pi*n2/lambda; F(1) =(((1/2.*besselj(0,z(2))-besselj(2,z(2)))./(z(2)/a).*(besselj(1,z(2))))) +(((1/2.*besselk(0,z(1))-besselk(2,z(1))))./((z(1)/a).*(besselk(1,z(1))))).*(((k1^2).*(((1/2.*besselj(0,z(2))-besselj(2,z(2))))./((z(2)/a).*(besselj(1,z(2))))) + ((k2^2).*(((1/2.*besselk(0,z(1))-besselk(2,z(1))))./((z(1)/a).*(besselk(1,z(1))))))));

F(2) = (((((k0^2*n1^2-(z(2)/a)^2))*1/(a^2))).*((z(2)/a).^-2+(z(1)/a).^-2).^2); F = sqrt(dot(F,F)); end

Please help. Thanks Chetna.

[Star Strider]

Your function does not seem to be doing what you believe it is doing. I cannot troubleshoot it, because the code obviously works. The code you wrote seems to be the problem.

Plot this to help you troubleshoot it:

zv = linspace(0.05,1.5,25);

for z1 = 1:length(zv)
    for z2 = 1:length(zv)
        z = [zv(z1) zv(z2)];
        F_fcn(z1,z2) = test5(z);
    end
end

figure(1)
meshc(zv,zv,F_fcn)
grid on
set(gca, 'ZScale','log')
xlabel('\bfz_1')
ylabel('\bfz_2')
zlabel('\bfF(z_1,z_2)')
view([-200, 25])

You can comment-out the set call that creates a logarithmic scale for the ‘F’ axis. I added it because it makes the plot easier to read.

[Chetna Sharma]

Hello Star, Before I read your comment I was able to perform Globalsearch to find global minimum,( as I know my function has various local minimas) which gave me zmin and fval near to... e-7 or.. e-6 always, then I used this zmin as z(2) and z(1) in my program and ran it. I got the expected graphs. My concern was, is this correct ? Or, should I need to perform fmincon (using zmin got from gs) and that result z should be put in the equation as z(2) and z(1) to get the final result.Also, one thing more, fmincon gives me local minimum with fval ..e+18 (saying 'it is a possible local min'.). If I use this as z(2) and z(1), i do not expected result.Below is my code for one of my results exactly as expected, which contain z(2) and z(1) obtained by gs.. lambda = 632*10^-9; k0=2*pi/lambda; n1 = 1.45; a = 52.5*10^-6; n2 = 1.43; u = 6.68583; w = 6.0539; N= 490; J = zeros(0,490); K=zeros(0,490); r = zeros(0,490); f = zeros(0,490); psi = zeros(0,490); r1 = zeros(0,490); rmax = 52.5*10^-6; rmin = 0*10^-6; r1min = 52.5*10^-6; r1max = 200*10^-6; for n=1:N; r(n)= ((n-1)*(rmax-rmin))/(N-1)+rmin; r1(n)= ((n-1)*(r1max-r1min))/(N-1)+r1min; end

if  (r<=a)
   f  = (1/besselj(0,u).*(besselj(0,u.*(r/a))));

elseif  (r1>=a)
    f1  = (1/besselk(0,w).*(besselk(0,w.*(r1/a))));
end

plot((r/a),(1/besselj(0,-1.8232).*(besselj(0,-1.8232.*(r/a)))),'b',(r1/a),(1/besselk(0,1.9620).*(besselk(0,1.9620*(r1/a)))),'b'); axis = ([ 0,rmax 0,r1max]); hold on grid on

....... Thanks Chetna.

[Star Strider]

My pleasure.

If GlobalSearch found the global minimum that works for you, and especially if it gives the desired plot, go with those results! I doubt fmincon or any other optimisation routine could improve on them.

[Chetna Sharma]

Hello Star, I still have a doubt. How would I know what x0 should I choose in globalsearch, and how to set lb and ub limits ? The result of course depends on the selection of x0. So, how do I know, that I have chosen the right one ? To me it seems like manipulating your result by changing x0 to get what you want? Please help Chetna.

[Star Strider]

Unless you have a legitimate reason for constraining the parameters (for example to ensure real or positive parameter estimates), I would not constrain them. Let it find the global minimum.

[Chetna Sharma]

Hello Star, I tired doing without constraints but then neither I get the desired result nor the fval was...e-6 or -7. It was e+13 instead. Which means it needs constraint. Then question still remains the same; what criteria to follow for its selection ?

Thanks Chetna.

[Star Strider]

I cannot tell you the criteria. You need to experiment until you get an acceptable result. The globalSearch function is likely the best.

Be sure to carefully proof-read your function to be certain that you typed in what you intended to type in. If you have the Symbolic Math Toolbox, the simplify and pretty functions can be quite helpful. Just declare the variables in the syms line, paste your expression (not the anonymous function but the contents of it), then run the functions and see the result.

[Chetna Sharma]

Hello Star, It was a great help from you! My program is working well now with globalsearch, with the intervals I put. Now, I want to generalize my program for any given range of numbers. I have three variables which should be arranged in such an equation which would give me the desired result. How can I model such an equation containing three variables ?

[Star Strider]

Thank you.

If you are using a function file (rather than an anonymous function), you would have to pass the variables as extra parameters, then use an anonymous function in your code to take only your parameter vector as a formal argument.

Example:

function f = my_fcn(z,a,b,c);
  ... CODE ...
end

with ‘a’, ‘b’, and ‘c’ already in your workspace, then in your code:

call_fcn = @(z) my_fcn(z,a,b,c);

and use ‘call_fcn’ in your optimisation routine call. (Obviously, name them whatever you want. The names here are only for illustration.)

[Chetna Sharma]

Hello Star, My program is not a function file. I need to build an equation for my interval in gs. Like if my interval is 'lb' [p,q] and 'ub' is [y,s]. I have fixed p and y and q =s. q depends on three parameters a,b and c. I want to build an equation for q containg a,b and c.Now how should I proceed ??

Thanks, Chetna.

[Star Strider]

I’m sorry. I’m lost.

I have absolutely no idea.

Answer 2

[Chetna Sharma]

Hello Star, I tried to find out the interval by using fmincon in Global optimization toolbox.The funtion doesn't return to near zero value.Instead it gives me ...e+18 value. And also to run this you need to put a start value. But how do I know, which start value I should put to get the expected result.

My code is below, I am expecting z(2) and z(1) to be somewhat like to near 1.3 and 0.7 respectively. function F = test5(z) a = 5*10-6; lambda = 632*10^-9; n1 = 1.455; n2 = 1.45; k0 = 2*pi/lambda; k1 = 2*pi*n1/lambda; k2 = 2*pi*n2/lambda; F(1) =(((1/2.*besselj(0,z(2))-besselj(2,z(2)))./(z(2)/a).*(besselj(1,z(2))))) +(((1/2.*besselk(0,z(1))-besselk(2,z(1))))./((z(1)/a).*(besselk(1,z(1))))).*(((k1^2).*(((1/2.*besselj(0,z(2))-besselj(2,z(2))))./((z(2)/a).*(besselj(1,z(2))))) + ((k2^2).*(((1/2.*besselk(0,z(1))-besselk(2,z(1))))./((z(1)/a).*(besselk(1,z(1))))))));

F(2) = (((((k0^2*n1^2-(z(2)/a)^2))*1/(a^2))).*((z(2)/a).^-2+(z(1)/a).^-2).^2); F = sqrt(dot(F,F)); end

Please help. Thanks Chetna.