(* Bor Plestenjak     *)
(* IMFM/TCS Ljubljana *)

(* Programi za laboratorijske vaje pri NLA,NM                          *)

(* PODROCJE : resevanje nelinearnih sistemov                           *)
(* SNOV : bisekcija, navadna iteracija, tangentna metoda, modificirana
          tangentna metoda, sekantna metoda, Laguerreova metoda        *)

(* =========================== KONSTANTE ================================== *)

dec=20;             (* racunanje na dec mest        *)
decizpis=10;        (* izpisovanje na decizpis mest *)

(* =========================== NAVODILA =================================== *)

Bisection::usage="Bisection[f,x0,x1,maxit,nat] je metoda bisekcije za iskanje 
  nicel funkcije f spremenljivke x z zacetnima priblizkoma x0 in x1.
  Bisection[f,x0,x1,maxit,nat,spr] je metoda bisekcije za iskanje nicel 
  funkcije f spremenljivke spr. Maksimalno stevilo iteracij je maxit, 
  zahtevana natancnost pa nat."

Bisection1::usage="Ukaz Bisection1 je enak ukazu Bisection, le da vmes 
  izpisuje se tekoce priblizke."

Iterate::usage="Iterate[f,x0,maxit,nat] iterira funkcijo f spremenljivke
  x z zacetnim priblizkom x0. Iterate [f,x0,maxit,nat,spr] iterira funkcijo
  f spremenljivke spr z zacetnim priblizkom x0. Iterate[{f,g},{x0,y0},maxit,nat] iterira par funkcij
  f,g spremenljivk x,y z zacetnim priblizkom {x0,y0}. Iterate[{f,g},{x0,y0},maxit,nat,{s1,s2}] iterira par funkcij
  f,g spremenljivk {s1,s2} z zacetnim priblizkom {x0,y0}."

Iterate1::usage="Ukaz Iterate1 je enak ukazu Iterate, le da vmes izpisuje se
  tekoce priblizke."

Newton::usage="Newton[f,x0,maxit,nat] je tangentna metoda za iskanje nicel
  funkcije f spremenljivke x z zacetnim priblizkom x0.
  Newton[f,x0,maxit,nat,spr] je tangentna metoda za iskanje nicel funkcije f
  spremenljivke spr. Newton[{f,g},{x0,y0},maxit,nat] je Newtonova metoda v
  dveh dimenzijah za iskanje nicel sistema funkcij f,g spremenljivk x,y.
  Newton[{f,g},{x0,y0},maxit,nat,{s1,s2}] je Newtonova metoda v dveh
  dimenzijah za iskanje nicel sistema funkcij f,g spremenljivk s1,s2.
  Maksimalno stevilo iteracij je maxit, zahtevana natancnost pa nat."

Newton1::usage="Ukaz Newton1 je enak ukazu Newton, le da se med racunanjem
  izpisujejo tekoci priblizki in vrednosti funkcije."

ModNewton::usage="ModNewton[f,x0,maxit,nat,k] je modificirana tangentna metoda
  za iskanje k-kratnih nicel funkcije f spremenljivke x z zacetnim priblizkom
  x0. ModNewton[f,x0,maxit,nat,k,spr] je modificirana tangentna metoda za
  iskanje k-kratnih nicel funkcije f spremenljivke spr z zacetnim priblizkom
  x0. Maksimalno stevilo iteracij je maxit, zahtevana natancnost pa nat."

ModNewton1::usage="ModNewton1 je enak ukazu ModNewton, le da se vmes
  izpisujejo vmesni priblizki in vrednosti funkcije."

Secant::usage="Secant[f,x0,x1,maxit,nat] je sekantna metoda za iskanje nicel
  funkcije f spremenljivke x z zacetnima priblizkoma x0 in x1.
  Secant[f,x0,x1,maxit,nat,spr] je sekantna metoda za iskanje nicel funkcije
  f spremenljivke spr. Maksimalno stevilo iteracij je maxit, zahtevana
  natancnost pa nat."

Secant1::usage="Ukaz Secant1 je enak ukazu Secant, le da se med racunanjem
  izpisujejo tekoci priblizki in vrednosti funkcije."

Laguerre::usage="Laguerre[f,x0,n,maxit,nat,izb] je Laguerrova metoda za
  iskanje nicel funkcije f spremenljivke x z zacetnim priblizkom x0 in
  konstanto n. Laguerre[f,x0,maxit,nat,izb,spr] je Laguerrova metoda za
  iskanje nicel funkcije f spremenljivke spr z zacetnim priblizkom x0 in
  konstanto n. Maksimalno stevilo iteracij je maxit, zahtevana natancnost
  pa nat. Izbira je lahko: 1 za minus, 2 za plus in 3 za stabilno varianto."

LaguerreM::usage="LaguerreM[f,x0,n,maxit,nat] je Laguerrova metoda za
  iskanje nicel funkcije f spremenljivke x z zacetnim priblizkom x0 in
  konstanto n. LaguerreM[f,x0,n,maxit,nat,spr] je Laguerrova metoda za
  iskanje nicel funkcije f spremenljivke spr z zacetnim priblizkom x0 in
  konstanto n. Maksimalno stevilo iteracij je maxit, zahtevana natancnost pa
  nat. V imenovalcu konstantno uporablja -."

LaguerreP::usage="LaguerreP[f,x0,n,maxit,nat] je Laguerrova metoda za
  iskanje nicel funkcije f spremenljivke x z zacetnim priblizkom x0 in
  konstanto n. LaguerreP[f,x0,n,maxit,nat,spr] je Laguerrova metoda za
  iskanje nicel funkcije f spremenljivke spr z zacetnim priblizkom x0 in
  konstanto n. Maksimalno stevilo iteracij je maxit, zahtevana natancnost pa
  nat. V imenovalcu konstantno uporablja +."

LaguerreS::usage="LaguerreS[f,x0,n,maxit,nat] je stabilna Laguerrova metoda za
  iskanje nicel funkcije f spremenljivke x z zacetnim priblizkom x0 in
  konstanto n. LaguerreS[f,x0,n,maxit,nat,spr] je Laguerrova metoda za
  iskanje nicel funkcije f spremenljivke spr z zacetnim priblizkom x0 in
  konstanto n. Maksimalno stevilo iteracij je maxit, zahtevana natancnost
  pa nat."

(* =========================== PROGRAMI =================================== *)

(* --------------------------- BISEKCIJA ---------------------------------- *)

Bisection[f_,x0_,x1_,maxit_,nat_,izpis_,spr_]:=Module[{a,b,c,fa,fb,fc,step,g},
  g[z_]=f /. spr->z;
  a=N[x0]; b=N[x1]; c=(a+b)/2;
  fa=g[a]; fb=g[b]; fc=g[c];
  If[fa fb>0,
    (* then *) Print["f(a) :",fa];
               Print["f(b) :",fb];
               Message[Bisection::wrdata],
    (* else *) step=1;
               If[izpis,Print[step,":  ",N[c,decizpis]]];
               While[fc!=0 && Abs[b-a]>nat && step<maxit,
                 If[N[fa fc]>0,
                   (* then *) {a,fa}={c,fc},
                   (* else *) {b,fb}={c,fc}];
                  c=(a+b)/2;
                  fc=g[c];
                  step=step+1;
                  If[izpis,Print[step,":  ",c,"   ",fc]] 
                  ]; (* end while *)
    If[step==maxit,Message[Bisection::slowitr]
    ]; (* end if *)
  c ] ];

Bisection[f_,x0_,x1_,maxit_,nat_,spr_:x]:=
  Bisection[f,x0,x1,maxit,nat,False,spr]; 
Bisection1[f_,x0_,x1_,maxit_,nat_,spr_:x]:=
  Bisection[f,x0,x1,maxit,nat,True,spr]; 

Bisection::wrdata="Funkcija je v obeh krajiscih intervala enako predznacena."

Bisection::slowitr="Premajhno stevilo iteracij ali prevelik zacetni interval!";

(* --------------------- NAVADNA ITERACIJA -------------------------------- *)

Iterate[f_,x0_,maxit_,nat_,izpis_,spr_]:=Module[{a,b,step,g},
  g[z_]=f /. spr->z;
  a=SetPrecision[N[x0,dec],dec];
  b=SetPrecision[N[g[a],dec],dec];
  step=1;
  If[izpis,Print[step,":  ",N[b,decizpis]]];
  While[Abs[b-a]>nat && step<maxit,
    a=SetPrecision[N[b,dec],dec];
    b=SetPrecision[N[g[a],dec],dec];
    step=step+1;
    If[izpis,Print[step,":  ",N[b,decizpis]]] ];
  If[step==maxit,
    Message[Iterate::slowitr]];
  b ];

Iterate2[{f_,g_},{x0_,y0_},maxit_,nat_,izpis_,spr_]:=Module[
  {a,b,step,h},
  h[z_,w_]={f,g} /. {spr[[1]]->z,spr[[2]]->w};
  a=SetPrecision[N[{x0,y0},dec],dec];
  b=SetPrecision[N[h[a[[1]],a[[2]]],dec],dec];
  step=1;
  If[izpis,Print[step,":  ",N[b,decizpis]]];
  While[Max[Abs[b-a]]>nat && step<maxit,
    a=SetPrecision[N[b,dec],dec];
    b=SetPrecision[N[h[a[[1]],a[[2]]],dec],dec];
    step=step+1;
    If[izpis,Print[step,":  ",N[b,decizpis]]] ];
  If[step==maxit,
    Message[Iterate::slowitr]];
  b ];

Iterate[f_,x0_,maxit_,nat_,spr_:x]:=
  Iterate[f,x0,maxit,nat,False,spr] /; Not[ListQ[f]];
Iterate1[f_,x0_,maxit_,nat_,spr_:x]:=
  Iterate[f,x0,maxit,nat,True,spr] /; Not[ListQ[f]];

Iterate[{f_,g_},{x0_,y0_},maxit_,nat_,spr_:{x,y}]:=
  Iterate2[{f,g},{x0,y0},maxit,nat,False,spr];
Iterate1[{f_,g_},{x0_,y0_},maxit_,nat_,spr_:{x,y}]:=
  Iterate2[{f,g},{x0,y0},maxit,nat,True,spr]; 

Iterate::slowitr = "Premajhno stevilo iteracij ali pa slab priblizek !";

(* -------------------- TANGENTNA METODA ---------------------------------- *)

NewtonIt[f_,x0_,maxit_,nat_,izpis_,spr_]:=Module[{Newt,a,b,step,g},
  g[z_]=f /. spr->z;
  Newt[z_]=z-g[z]/D[g[z],z];
  a=SetPrecision[N[x0,dec],dec];
  b=SetPrecision[N[Newt[a],dec],dec];
  step=1;
  If[izpis,Print[step,":  ",N[b,decizpis],"   ",N[g[b],decizpis]]];
  While[N[g[b]]!=0 && Abs[b-a]>nat && step<maxit,
    a=SetPrecision[N[b,dec],dec];
    b=SetPrecision[N[Newt[a],dec],dec];
    step=step+1;
    If[izpis,Print[step,":  ",N[b,decizpis],"    ",N[g[b],decizpis]]] ];
  If[step==maxit,
    Message[NewtonIt::slowitr]];
  b ];

Newton2[{f_,g_},{x0_,y0_},maxit_,nat_,izpis_,{s1_,s2_}]:=Module[
  {Newt,a,b,step,h1,h2,Jac,Inv},
  h1[{z_,w_}]=f /. {s1->z,s2->w};
  h2[{z_,w_}]=g /. {s1->z,s2->w};
  Jac[{z_,w_}]={{D[f,s1],D[f,s2]},{D[g,s1],D[g,s2]}} /. {s1->z,s2->w};
  Inv[par_]:=Inverse[Jac[par]];
  Newt[par_]=par-Inv[par].{h1[par],h2[par]};
  a=SetPrecision[N[{x0,y0},dec],dec];
  b=SetPrecision[N[Newt[a],dec],dec];
  step=1;
  If[izpis,Print[step,":  ",N[b,decizpis],"   ",N[{h1[b],h2[b]},decizpis]]];
  While[N[{h1[b],h2[b]}]!={0,0} && Max[Abs[b-a]]>nat && step<maxit,
    a=SetPrecision[N[b,dec],dec];
    b=SetPrecision[N[Newt[a],dec],dec];
    step=step+1;
    If[izpis,Print[step,":  ",N[b,decizpis],"    ",N[{h1[b],h2[b]},decizpis]]] ];
  If[step==maxit,
    Message[NewtonIt::slowitr]];
  b ];

NewtonIt[f_,x0_,maxit_,nat_,spr_:x]:=
  NewtonIt[f,x0,maxit,nat,False,spr] /; Not[ListQ[f]];

NewtonIt[{f_,g_},{x0_,y0_},maxit_,nat_,spr_:{x,y}]:=
  Newton2[{f,g},{x0,y0},maxit,nat,False,spr];

Newton1[f_,x0_,maxit_,nat_,spr_:x]:=	
  NewtonIt[f,x0,maxit,nat,True,spr] /; Not[ListQ[f]];

Newton1[{f_,g_},{x0_,y0_},maxit_,nat_,spr_:{x,y}]:=
  Newton2[{f,g},{x0,y0},maxit,nat,True,spr];

NewtonIt::slowitr = "Premajhno stevilo iteracij ali pa slab priblizek !";

(* ----------------- MODIFICIRANA TANGENTNA METODA ------------------------ *)

ModNewton[f_,x0_,maxit_,nat_,k_,spr_,izpis_]:=Module[{Newt,a,b,step,g},
  g[z_]=f /. spr->z;
  Newt[z_]=z-k g[z]/D[g[z],z];
  a=SetPrecision[N[x0,dec],dec];
  b=SetPrecision[N[Newt[a],dec],dec];
  step=1;
  If[izpis,Print[step,":  ",N[b,decizpis],"   ",N[g[b],decizpis]]];
  While[N[g[b]]!=0 && Abs[b-a]>nat && step<maxit,
    a=SetPrecision[N[b,dec],dec];
    b=SetPrecision[N[Newt[a],dec],dec];
    step=step+1;
    If[izpis,Print[step,":  ",N[b,decizpis],"    ",N[g[b],decizpis]]] ];
  If[step==maxit,
    Message[Iteracija::slowitr]];
  b ];

ModNewton[f_,x0_,maxit_,nat_,k_,spr_:x]:=
  ModNewton[f,x0,maxit,nat,k,spr,False];

ModNewton1[f_,x0_,maxit_,nat_,k_,spr_:x]:=
  ModNewton[f,x0,maxit,nat,k,spr,True];

NewtonIt::slowitr = "Premajhno stevilo iteracij ali pa slab priblizek !";

(* -------------------- SEKANTNA METODA ----------------------------------- *)

Secant[f_,x0_,x1_,maxit_,nat_,spr_,izpis_]:=Module[{Sec,a,b,c,step,g},
  g[z_]=f /. spr->z;
  Sec[z_,w_]=z-g[z] (z-w)/(g[z]-g[w]);
  a=SetPrecision[N[x0,dec],dec];
  b=SetPrecision[N[x1,dec],dec];
  c=SetPrecision[N[Sec[b,a],dec],dec];
  step=1;
  If[izpis,Print[step,":  ",N[c,decizpis],"   ",N[g[c],decizpis]]];
  While[N[g[c]]!=0 && Abs[c-b]>nat && step<maxit,
    a=SetPrecision[N[b,dec],dec];
    b=SetPrecision[N[c,dec],dec];
    c=SetPrecision[N[Sec[b,a],dec],dec];
    step=step+1;
    If[izpis,Print[step,":  ",N[c,decizpis],"    ",N[g[c],decizpis]]] ];
  If[step==maxit,
    Message[Secant::slowitr]];
  c ];

Secant[f_,x0_,x1_,maxit_,nat_,spr_:x]:=Secant[f,x0,x1,maxit,nat,spr,False];

Secant1[f_,x0_,x1_,maxit_,nat_,spr_:x]:=Secant[f,x0,x1,maxit,nat,spr,True];

Secant::slowitr="Premajhno stevilo iteracij ali pa slaba zacetna priblizka !";

(* -------------------- LAGUERROVA METODA --------------------------------- *)

Laguerre[f_,x0_,n_,maxit_,nat_,izbira_,spr_]:=Module[{g,Lag,a,b,step},
  g[z_]=f /. spr->z;
  Lag[x_]:=Module[{ima,imb},
    If[g[x]!=0,
    ima=SetPrecision[N[g'[x]/g[x],dec],dec];
    imb=SetPrecision[N[Sqrt[(n-1)((n-1) (g'[x]/g[x])^2-n g''[x]/g[x])],dec],dec];
    If[izbira==1 || (izbira==3 && Abs[ima-imb]>Abs[ima+imb]),
      SetPrecision[N[x-n /(ima-imb),dec],dec],
      SetPrecision[N[x-n /(ima+imb),dec],dec] ],
      x] ];
  a=SetPrecision[N[x0,dec],dec];
  b=SetPrecision[N[Lag[a],dec],dec];
  step=1;
  Print[step,":  ",N[b,decizpis],"   ",N[g[b],decizpis]];
  While[N[g[b]]!=0 && Abs[b-a]>nat && step<maxit,
    a=SetPrecision[N[b,dec],dec];
    b=SetPrecision[N[Lag[a],dec],dec];
    step=step+1;
    Print[step,":  ",N[b,decizpis],"    ",N[g[b],decizpis]] ];
  If[step==maxit,
    Message[Laguerre::slowitr]];
  b ];

LaguerreM[f_,x0_,n_,maxit_,nat_,spr_:x]:=Laguerre[f,x0,n,maxit,nat,1,spr];

LaguerreP[f_,x0_,n_,maxit_,nat_,spr_:x]:=Laguerre[f,x0,n,maxit,nat,2,spr];

LaguerreS[f_,x0_,n_,maxit_,nat_,spr_:x]:=Laguerre[f,x0,n,maxit,nat,3,spr];

Laguerre::slowitr = "Premajhno stevilo iteracij ali pa slab priblizek !";

(* ============================ SPISEK UKAZOV ============================= *)

Print["Ukazi z metodami za iskanje nicel nelinearnih enacb so:"]
Print[""]
Print["   - metoda bisekcije"]
Print["         - Bisection, Bisection1"]
Print[""]
Print["   - navadna iteracija v eni in dveh spremenljivkah"]
Print["         - Iterate, Iterate1"]
Print[""]
Print["   - tangentna metoda v eni in dveh spremenljivkah"]
Print["         - NewtonIt, Newton1, ModNewton, ModNewton1"]
Print[""]
Print["   - sekantna metoda"]
Print["         - Secant, Secant1"]
Print[""]
Print["   - Laguerrova metoda"]
Print["         - Laguerre, LaguerreM, LaguerreP, LaguerreS"]