fsolve stopped because the problem appears regular

Dear all,

I tried to solve an equation as below, but "fsolve" failed.

since ω and data are known, only two variables x(1) and x(2) are required to be solved. However, the error message is shown as "fsolve stopped because the problem appears regular". How to resovle this issue?

clear all;

omega0=2*pi*599.585e12;
data=(2+0.5i)^2;

options=optimoptions('fsolve','Display','iter');
x=fsolve(@(x)rfpnk(x,omega0,data),[2*omega0,2*omega0],options);

function F=rfpnk(x,omega0,data)
    F(1)=1-x(1)*x(2)/(omega0^2+x(1)^2)-real(data);
    F(2)=omega0*x(2)/(omega0^2+x(1)^2)+imag(data);

end

[Yuanhao Zhu]

It is likely that your optimization starting point was not chosen properly. The optimization process is guided by gradient estimation. Thus, a more reasonable starting point may solve the issue. Try a random starting point if you have no idea about what the solution would be. Or if you have a legitimate guess for your solution, you can start with a number that is close to the real solution .

[Jiali]

Initially, I doubt whether starting point is proper. I only guess that the solution should be [-1000*omega0, 1000*omega0]. But even I try to use rand(1,2)*omega0 as the starting point, the "fsolve" still failed.

Answer 1

[Walter Roberson]

format long g

omega0=2*pi*599.585e12;
data=(2+0.5i)^2;

syms x [1 2]
rf = rfpnk(x,omega0,data)

rf = 

sol = solve(rf)

sol = struct with fields:
    x1: [1×1 sym]
    x2: [1×1 sym]

[sol.x1, sol.x2]

ans = 

soln = double(ans)

soln = 1×2

5.18004253580725e+15 -2.17797242982805e+16

solv = vpasolve(rf, [-1000*omega0; 1000*omega0])

solv = struct with fields:
    x1: [3×1 sym]
    x2: [3×1 sym]

[solv.x1, solv.x2]

ans = 

double(subs(rf, sol))

ans = 1×2

0 0

double(subs(rf, solv))

ans = 3×2

-2.75 2 -2.75 2 0 0

options=optimoptions('fsolve','Display','iter');
x = fsolve(@(x)rfpnk(x,omega0,data), soln, options)

                                         Norm of      First-order   Trust-region
 Iteration  Func-count     f(x)          step         optimality    radius
     0          3     1.97215e-31                      7.26e-32               1

Equation solved at initial point.

fsolve completed because the vector of function values at the initial point
is near zero as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

x = 1×2

5.18004253580725e+15 -2.17797242982805e+16

xagain = fsolve(@(x)rfpnk(x,omega0,data), [5.18e15 -2.17e16], options)

                                         Norm of      First-order   Trust-region
 Iteration  Func-count     f(x)          step         optimality    radius
     0          3     0.000154475                      5.31e-18               1

Equation solved at initial point.

fsolve completed because the vector of function values at the initial point
is near zero as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

xagain = 1×2

5.18e+15 -2.17e+16

But initial values 5.1e+15 -2.1e+16 were not good enough for fsolve() to find something it liked. And notice that even though with the 3 digits of precision of the input solution, that fsolve's f(x) is not great compared to what it is when using the full precision found by solve(): if you were to tighten the tolerances you might need more input digits.

To clarify: I use solve() and vpasolve() here to show that there are solutions and to give us an idea of where they are so that we can explore what is needed in order to get fsolve() to work. solve() and vpasolve() are not intended to be part of the permanent solution (though if you have access to the Symbolic Toolbox, then solve() finds the exact solution easily)

It turns out that you need to be somewhat precise in order for fsolve() to be able to say it is happy, indicating that the equations are numerically sensitive.

function F=rfpnk(x,omega0,data)
    F(1)=1-x(1)*x(2)/(omega0^2+x(1)^2)-real(data);
    F(2)=omega0*x(2)/(omega0^2+x(1)^2)+imag(data);
end

[Jiali]

Dear Walter,

Thank you very much for your detailed illustrations. Now I got what you mean that the equations are too numerically sensitive to be solved perfectly by 'fsolve'. Big thanks again.

Regards,

Jiali

Answer 2

[Matt J]

If you are going to solve for x(i) that are expected to be on the order of 1e15, you need to adjust all of fsolve's tolerance parameters (StepTolerance, FunctionTolerance, OptimalityTolerance, etc...) to reflect that. The default tolerance values expect x and f(x) to be of a much lower order of magnitude.

An easier way to fix it is to change the units of x:

clear all;
omega0=2*pi*599.585e12;
data=(2+0.5i)^2;
options=optimoptions('fsolve','Display','iter');
x=fsolve(@(x)rfpnk(1e12*x,omega0,data),[2,2],options)*1e12;

                                         Norm of      First-order   Trust-region
 Iteration  Func-count     f(x)          step         optimality    radius
     0          3         11.5646                      0.000531               1
     1          6         11.5636              1       0.000531               1
     2          9         11.5609            2.5       0.000531             2.5
     3         12         11.5543           6.25       0.000531            6.25
     4         15         11.5377         15.625        0.00053            15.6
     5         18         11.4964        39.0625       0.000527            39.1
     6         21         11.3941        97.6562        0.00052            97.7
     7         24         11.1423        244.141       0.000506             244
     8         27          10.522        610.352       0.000481             610
     9         30         8.85775        1525.88       0.000474        1.53e+03
    10         33         5.30144         3814.7        0.00044        3.81e+03
    11         36        0.581522        9536.74       8.53e-05        9.54e+03
    12         39       0.0340423        5547.94       6.65e-05        2.38e+04
    13         42     6.47837e-06        1388.45       1.34e-06        2.38e+04
    14         45     2.34472e-12         5.4031       7.17e-10        2.38e+04

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

function F=rfpnk(x,omega0,data)

    F(1)=1-x(1)*x(2)/(omega0^2+x(1)^2)-real(data);
    F(2)=omega0*x(2)/(omega0^2+x(1)^2)+imag(data);
end

[Jiali]

I must confess that I like your answer more although I had accept Mr. Walter's answer as the best. Thank you for your clarification.