Matlab - find intersection points between two elipses - newtons methods - strange result

I am trying to write a Matlab program that needs to find the intersection point between one elipse and one elipse curve.

((x − 4)^2/a^2)+((y-2)^2/2)=1 - equation for elipse

0.4x^2+y^2-xy = 10 - equation for crooked elipse
r = sqrt(10/((0.4*cos(x)^2)+sin(x)^2-cos(x)*sin(x)) - crooked elipse in polar form

fi =-pi:pi/100:pi;
a=4;b=6; % Half axis's
xc=4;yc=2;
xx=xc+a*cos(fi); yy=yc+b*sin(fi);
plot(xx,yy,xc,yc,'x')

grid;
hold on
% Non polar form x^2+y^2-xy = 10
y = sqrt(10./((0.4.*cos(fi).^2)+(sin(fi).^2)-(cos(fi).*sin(fi))));
polar(fi,y)
xstart = [-4 -4]'; % This is just an example, ive tried 100 of start  values
iter=0;
x = xstart;  
dx = [1 1]';
fel=1e-6;
while norm(dx)>fel & iter<30
    f = [(0.4*x(1)).^2+x(2).^2-x(1).*x(2)-10
        (((x(1)-4).^2)/a.^2) + (((x(2)-2).^2)/b.^2)-1];
    j = [0.8*x(1)-x(2) 2*x(2)-x(1) % Jacobian
        2*((x(1)-4)/a.^2) 2*((x(2)-2)/b.^2)];
        dx = -j\f;
        x=x+dx;
        iter=iter+1;
        disp([x' dx']);

end

iter
x
plot(x(1),x(2),'o')

      

results

The circles in the picture show my approximate points. As you can see that two points are correct and the other two are not. Does anyone have an explanation why the values ​​appear where the ellipses don't intersect with each other? Ive tried to fix this problem for hours to no avail. Note that the four points shown in the figure are the only results regardless of which starting value I choose.

+3


source to share


1 answer


I haven't gone through all the code, but at least the first Jacobian coefficient should be 0.32

instead 0.8

to match the function defined one line above. Having an incorrect Jacobian can lead to subquadratic convergence, although convergence is still achieved in some cases.



+4


source







All Articles