Illustration des Newton-Verfahrens fuer Systeme

clear all;
close all;

f1 = @(x,y) x.^2 + y.^2 - 6;
f2 = @(x,y) x.^3 - y.^2;

ezsurf('0',[-3 3 -3 3])
hold on;
ezsurf(f1,[-3 3 -3 3])
figure
ezsurf('0',[-3 3 -3 3])
hold on;
ezsurf(f2,[-3 3 -3 3])
figure
[X,Y] = meshgrid(-3:0.1:3,-3:0.1:3);
contour(X,Y,f1(X,Y),[0 0],'k')
hold on
contour(X,Y,f2(X,Y),[0 0],'k')
df1x = @(x) 2*x;
df1y = @(y) 2*y;
df2x = @(x) 3*x.^2;
df2y = @(y) -2.*y;
xi = 1;
yi = 1;
plot(xi,yi,'k+')

maxit = 100;
eps = 1.d0;
i = 0;
while eps > 1.d-8 && i <=maxit
    DF(1,1) = df1x(xi);
    DF(1,2) = df1y(yi);
    DF(2,1) = df2x(xi);
    DF(2,2) = df2y(yi);
    f = [-f1(xi,yi);-f2(xi,yi)];
    delta = DF\f;
    xi = xi + delta(1);
    yi = yi + delta(2);
    eps = max(abs(f1(xi,yi)),abs(f2(xi,yi)));
    plot(xi,yi,'k+')
    i = i+1;
    fprintf('i= %d, eps = %d\n', i,eps)
end
fprintf('x = %d, y = %d\n',xi,yi)

plot(xi,yi,'r+')
i= 1, eps = 2.080000e+00
i= 2, eps = 1.942885e-01
i= 3, eps = 4.455845e-03
i= 4, eps = 1.330651e-06
i= 5, eps = 1.203482e-13
x = 1.537656e+00, y = 1.906728e+00