(***********************************************************************)
(* NA4.M                                                               *)
(* Roman Maurer <roman.maurer@fmf.uni-lj.si>, 29-Aug-1996              *)
(* 4. domača naloga iz Numerične analize:                              *)
(* Naloga 1. c:   u_t = u_xx,  0 < x < 1,  t > 0                       *)
(*      u(x,0) = sin(pi*x),   0 <= x <= 1                              *)
(*      u(0,t) = u(1,t) = 0,  t >= 0                                   *)
(* Uporabimo:                                                          *)
(* - eksplicitno metodo (dx=0.1, lambda=1/4),                          *)
(* - implicitno metodo (dx=0.1, lambda=1/4),                           *)
(* - Crank-Nicholsonovo metodo (dx=0.1, lambda=1).                     *)
(*                                                                     *)
(***********************************************************************)

(* u(x,0)=f(x) *)
f[x_] :=  Sin[Pi*x]

(* u(0,t)=g(t), u(1,t)=h(t) *)
(* g[t_] := 0;
h[t_] := 0 *)


(* x teče s korakom dx *)
dx=1/10

(* lambda = dt/dx^2 *)
lambda=1/4;
(* t teče od 0 do 1 s K koraki po dt *)
dt = lambda*dx^2;
K=Round[1/dt]

(* točke za x so {j*dx; j=0,1,...,J+1} *)
J = 1/dx - 1


u0::usage = "u0 je vrednost funkcije u(x,t) v notranjih točkah za t=0,
x=j*dx, j=1,..,J"

u0 = N[ Table[ f[j*dx], {j, J} ] ]


u1::usage = "u1 je analitično izračunana vrednost konkretne funkcije u(x,t)
v notranjih točkah za t=1, x=j*dx, j=1,..,J
(potrebujemo za končno primerjavo primernosti metod)"

u1 = N[u0*Exp[-Pi^2]]


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


InverseTridiagForImpl::usage = "InverseTridiagForImpl je inverz
tridiagonalne matrike, ki ga potrebujemo pri implicitni metodi.
Izračuna ga funkcija ImplicitAt1."


TridiagForExpl::usage = "TridiagForExpl je tridiagonalna matrika,
ki jo potrebujemo za eksplicitno metodo. Izračuna jo funkcija
ExplicitAt1."


CrankNicholsonMatrix::usage = "CrankNicholsonMatrix je matrika 9x9
za Crank-Nicholsonovo iteracijo za lambda=1."



ImplicitStep::usage = "ImplicitStep[u_List] predpostavi, da vsebuje u
vrednosti u(j*dx,t), j=1,...,J. Po implicitni metodi določi
vrednosti u(j*dx,t+dt), j=1,...,J. Bistveno upošteva, da je g(t)=h(t)=0."
(* V konkretnem primeru za g(t)=h(t)=0 je dovolj rešiti linearni sistem:
        T.u_(n+1) = u_n
   V splošnem je treba na desni strani enačbe dodati še
        + lambda*[ g(t+dt), 0, 0, ... ,0, 0, h(t+dt) ]  *)

ImplicitStep[u_List] := N[ InverseTridiagForImpl . u]


ImplicitAt1::usage = "ImplicitAt1 po implicitni metodi določi notranje
vrednosti u(j*dx,t=1), j=1,...,J. Vrne seznam s temi vrednostmi.
Če je metoda dobra, mora biti seznam blizu u1."

ImplicitAt1 := Module[ {},
   InverseTridiagForImpl = N[ Inverse[ Table[
      Switch[ i-j, -1, -lambda, 0, 1+2*lambda, 1, -lambda, _, 0],
      {i, J}, {j, J} ] ] ];
   Nest[ ImplicitStep, u0, K ]
]


ExplicitStep::usage = "ExplicitStep[u_List] predpostavi, da vsebuje u
vrednosti u(j*dx,t), j=1,...,J. Po eksplicitni metodi določi
vrednosti u(j*dx,t+dt), j=1,...,J. Bistveno upošteva, da je g(t)=h(t)=0."

ExplicitStep[u_List] := TridiagForExpl . u


ExplicitAt1::usage = "ExplicitAt1 po eksplicitni metodi določi notranje
vrednosti u(j*dx,t=1), j=1,...,J. Vrne seznam s temi vrednostmi.
Če je metoda dobra, mora biti seznam blizu u1."

ExplicitAt1 := Module[ {},
   TridiagForExpl = N[ Table[
      Switch[ i-j, -1, lambda, 0, 1-2*lambda, 1, lambda, _, 0],
      {i, J}, {j, J} ] ];
   Nest[ ExplicitStep, u0, K ]
]


CrankNicholsonStep::usage = "CrankNicholsonStep[u_List] predpostavi,
da vsebuje u vrednosti u(j*dx,t), j=1,...,J. Po Crank-Nicholsonovi metodi
določi vrednosti u(j*dx,t+dt), j=1,...,J. Bistveno upošteva, da je
g(t)=h(t)=0."

CrankNicholsonStep[u_List] := CrankNicholsonMatrix . u

CrankNicholsonAt1::usage = "CrankNicholsonAt1 po Crank-Nicholsonovi metodi
določi notranje vrednosti u(j*dx,t=1), j=1,...,J. Vrne seznam s temi
vrednostmi. Če je metoda dobra, mora biti seznam blizu u1."

CrankNicholsonAt1 := Module[ {},
   CrankNicholsonMatrix = Inverse[ Table[
         Switch[ i-j, -1, -1, 0, 4, 1, -1, _, 0],
         {i, J}, {j, J} ] ] .
      Table[ Switch[ i-j, -1, 1, 1, 1, _, 0],
         {i, J}, {j, J} ];
   Nest[ CrankNicholsonStep, u0, Round[1/(dx)^2]]
]

Print["Roman Maurer, uporabna matematika, 29-Aug-1996"]
Print["4. domača naloga iz Numerične analize:"];
Print["Naloga 1. c; rešitev na [0,1]x[0,1] PDE  u_t = u_xx,"];
Print["               u(x,0) = sin(pi*x),"];
Print["               u(0,t) = u(1,t) = 0."];
Print["- z eksplicitno metodo (dx=0.1, lambda=1/4),"];
Print["- z implicitno metodo (dx=0.1, lambda=1/4),"];
Print["- s Crank-Nicholsonovo metodo (dx=0.1, lambda=1)."];
Print[" "];
Print["metoda                   čas      rel. napaka [%]"];
Print["eksplicitna:        ", Timing[100*((ExplicitAt1/u1)[[1]]-1) ]];
Print["implicitna:         ", Timing[100*((ImplicitAt1/u1)[[1]]-1) ]];
Print["Crank-Nicholson:    ", Timing[100*((CrankNicholsonAt1/u1)[[1]]-1) ]];

(* konec datoteke NA4.M *)