Matlab - 找到两个椭圆之间的交点 - 牛顿法 - 奇怪的结果

Matlab - find intersection points between two elipses - newtons method - weird result

我正在尝试编写一个 matlab 程序,该程序应该在一个椭圆和一个弯曲的椭圆之间找到交点。

((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's 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')

图片上的圆圈是我大概的点数。如您所见,其中两点是正确的,但另外两点不正确。有谁能解释为什么值出现在椭圆不相交的地方?我试图解决这个问题几个小时而没有结果。请注意,无论我选择什么起始值,图中显示的四个点都是唯一的结果。

我没有遍历所有代码,但至少 Jacobian 的第一个系数应该是 0.32 而不是 0.8 以匹配上面一行定义的函数。 错误的雅可比行列式会导致次二次收敛,尽管在某些情况下仍能收敛。