MATLAB Answers

fsolve stopped because the problem appears regular

조회 수: 6(최근 30일)
Jiali
Jiali 2021년 7월 30일
댓글: Jiali 2021년 8월 13일
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
  댓글 수: 2
Jiali
Jiali 2021년 7월 30일
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.

댓글을 달려면 로그인하십시오.

채택된 답변

Walter Roberson
Walter Roberson 2021년 7월 30일
편집: Walter Roberson 2021년 7월 30일
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
  댓글 수: 1
Jiali
Jiali 2021년 7월 31일
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

댓글을 달려면 로그인하십시오.

추가 답변(1개)

Matt J
Matt J 2021년 7월 31일
편집: Matt J 2021년 7월 31일
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
  댓글 수: 1
Jiali
Jiali 2021년 8월 13일
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.

댓글을 달려면 로그인하십시오.

태그

제품


릴리스

R2018a

Community Treasure Hunt

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

Start Hunting!

Translated by