Need help in finding all the solutions of non linear ellipse equations using Newtons Method using following code.
조회 수: 3 (최근 30일)
이전 댓글 표시
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
댓글 수: 0
채택된 답변
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
Walter Roberson
2022년 2월 22일
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.x, sol.y]
double(ans)
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
추가 답변 (2개)
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
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])
댓글 수: 0
참고 항목
카테고리
Help Center 및 File Exchange에서 Symbolic Math Toolbox에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!