6 views (last 30 days)

Show older comments

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

Walter Roberson
on 30 Jul 2021

Edited: Walter Roberson
on 30 Jul 2021

format long g

omega0=2*pi*599.585e12;

data=(2+0.5i)^2;

syms x [1 2]

rf = rfpnk(x,omega0,data)

sol = solve(rf)

[sol.x1, sol.x2]

soln = double(ans)

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

[solv.x1, solv.x2]

double(subs(rf, sol))

double(subs(rf, solv))

options=optimoptions('fsolve','Display','iter');

x = fsolve(@(x)rfpnk(x,omega0,data), soln, options)

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

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

Matt J
on 31 Jul 2021

Edited: Matt J
on 31 Jul 2021

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;

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

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!