gistfile1.matlab

% funkcja pierwotna
f = @(x1,x2) 10*(2*x2.^2-x1).^2+(x2-2).^2;

% gradient
df = @(x) [
  20*x(1) - 40*x(2)^2;
   2*x(2) - 80*x(2)*(x(1) - 2*x(2)^2) - 4
]

% macierz drugich pochodnych
d2f = @(x) [
        20,                  -80*x(2);
  -80*x(2),  480*x(2)^2 - 80*x(1) + 2
]

% aproksymacja kwadratowa funkcji w punkcie X0
af_old = @(x1,x2) (x2-2).^2+10*x1.^2;
X0 = [0;0];
F = f(X0(1), X0(2));
DF = df(X0);
D2F = d2f(X0);
af = @(x) F + DF' * (x-X0) + (x-X0)' * D2F * (x-X0) / 2

[X1, X2] = meshgrid(-5:0.1:5);
VAF = zeros(101)
for x1 = 1:101
    for x2 = 1:101
        %VAF(x2, x1) = af([-5 + 0.1*x1; -5 + 0.1*x2]);
        VAF(x1, x2) = af([X1(x1, x2); X2(x1, x2)]);
        %VAF(x2, x1) = af_old(-5 + 0.1*x1, -5 + 0.1*x2);
    end
end

mesh(X1, X2, f(X1, X2));
hold on;
mesh(X1, X2, VAF)
hold off;