Brent[f_, AX_, BX_, CX_, TOL_, ITMAX_]:=

(* \(Treba\ \(zadati : 
Funkciju\ "f"\)\),
Tri početne aproksimacije "AX<BX<CX",ali tako da bude "f[BX]<Min[AX,CX]",
Točnost "TOL",
Maksimalni broj iteracija "ITMAX",
\(Rjeenje\ je\ XMIN\) *)

Module[{ CGOLD =  .5 (3 - Sqrt[5]),  
              ZEPS = N[10^(-12)],  A, B, V, W, X, U, e, FX, FV, FW, 
              FU, XM, R, Q, P, d, n1, n2},
A = Min[AX, CX];  
          B = Max[AX, CX];  V = BX;  W = V;   X = V; e = 0. ; 
          FX = f[X];  FV = FX;  FW = FX; 
          ITER = 0; 
While[XM = (A + B)/2; 
            TOL1 = TOL Abs[X] + ZEPS; TOL2 = 2  TOL1;  
            Abs[X - XM] > TOL2 - (B - A)/2 && ITER < ITMAX,  
            ITER = ITER + 
                1; 
(
If[
                  Abs[e] > TOL1,  R = (X - W)  (FX - FV); 
                  Q = (X - V) (FX - FW); 
                  P = (X - V)  Q - (X - W)  R; 
                  Q = 2  (Q - R); 
   If[Q > 0, 
                    P = -P]; Q = Abs[Q]; ETEMP = e; 
                  e = d; 
   If[
                    Abs[P] > Abs[ .5 Q ETEMP] || 
                      P ≤ Q  (A - X) || 
                      P ≥ Q (B - X), 
                    Goto[n1]]; 
d = P/Q;  
                  U = X + d; 
     If[
                    U - A < TOL2 || B - U < TOL2, 
                    If[XM == X, d = A bs[TOL1], 
                      d = Sign[XM - X] A bs[TOL1]]]; 
                  Goto[n2];
]; 
Label[
                  n1]; 
If[X ≥ XM, e = A - X, 
                  e = B - X]; 
                d = CGOLD  e; 
Label[
                  n2]; 
If[d == 0, XX = Abs[TOL1], 
                  XX = Abs[TOL1] Sign[d]]; 
If[
                  Abs[d] > TOL1,  U = X + d, 
                  U = X + XX]; 
FU = 
                  f[U]; 
)
If[
                FU ≤ FX, 
    If[
                  U ≥ X, A = X, B = X]; V = W; FV = FW; W = X; 
                FW = FX; X = U; 
                FX = FU, 
    If[U < X, A = U, B = U]; 
                If[FU ≤ FW || W == X, V = W; FV = FW; W = U; 
                  FW = FU, 
                  If[FU ≤ FV || V == X || V == W, 
                    V = U; FV = 
                      FU]]
    ]
]; 
          XMIN = X; 
If[ITER ≥ ITMAX, 
            Print["Brent Exceed Maximum Iteration"]]
]

NEWTONOVA METODA ZA JEDNODIMENZIONALNU MINIMIZACIJU

\(\(\(
\)\(NewtonMin1[f_, x0_, eps_, it_] := 
      Module[{a = N[x0], 
          k = 0}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ While[\ 
          Abs[\(f'\)[a]] > eps\  && \ 
            k < it, \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Print["x", k, "=", a, "   f", k, "=", 
            f[a]]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ a = 
            a - \(f'\)[a]/\(f''\)[
                  a]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ k = 
            k + 1;
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ a]\)\)\)

f[x_] := 
        10 - 2  x + 4  x^2 - 2  x^3 + 
          0.25  x^4; 
x0 = 4;  eps =  .000005;  it = 10;  
    NewtonMin1[f, x0, eps, it]

FindMinimum[f[x], {x, x0}]

MINIMIZACIJA funkcija vise varijabli

f[x_] := 
      Exp[x[1]] (4  x[1]^2 + 2  x[2]^2 + 4  x[1] x[2] + 2  x[2] + 1);  
    n = 2;
g[x_]:=Table[D[f[x], x[i]],{i, n}]
H[x_]:=Table[D[g[x][[i]], x[j]],{i, n},{j, n}]
Print["grad f=",  
      Simplify[g[x] /. {x[1] -> x1, x[2] -> x2}] // 
        MatrixForm,  "   H(f)=",  
      Simplify[H[x] /. {x[1] -> x1, x[2] -> x2}] // 
        MatrixForm]

FindMinimum[
      f[x] /. {x[1] -> x1, x[2] -> x2}, {x1, 0}, {x2, 
        0}]

g[x] /. {x[1] ->  .5, 
        x[2] -> -1}
H[x] /. {x[1] ->  .5, 
        x[2] -> -1}
Det[H[x] /. {x[1] ->  .5, x[2] -> -1}]

ContourPlot[
        f[x] /. {x[1] -> x1, x[2] -> x2}, {x1, -2.5, 
          1}, {x2, -2, 1.5}, Contours -> 40,  
        PlotPoints -> 30, ContourShading -> None];

GRADIJENTNA METODA

Clear[f, a, x1, x2];  n = 2;
f[x_] := 
        100 (x[2] - x[1]^2)^2 + (1 - 
                x[1])^2 

        
  
 d[x_]:=Table[D[f[x], x[i]],{i, n}]
d[x] /. {x[1] -> x1, x[2] -> x2} 
ContourPlot[
        f[x] /. {x[1] -> x1, x[2] -> x2}, {x1, -2, 
          2}, {x2, -2, 4}, Contours -> 20, 
        ContourShading -> None];

Korak[n_, f_, pk_, xk_]:=Module[{epsilon =  .5, h =  .5, limit =  .000001, f0, f1, fd},
al = 1; w = Table[x[i] -> xk[[i]], {i, n}];
                                     d[x_]:=Table[D[f[x], x[i]],{i, n}];
                                     f0 =  f[x] /. w;
                                     fd =  d[x] /. w;
                    While[
              f1 =  f[x] /. 
                  Table[x[i] -> xk[[i]] + al pk[[i]], {i, 
                      n}];       
         f1 - f0 > epsilon al Apply[Plus,  fd pk]  &&  al > limit, al = h al];
If[al < limit, al = 1];
al]

MinGrad[n_, f_, x0_, eps_, it_]:=Module[{k = 0, xs = x0, xn, gr, w},w = Table[x[i] -> x0[[i]], {i, n}];
                d[x_]:=Table[D[f[x], x[i]],{i, n}];
While[ 
              w = Table[x[i] -> xs[[i]], {i, n}];   grad = 
                Max[Abs[d[x] /. w]]; 
grad > eps  &&  
                k <  it, 
pk = -d[x] /. 
                  w;  
alpha = 
                Korak[n, f, pk, 
                  xs]; 
  
xn = 
                xs + alpha pk; 
xs = xn; 
              k = k + 1;
];
{ k, alpha, xn, 
              f[x] /. Table[x[i] -> xn[[i]], {i, n}], grad}]

x0 = { .5,  .5}; MinGrad[n, f, x0,  .005, 1150]

METODA NAJBRZEG SPUSTA

f[x_] := 100 (x[2] - x[1]^2)^2 + (1 - x[1])^2 ;  
    n = 2;
f[x] /. {x[1] -> x1, 
        x[2] -> x2}
p[x_]:=-Table[D[f[x], x[i]],{i, n}]
p[x] /. {x[1] -> x1, 
        x[2] -> x2}

x0 = { .5,  .5};
Print["x0= ", MatrixForm[x0], ",  f[x]= ", 
        f[x] /. Table[
            x[i] -> x0[[i]], {i, n}]];
Do[
korak[al_] := 
        Module[{pp}, 
          pp = p[x] /. Table[x[i] -> x0[[i]], {i, n}];     f[x] /. 
            Table[x[i] -> x0[[i]] + al pp[[i]]  , {i, 
                n}] ]; 
Brent[korak,  .0001, 1,  .5,  .01, 
        50]; alpha = XMIN;  
pp = 
        p[x] /. Table[x[i] -> x0[[i]], {i, n}]; 
      npp = Sqrt[Sum[pp[[i]]^2, {i, n}]]; 
x1 = 
        x0 + alpha pp/npp; 
      x0 = x1;  
w = 
        Table[x[i] -> x0[[i]], {i, n}];  
      fx = f[x] /. w;   grad = Max[Abs[p[x] /. w]], {k, 
        1000}]
Print["x", k, "= ", x1, ",  f[x]= ", fx, "  grad=", 
      grad, "  alpha=", alpha]

Newtonova metoda

MinNewton[n_, f_, x0_, eps_, it_]:=Module[{k = 0, xs = x0, xn, gr, w},
                     w = Table[x[i] -> x0[[i]], {i, n}];
                     d[x_]:=Table[D[f[x], x[i]],{i, n}];
                     H[x_]:=Table[D[d[x][[i]], x[j]],{i, n},{j, n}];
                     Print["x0= ", x0, "   f0= ", 
              f[x] /. w, "  grad=", Max[Abs[d[x] /. w]]];
       \(While[\ 
              w = Table[x[i] -> xs[\([i]\)], {i, n}]; \ \ grad = 
                Max[Abs[
                    d[x] /. 
                      w]]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ grad > eps\  && \ k < \ 
                  it, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ pk = \(-LinearSolve[H[x] /. w, \ 
                    d[x] /. 
                      w]\); 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ alpha = Korak[n, f, pk, 
                  xs]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ xn = 
                xs + alpha\ pk; xs = xn; 
              k = k + 1; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ w = 
                Table[x[i] -> xn[\([i]\)], {i, n}]; \ \ grad = 
                Max[Abs[d[x] /. w]]; 
Print["x", 
                k, "= ", xn, ",  f[x]= ", f[x] /. w, "  grad=", 
                grad, "  alpha=", alpha // N]
\ ]\)
]

n = 2;
 f[x_] := 100 (x[2] - x[1]^2)^2 + (1 - x[1])^2
  
f[x] /. {x[1] -> x1, x[2] -> x2}
d[x_]:=Table[D[f[x], x[i]],{i, n}]
d[x] /. {x[1] -> x1, x[2] -> x2}
H[x_]:=Table[D[d[x][[i]], x[j]],{i, n},{j, n}]
H[x] /. {x[1] -> x1, x[2] -> x2}
x0 = { .5,  .5}; 
MinNewton[n, f, x0,  .01, 20]

Quasi - Newtonove metode

24.  05.  2005.

Duljina koraka

Korak[n_, f_, pk_, xk_]:=Module[{epsilon =  .5, h =  .5, limit =  .000001, f0, f1, fd},
al = 1; w = Table[x[i] -> xk[[i]], {i, n}];
                                     d[x_]:=Table[D[f[x], x[i]],{i, n}];
                                     f0 =  f[x] /. w;
                                     fd =  d[x] /. w;
                    While[
              f1 =  f[x] /. 
                  Table[x[i] -> xk[[i]] + al pk[[i]], {i, 
                      n}];       
         f1 - f0 > epsilon al Apply[Plus,  fd pk]  &&  al > limit, al = h al];
If[al < limit, al = 1];
al]

f[x_] := 
      Exp[x[1]] (4  x[1]^2 + 2  x[2]^2 + 4  x[1] x[2] + 2  x[2] + 1);  
    n = 2;
g[x_]:=Table[D[f[x], x[i]],{i, n}]
H[x_]:=Table[D[g[x][[i]], x[j]],{i, n},{j, n}]
Print["grad f=",  
      Simplify[g[x] /. {x[1] -> x1, x[2] -> x2}] // 
        MatrixForm,  "   H(f)=",  
      Simplify[H[x] /. {x[1] -> x1, x[2] -> x2}] // 
        MatrixForm]
ContourPlot[
        f[x] /. {x[1] -> x1, x[2] -> x2}, {x1, -2.5, 
          1}, {x2, -2, 1.5}, Contours -> 40,  
        PlotPoints -> 30, ContourShading -> None];

FindMinimum[
      f[x] /. {x[1] -> x1, x[2] -> x2}, {x1, 0}, {x2, 
        0}]

n = 2; x0 = {1. , 1};  k = 0; eps = 0.005;
 w = Table[x[i] -> x0[[i]], {i, n}]; 
    B0 = IdentityMatrix[n];
B0 = H[x]/.w 
g0 = g[x] /. w

Direktna QN metoda

Do[
p0 = -LinearSolve[B0, 
            g0]; 
alpha = 
        Korak[n, f, p0, x0]; 
x1 = 
        x0 + alpha p0; 
 
      w1 = Table[x[i] -> x1[[i]], {i, n}];  
      g1 = g[x] /. w1; 
y0 = 
        g1 - g0; 
  
B1 = 
        B0 + Outer[
              Times, (y0 - B0 . p0), (y0 - 
                  B0 . p0)]/((y0 - B0 . p0)  .  
                p0); 
B0 = B1; 
      g0 = g1; 
x0 = x1; 
      k = k + 1; 
Print["x", k, "=", 
        x1, "  fx=", f[x] /. w1, "  grad=", g1, " ", alpha], {k, 
        3}]

While[
p0 = -LinearSolve[B0, 
              g0]; 
alpha = 
          Korak[n, f, p0, x0]; 
x1 = 
          x0 + alpha p0; 
 
        w1 = Table[x[i] -> x1[[i]], {i, n}];  
        g1 = g[x] /. w1; 
y0 = 
          g1 - g0; 
  
B1 = 
          B0 + Outer[
                Times, (y0 - B0 . p0), (y0 - 
                    B0 . p0)]/((y0 - B0 . p0)  .  
                  p0); 
f1 = 
          f[x] /. w1; 
Abs[f1] > 
          eps, 
B0 = B1; 
        g0 = g1; 
x0 = x1; 
        k = k + 1; 
Print["x", k, "=", 
          x1, "  fx=", f[x] /. w1, "  grad=", g1, "  alpha=", 
          alpha]
];
Print["x", k + 1, "=", x1, "  fx=", 
      f[x] /. w1, "  grad=", g1, "  alpha=", alpha]

g[x] /. w1
h = H[x] /. w1
Det[h]