fsolve cant find a solution for some points

14 ビュー (過去 30 日間)
Robin Kirsch
Robin Kirsch 2020 年 7 月 13 日
コメント済み: Alex Sha 2020 年 7 月 14 日
Hello,
I am using fsolve to solve an equation. In some operating points, I can find a solution, whereas in others it says "No solution found". In Simulink, the system always has a stationary solution (and it matches the solution from fsolve for the operation point where it finds a solution).
For working point u = [300; 1/60] -> no solution found, whereas for u =[300;4.83/60] -> solution found and matches with the stationary solution in simulink.
Appreciate any help.
edit: I first didnt have the "dot-operator", but added it when I saw it on other questions regarding this problem that it might solve the problem. in my case it doesnt :(
EDIT: when i solve them different (single over solve) I have big "roots(..)" solutions which I solve with vpa() and then imaginary solutions which are irrelevant which I filter then aswell and I get the right solution for every working point. --> But I assume it should be possible with fsolve aswell instead of manually calculating the stationary point with solve for each equation
function F = solveme(x,u)
F(1) = (618830617597725.*x(2))./18014398509481984 - (21108551597470609.*x(1))./90071992547409920 + u(1)/5;
F(2) = ((2382265234661919.*x(1))./4722366482869645213696 - (2382265234661919.*x(2))./4722366482869645213696 + u(2).*((- x(2).^4 + (6707.*x(2).^3)./100 - 2121.*x(2).^2 + 8532.*x(2))./(200.*(x(2).^4 - (671.*x(2).^3)./10 + 2396.*x(2).^2 - 31350.*x(2) + 936000)) - (6375194751874021.*x(2))./4722366482869645213696 + 37./10000))./((4.*x(2).^3 - (20121.*x(2).^2)./100 + 4242.*x(2) - 8532)./(800.*(x(2).^4 - (671.*x(2).^3)./10 + 2396.*x(2).^2 - 31350.*x(2) + 936000)) + ((4.*x(2).^3 - (2013.*x(2).^2)./10 + 4792.*x(2) - 31350).*(- x(2).^4 + (6707.*x(2).^3)./100 - 2121.*x(2).^2 + 8532.*x(2)))./(800.*(x(2).^4 - (671.*x(2).^3)./10 + 2396.*x(2).^2 - 31350.*x(2) + 936000).^2) + 6375194751874021./18889465931478580854784);
end
u = [300; 1/60];
x0 = [90; 40];
fun = @solveme;
result = fsolve(@(x) fun(x,u),x0);
  1 件のコメント
Alex Sha
Alex Sha 2020 年 7 月 13 日
Hi, it is easy to get results like below;
x1: 290.923491047161
x2: 238.079354153592

サインインしてコメントする。

採用された回答

Matt J
Matt J 2020 年 7 月 13 日
編集済み: Matt J 2020 年 7 月 13 日
I'm not sure how Alex arrived at his results, however, his solution is reachable with even very modest improvements in the initial guess x0:
u = [300; 1/60];
x0 = [350,150];
fun = @solveme;
[result,Fopt] = fsolve(@(x) fun(x,u),x0)
result =
290.9235 238.0794
Fopt =
1.0e-08 *
-0.0000 -0.1272
  5 件のコメント
Robin Kirsch
Robin Kirsch 2020 年 7 月 14 日
which one?
Alex Sha
Alex Sha 2020 年 7 月 14 日
a package named 1stOpt

サインインしてコメントする。

その他の回答 (1 件)

John D'Errico
John D'Errico 2020 年 7 月 14 日
編集済み: John D'Errico 2020 年 7 月 14 日
It is actually pretty simple to solve this. Remember that fsolve uses a basic algorithm. Start the problem out in the wrong place though, and it will fail. I compare it to the idea of starting a blind person out at some random place on the surface of the earth, and then asking them to find some spot, using only a cane. (Disclaimer: No blind people were harmed in this thought experiment!) If you start the person out in a difficult to escape place, they will not succeed in your task.
I'll do part of this problem symbolically, because it will make things simpler to show you.
syms x y
u = [300; 1/60];
F(1) = (618830617597725.*y)./18014398509481984 - (21108551597470609.*x)./90071992547409920 + u(1)/5;
F(2) = ((2382265234661919.*x)./4722366482869645213696 - (2382265234661919.*y)./4722366482869645213696 + u(2).*((- y.^4 + (6707.*y.^3)./100 - 2121.*y.^2 + 8532.*y)./(200.*(y.^4 - (671.*y.^3)./10 + 2396.*y.^2 - 31350.*y + 936000)) - (6375194751874021.*y)./4722366482869645213696 + 37./10000))./((4.*y.^3 - (20121.*y.^2)./100 + 4242.*y - 8532)./(800.*(y.^4 - (671.*y.^3)./10 + 2396.*y.^2 - 31350.*y + 936000)) + ((4.*y.^3 - (2013.*y.^2)./10 + 4792.*y - 31350).*(- y.^4 + (6707.*y.^3)./100 - 2121.*y.^2 + 8532.*y))./(800.*(y.^4 - (671.*y.^3)./10 + 2396.*y.^2 - 31350.*y + 936000).^2) + 6375194751874021./18889465931478580854784);
First, solve the first equation for x, as a function of y. This one is just linear, so easy to solve.
x1 = solve(F(1),x)
x1 =
(1031384362662875*y)/7036183865823536 + 112589990684262400/439761491613971
F2 = subs(F(2),x,x1);
substitute that into the second equation, and plot. I'm doing this for a reason, so that you will understand how a solver gets hung up, NOT how you should necessarily have solved the problem. It would have been a great deal of help regardless.
fplot(F2,[0,400])
yline(0);
The curve seems to go to -inf as y grows large. with one root at
y0 = fzero(matlabFunction(F2),300)
y0 =
238.079354153593
But there is also a flat, shallow valley on the left, but it never crosses zero there.
The problem is though, IF you start the iterations out below around y at 150? Then the solver will get hung up in that valley to the left. This is as if you had started the blind sercher out in the vicinity of the Dead Sea, and asked him to find the lowest spot on Earth. Is it even conceivable they will manage to climb uphill to get out of the local depression you put them in, and then manage to wander to the other side of the Earth to find the Marianas trench?
Now, where did you start out fsolve?
x0 = [90; 40];
Yep. You started things out in a bad place. fsolve never had a chance to achieve success.
Could you make things simple enough to solve? Well, yes. Reduce the problem to one variable as I did. One variable problems are MUCH easier to solve than are two variable problems. Now, with a little effort, you will be able to make fzero find the root consistently.
  1 件のコメント
Robin Kirsch
Robin Kirsch 2020 年 7 月 14 日
very nice explanation, thank you

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeCalculus についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by