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

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

Answer 1

[Walter Roberson]

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.

[Haseeb Hashim]

Yeah you are right it gives 1 pair of answers. one of my jacobian equations was slighly wrong and also the convergence test needed modification,

Thanks man

[Walter Roberson]

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]

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.

Answer 2

[Alan Stevens]

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;

[Haseeb Hashim]

Thanks man you just saved my life. Now it gives the exact answers

Live long

Answer 3

[Torsten]

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])