Need help in finding all the solutions of non linear ellipse equations using Newtons Method using following code.

2 次查看（过去 30 天）
Haseeb Hashim
Haseeb Hashim 2022-2-22
评论： Haseeb Hashim 2022-2-26
I am trying to find all the solutions of non linear ellipse equations using Newtons Method using the following code but it only gives me 2 sets of answers no matter the initial guess value.
How can I find all the coordinates of intersection using the following code
func = @(x) [(x(1)+x(2)+2)^2+(x(1)+3)^2-5; 2*(x(1)+3)^2+(x(2)/3)^2-4];
% Jacobian of the Above system calculated
jacobi_func = @(x) [4*x(1)+2*x(2)+10,2*x(1)+2*x(1)+4;4*x(1)+12,2*x(2)/9];
%Define the stopping criteria i-e the tolerance value
tol = 10^-4;
% Define the initial guess
% Note that this effects our final value
% First Coordinate
x = [-1.62;1.38] ;
no_itr = 10 ;
x1 = x;
fnx1 = feval(func,x1);
i = 0;
fprintf(' Iteration| x | y | Error | \n')
while true
norm1 = norm(fnx1);
fprintf('%10d |%10.4f| %10.4f | %10.4d |\n',i,x1(1),x1(2),norm1)
jacob_fnx1 = feval(jacobi_func,x1);
H = jacob_fnx1\fnx1;
x1 = x1 - H;
fnx1 = feval(func,x1);
i = i + 1 ;
norm1 = norm(fnx1);
if i > no_itr && norm1 < tol, break , end
%if norm(fnv1) < error , break , end
end
It only gives me two sets of coordinates but these elllipses have four intersection points
The equations of ellipses are
(x+y+2)^2 + (x+3)^2 = 5
2*(x+3)^2 + (y/3)^2 = 4

请先登录，再回答此问题。

采纳的回答

Walter Roberson
Walter Roberson 2022-2-22
Look at the dimensionality of the problem. At any one time you have one pair of values in "x" (but which you treat as an x-y pair). You loop through using information from a function that has 2 outputs and 2 x 2 jacobian, so using Newton's Method, you get out 2 results -- an updated x/y pair. At the end of the loop you have at most one x/y pair; you cannot possibly get four x/y pairs out of this code.
Or are you trying to say that no matter what initial coordinates you use, that the result always settles on one of two locations, and your code does not seem to be able to reach the other locations?
If so.. what happens if you use the other locations as the initial values? Does it not decide the point is within tolerance? How about if you put in something very "close" to the other locations: does it move the small distance, or does it always diverge away to one of the two locations?
If passing in the known other locations does not decide that the point is a match, then you have an error in your calculations or an error in the convergence test.
  3 个评论
显示 1更早的评论
Walter Roberson
Walter Roberson 2022-2-22
Ran in:
Your code does not give only two possible answers. If you sweep the starting point over a grid of values, and give it enough iterations, then your x values settle down to one of four choices but your y values vary a lot.
The following graph shows the locations found as "solutions" when sweeping a number of different starting points.
I also show finding the solution symbolically. You can see that your solutions (from the graph) do not match the actual solutions well at all.
func = @(x) [(x(1)+x(2)+2)^2+(x(1)+3)^2-5; 2*(x(1)+3)^2+(x(2)/3)^2-4];
% Jacobian of the Above system calculated
jacobi_func = @(x) [4*x(1)+2*x(2)+10,2*x(1)+2*x(1)+4;4*x(1)+12,2*x(2)/9];
% Define the initial guess
% Note that this effects our final value
% First Coordinate
xv = linspace(-20,20);
yv = linspace(-20,20);
[X, Y] = ndgrid(xv, yv);
results = arrayfun(@(x,y) NM(x,y,func,jacobi_func), X, Y, 'uniform', 0);
resultsXY = cell2mat(results(:));
scatter(resultsXY(:,1), resultsXY(:,2));
syms x y
sol = solve(func([x,y]))
sol = struct with fields:
x: [4×1 sym] y: [4×1 sym]
[sol.x, sol.y]
ans = 
double(ans)
ans = 4×2
-4.4055 0.6663 -1.6241 1.3867 -1.6775 -2.1256 -4.0491 4.0238
function bestxy = NM(x, y, func, jacobi_func)
%Define the stopping criteria i-e the tolerance value
tol = 10^-4;
no_itr = 150 ;
x1 = [x, y];
fnx1 = feval(func,x1);
i = 0;
%fprintf(' Iteration| x | y | Error | \n')
while true
norm1 = norm(fnx1);
%fprintf('%10d |%10.4f| %10.4f | %10.4d |\n',i,x1(1),x1(2),norm1)
jacob_fnx1 = feval(jacobi_func,x1);
H = jacob_fnx1\fnx1;
x1 = x1 - H;
fnx1 = feval(func,x1);
i = i + 1 ;
norm1 = norm(fnx1);
if i > no_itr || norm1 < tol, break , end
%if norm(fnv1) < error , break , end
end
bestxy = x1;
end
Haseeb Hashim
Haseeb Hashim 2022-2-26
Thanks again for detailed explanation but as I said before the Jacobian at 1st row and second column is
2*x(1)+2*x(2)+4;
in place of
2*x(1)+2*x(1)+4
thats why my other answers were not converging but anyways I learned lot of other things form your answer.

请先登录，再进行评论。

更多回答（2 个）

Alan Stevens
Alan Stevens 2022-2-22
Check your Jacobian equations, especially the term:
2*x(1)+2*x(1)+4;
I think this should be
2*x(1)+2*x(2)+4;
Torsten
Torsten 2022-2-22
Maybe a symbolic approach is want you want:
syms x y
res1 = (x+y+2)^2 + (x+3)^2 - 5 == 0;
res2 = 2*(x+3)^2 + (y/3)^2 - 4 == 0;
s = solve([res1,res2],[x,y])

类别

在 Help Center 和 File Exchange 中查找有关 Symbolic Math Toolbox 的更多信息

标签

产品


版本

R2017a

Community Treasure Hunt

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

Start Hunting!

Translated by