% EX 3.9
% Skärningspunkter parabel-lemniskata
clear

% Lemniskatan: L(x,y) = (x^2 + y^2/2)^2 - x^2 + y^2/2 = 0
% Parabeln: P(x,y) = y - x^2 + 1/2

% Parabeln är enkel att plotta
xp = -1:0.1:1;
yp = xp.^2 - 1/2;
% Gå över till radiella koordinater för plottning av lemniskatan
% x = r*cos(fi)
% y = r*sin(fi)
% (r^2*cos^2(fi) + r^2*sin^2(fi)/2)^2 - r^2*cos^2(fi) + r^2*sin^2(fi)/2 = 0, r >= 0
% För varje fi = 0:(pi/10):(2*pi), lös ekvationen mha fzero
xl = [];
yl = [];
for fi = 0:(pi/100):(2*pi)
	f = @(r) (r^2*cos(fi)^2 + r^2*sin(fi)^2/2)^2 - r^2*cos(fi)^2 + r^2*sin(fi)^2/2;
	if (abs(fi - pi/2) < 0.2*pi) || (abs(fi - 3*pi/2) < 0.2*pi)
		r = 0; % r = 0
	else
		r = fzero(f, [eps 1]);  % r > 0, r <= 1
	end
	xl = [xl, r*cos(fi)];
	yl = [yl, r*sin(fi)];
end

figure(1)
clf
hold on
plot(xp, yp, 'r', 'LineWidth', 2);
plot(xl, yl, 'g', 'LineWidth', 2);
xlabel('X', 'FontSize', 14)
ylabel('Y', 'FontSize', 14)
axis([-1 1 -1 1])


% Skärning: Både L(x,y) = 0 och P(x,y) = 0
% Olinjärt ekvationssystem med två variabler u = [x y] och två ekvationer.
% Kan lösas med Newton-Raphsons metod. Startgissningar från plot:
ustart = [-0.9 -0.3  0.3  0.9
	         0.3 -0.4 -0.4  0.3];
plot(ustart(1, :), ustart(2, :), 'ko', 'LineWidth', 2)
						
% Ekvationssystem, jacobian
F = @(u) [(u(1)^2 + u(2)^2/2)^2 - u(1)^2 + u(2)^2/2; ...
	        u(2) - u(1)^2 + 1/2];
J_F = @(u) [2*(u(1)^2 + u(2)^2/2)*2*u(1) - 2*u(1), ...
	          2*(u(1)^2 + u(2)^2/2)*u(2) + u(2); ...
	          -2*u(1), 1];
% Gör en sökning för var och en av nollställena
for n = 1:4
	pause
	
	u = ustart(:, n);
	h = [1; 1];
	while norm(h) > 5e-9
		h = J_F(u)\F(u);
		u = u - h;
	end
	
	plot(u(1), u(2), 'k*', 'LineWidth', 2)

end
	