(***********************************************************************)
(* NA3.M                                                               *)
(* Roman Maurer <roman.maurer@fmf.uni-lj.si>, 22-Aug-1996              *)
(* 3. domača naloga iz Numerične analize:                              *)
(* Naloga 2. f: y'=e^x+y, y(0)=1, y(10)=?                              *)
(* Uporabimo:                                                          *)
(* - eksplicitno Adams-Bashforthovo pravilo reda 4,                    *)
(* - implicitno Adams-Moultonovo pravilo reda 4,                       *)
(* - metodo prediktor-korektor s kombinacijo zgornjih dveh pravil.     *)
(*                                                                     *)
(***********************************************************************)

(* y'=f(x,y) *)
f[x_,y_] := N[ Exp[x] + y ]
(* Pozor! Definicija f se uporablja tudi v AdamsMoultonImplicit. *)

(* eksplicitna rešitev *)
g[x_] := N[ (x+1)*Exp[x] ]

(* interval na katerem računamo y *)
l=0;
r=10

(* začetni pogoj y(l)=y0 *)
y0 = 1

(* število delilnih točk maxn = število podintervalov + 1 *)
maxn=200;
h=N[(r-l)/(maxn-1)]


(***********************************************************************)

Print["Roman Maurer, uporabna matematika, 22-Aug-1996"];
Print["3. domača naloga iz Numerične analize:"];
Print["Naloga 2. f: Numerično izračunaj vrednost y(10), če veš, da je"];
Print["y'=e^x+y in y(0)=1."];
Print["Uporabi Adams-Bashforthovo metodo, Adams-Moultonovo metodo in"];
Print["kombinacijo obeh kot prediktor-korektor."];
Print[" "];

RungeKutta::usage = "RungeKutta[xn,yn] vrne približek za naslednji y_{n+1}
z metodo Runge-Kutta reda 4."

RungeKutta[ xn_, yn_ ] := Module[ {k1, k2, k3, k4},
   k1 = h*f[ xn, yn ];
   k2 = h*f[ xn+h/2, yn+k1/2 ];
   k3 = h*f[ xn+h/2, yn+k2/2 ];
   k4 = h*f[ xn+h, yn+k3];
   Return[ yn+(k1+2*k2+2*k3+k4)/6 ] ]


AdamsBashforthExplicit::usage = "AdamsBashforthExplicit[xn, y_List] vrne
naslednji y_{n+1} z eksplicitno Adams-Bashforthovo metodo reda 4.
V seznamu y so zadnje štiri vrednosti y."

AdamsBashforthExplicit[xn_,y_List] :=
   Return[ y[[4]] + h/24 * ( 55*f[ xn, y[[4]] ] - 59*f[ xn-h, y[[3]] ] +
      37*f[ xn-2*h, y[[2]] ] - 9*f[ xn-3*h, y[[1]] ] ) ]


AdamsBashforthAtR::usage = "AdamsBashforthAtR izračuna vrednost y(10)
z Adams-Bashforthovo metodo. Prve štiri približke za y določi z metodo
Runge-Kutta."

AdamsBashforthAtR := Module[ {x, y, ylist, i},
   x=l;  y=y0;
   (* prve štiri y_i naračunamo posebej, i=1,2,3,4 *)
   ylist = {y};
   For[ i=2, i<=4, i++,
      y=RungeKutta[x,y];
      ylist = Append[ ylist, y ];
      x=x+h;
   ];
   (* potem se lotimo metode A-B, i=5,6,7,... *)
   While[ i<=maxn,
      y=AdamsBashforthExplicit[x,ylist];
      x=x+h;
      ylist = ReplacePart[ RotateLeft[ ylist, 1 ], y, 4];
      i++;
   ];
   Return[y]
]


AdamsMoultonImplicit::usage = "AdamsMoultonImplicit[xn,yl_List] vrne
naslednji y, dobljen z implicitno Adams-Moultonovo metodo reda 4. Upošteva,
da velja enačba y'(x)=e^x+y(x). V seznamu yl so zadnje tri vrednosti za y."

AdamsMoultonImplicit[xn_,y_List] := Return[
   (24*y[[3]] + h * (9*Exp[xn+h] + 19*f[ xn, y[[3]] ] - 5*f[ xn-h, y[[2]] ] +
     f[ xn-2*h, y[[1]] ]) ) / (24-9*h) ]


AdamsMoultonAtR::usage = "AdamsMoultonAtR izračuna vrednost y(10)
z implicitno Adams-Moultonovo metodo. Prve tri približke za y določi z metodo
Runge-Kutta."

AdamsMoultonAtR := Module[ {x, y, ylist, i},
   x=l;  y=y0;
   (* prve tri y_i naračunamo posebej, i=1,2,3,4 *)
   ylist = {y};
   For[ i=2, i<=3, i++,
      y=RungeKutta[x,y];
      ylist = Append[ ylist, y ];
      x=x+h;
   ];
   (* potem se lotimo metode A-M, i=4,5,6,... *)
   While[ i<=maxn,
      y=AdamsMoultonImplicit[x,ylist];
      x=x+h;
      ylist = ReplacePart[ RotateLeft[ ylist, 1 ], y, 3];
      i++;
   ];
   Return[y]
]


AdamsMoultonExplicit::usage = "AdamsMoultonExplicit[xn,ylist] daje
vrednost y_{n+1}, dobljeno po metodi korekcije. Predpostavi, da
ima seznam ylist elemente od y_{n-2} do vključno cca. y_{n+1}."

AdamsMoultonExplicit[xn_,y_List] :=
   Return[ y[[3]] + h/24 * ( 9*f[ xn+h, y[[4]] ] +19*f[ xn, y[[3]] ] -
      5*f[ xn-h, y[[2]] ] + f[ xn-2*h, y[[1]] ] ) ]


PredictorCorrector::usage = "PredictorCorrector izračuna vrednost y(10)
z Adams-Bashforth-Moultnovo metodo prediktor-korektor. Prediktor je
Adams-Bashforthovo pravilo reda 4. Korektor je Adams-Moultnovo pravilo
reda 4."

PredictorCorrector := Module[ {x, y, ylist, yl2, i},
   x=l;  y=y0;
   (* prve štiri y_i naračunamo posebej, i=1,2,3,4 *)
   ylist = {y};
   For[ i=2, i<=4, i++,
      y=RungeKutta[x,y];
      ylist = Append[ ylist, y ];
      x=x+h;
   ];
   (* potem se lotimo metode A-B-M prediktor-korektor, i=5,6,7,... *)
   While[ i<=maxn,
      y=AdamsMoultonExplicit[x, ReplacePart[ RotateLeft[ ylist, 1 ],
          AdamsBashforthExplicit[x,ylist], 4] ];
      ylist = ReplacePart[ RotateLeft[ ylist, 1 ], y, 4];
      x=x+h;
      i++;
   ];
   Return[y]
]


Print["Adams-Bashforth:    ", Reverse[ AdamsBashforthAtR // Timing ]];
Print["Adams-Moulton:      ", Reverse[ AdamsMoultonAtR // Timing ]];
Print["prediktor-korektor: ", Reverse[ PredictorCorrector // Timing ]];
Print[" "];
Print["eksplicitna rešitev: y(10)=",g[10]];

(* konec datoteke NA3.M *)