Matlab - найти точки пересечения между двумя эллипсами - метод ньютонов - странный результат

Я пытаюсь написать программу 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')

результаты

Кружками на картинке показаны мои приблизительные точки. Как видите, два пункта верны, а два других нет. Есть ли у кого-нибудь объяснение, почему значения появляются там, где эллипсы не пересекаются друг с другом? Я пытался решить эту проблему в течение нескольких часов безрезультатно. Обратите внимание, что четыре точки, показанные на рисунке, являются единственными результатами, независимо от того, какое начальное значение я выберу.


person Jakob Svenningsson    schedule 29.04.2015    source источник


Ответы (1)


Я не просматривал весь код, но по крайней мере первый коэффициент якобиана должен быть 0.32 вместо 0.8, чтобы соответствовать функции, определенной строкой выше. Наличие неправильного якобиана может привести к субквадратичной сходимости, хотя в некоторых случаях сходимость все же достигается.

person Iban Cereijo    schedule 29.04.2015
comment
Почему это? 0,4*2 = 0,8? Протестировано, чтобы перейти на 0,32, и у меня все еще есть та же проблема. - person Jakob Svenningsson; 30.04.2015
comment
Должно быть 0,4^2*2 = 0,32, проверьте круглые скобки при определении функции f = [(0.4*x(1)).^2 ..... Либо якобиан, либо сама функция неверна. - person Iban Cereijo; 30.04.2015
comment
Производная 0,4x^2 = 0,8x. Откуда у вас 2*2? Кстати, я использую неявную производную. - person Jakob Svenningsson; 30.04.2015
comment
Извините, отредактировал мой комментарий, проверьте скобки в f = [(0.4*x(1)).^2 ..., возможно, функция неверна, а не якобиан. Производная от (0,4x)^2 = 0,32x - person Iban Cereijo; 30.04.2015
comment
Спасибо, приятель, что проблема была в двух парантсисах. Я удалил скобки около 0,4 * x (1), и теперь это работает. Теперь я понимаю, почему это было неправильно. Я слишком долго работал над исправлением этой проблемы, и у меня все замедлилось. Вы спасатель жизни! - person Jakob Svenningsson; 30.04.2015