3

Im trying to write a matlab program that is supposed to find intersection point between one elipse and one crooked elipse.

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

results

The circles on the picture shows my approximated points. As you can see two of the points are correct but the other two is not. Does anyone have an explanation why the the values appear where the ellipses are not intersecting eachother? Ive tried to solve this problem for hours without result. Note that the four points that are shown in the picture are the only results no matter what start value i choose.

Autonomous
  • 8,478
  • 1
  • 32
  • 71
Jakob Svenningsson
  • 3,625
  • 5
  • 19
  • 26

1 Answers1

4

I didn't go through all the code, but at least the first coefficient of the Jacobian should be 0.32 instead of 0.8 to match the function defined one line above. Having the wrong Jacobian can lead to sub-quadratic convergence, although still achieving convergence in some cases.

Iban Cereijo
  • 1,520
  • 1
  • 14
  • 25
  • Why is that? 0.4*2 = 0.8? Tested to change to 0.32 and i still got the same problem. – Jakob Svenningsson Apr 29 '15 at 22:06
  • It should be 0.4^2*2 = 0.32, check the parenthesis when you define the function `f = [(0.4*x(1)).^2 ....`. It's either the Jacobian or the function itself which is wrong. – Iban Cereijo Apr 29 '15 at 22:07
  • The derivative of 0.4x^2 = 0.8x. Where do you get 2*2 from? Im using implicit derivative by the way. – Jakob Svenningsson Apr 29 '15 at 22:09
  • Sorry, edited my coment, check the parenthesis in `f = [(0.4*x(1)).^2 ...`, maybe the function is wrong and not the Jacobian. The derivative of (0.4x)^2 = 0.32x – Iban Cereijo Apr 29 '15 at 22:12
  • 1
    Thanks mate that problem was the two paranthsis's. I removed the parenthesis's around 0.4*x(1) and now it works. I understand now why it was wrong. I been working with fixing this problem for too long and got all slow in my mind. You are a life saver! – Jakob Svenningsson Apr 29 '15 at 22:20