ex08c.m

% Lösung von Aufgabe 8c

% Parameter
c  = 70;
m = 1;
L = 0.3;
g = 9.81;
L1 = 0.50;
L2 = 0.50;

%% Definition der Funktion F = 2d-Funktion mit 2d-Argument
wurzp = @(x, y) sqrt((x + L).^2 + y.^2);
wurzm = @(x, y) sqrt((x - L).^2 + y.^2);
f1 = @(x,y) (1 - L1./wurzp(x, y)).*(x + L) ...
          + (1 - L2./wurzm(x, y)).*(x - L);
f2 = @(x,y) (1 - L1./wurzp(x, y)).*y ...
          + (1 - L2./wurzm(x, y)).*y + m*g/c; 
F = @(X) [f1(X(1), X(2)); f2(X(1), X(2))];

% Jacobimatrix DF = 2x2-Matrix-Funktion mit 2d-Argument
f1x = @(x,y) 2 - L1*y.^2./(wurzp(x, y).^3) ...
               - L2*y.^2./(wurzm(x, y).^3);
f1y = @(x,y) L1*(x + L).*y./(wurzp(x, y).^3) ...
           + L2*(x - L).*y./(wurzm(x, y).^3);
f2x = @(x,y) f1y(x,y);
f2y = @(x,y) 2 - L1*(x + L).^2./(wurzp(x, y).^3) ...
               - L2*(x - L).^2./(wurzm(x, y).^3);
DF = @(X) [f1x(X(1),X(2)), f1y(X(1),X(2));...
           f2x(X(1),X(2)), f2y(X(1),X(2))];

%% Plotten
xvals = -0.6:0.035:0.6;
yvals = -0.6:0.035:0.4;
plotZeros(f1, f2, xvals, yvals)
frame = getframe(gcf);
imwrite(frame.cdata, 'bild88.png');

%% Startwerte ablesen (mit Lupe Genauigkeit erhöhen!)
% auf der y-Achse
x01 = [0; -0.5];
x02 = [0;  0.13];
x03 = [0;  0.25];
% neben der Achse
x04 = [0.45; 0.10];
x05 = [-0.45; 0.10];
x06 = [0.3; 0.01];
x07 = [-0.3; 0.01];

% Lösen mit Newton-Verfahren
x1 = solveNewton(F, DF, x01, 1e-6)
x2 = solveNewton(F, DF, x02, 1e-6)
x3 = solveNewton(F, DF, x03, 1e-6)
x4 = solveNewton(F, DF, x04, 1e-6)
x5 = solveNewton(F, DF, x05, 1e-6)
x6 = solveNewton(F, DF, x06, 1e-6)  % klappt nicht!
x7 = solveNewton(F, DF, x07, 1e-6)  % klappt nicht!

%% untersuchen an der Stelle  (0.3, 0.0)
% Kurven für x = L-e, L, L+e
e = 0.01;
f1L = @(y) f1(L-e,y);
f2L = @(y) f2(L-e,y);
y = -1:0.001:1;
figure('Position',[0 0 400 600]);
subplot(3,1,1);
plot(y, f1L(y), y, f2L(y), y, 0*y, 'k')
legend('F_x', 'F_y', 'Location', 'Best')
title('x = 0.29')

f1L = @(y) f1(L,y);
f2L = @(y) f2(L,y);
y = -1:0.001:1;
subplot(3,1,2);
plot(y, f1L(y), y, f2L(y), y, 0*y, 'k')
legend('F_x', 'F_y', 'Location', 'Best')
title('x = 0.30')

f1L = @(y) f1(L+e,y);
f2L = @(y) f2(L+e,y);
y = -1:0.001:1;
subplot(3,1,3);
plot(y, f1L(y), y, f2L(y), y, 0*y, 'k')
legend('F_x', 'F_y', 'Location', 'Best')
title('x = 0.31')

frame = getframe(gcf);
imwrite(frame.cdata, 'bild89.png');