Algorithmen des wissenschaftlichen Rechnens II
Übungsblatt 5, 16.01.2007
Normen und Interpolation mit dünnen Gittern bezüglich der Energienorm

>    restart:
with(plots):
with(plottools):
with(LinearAlgebra):

Warning, the name changecoords has been redefined

Warning, the assigned name arrow now has a global binding

Ein paar generelle Optionen zum Plotten:

>    poptions := thickness=2, scaling=constrained;
poptions2d := thickness=2, scaling=constrained, axes=boxed;

poptions := thickness = 2, scaling = constrained

poptions2d := thickness = 2, scaling = constrained, axes = boxed

>   

Hilfsmittel vom letzten Blatt

Unsere Basisfunktionen phi[l,i](x) , die uns unsere Hüte in Abhängigkeit von Level l  und Index i  zurückliefern:

>    phi_li := (l,i) -> ((x) -> piecewise(x <= (i-1)/2^l, 0,
                                     x <= i/2^l,     2^l*x+1-i,
                                     x <= (i+1)/2^l, i+1-2^l*x,
                                     0));

phi_li := proc (l, i) options operator, arrow; proc (x) options operator, arrow; piecewise(x <= (i-1)/(2^l),0,x <= i/(2^l),2^l*x+1-i,x <= (i+1)/(2^l),i+1-2^l*x,0) end proc end proc

Für ein gegebenes Gitter können wir uns die Basisfunktionen gleich mit Hilfe einer Funktion ausgeben lassen:

>    plot_basisfunktionen := proc(gitter::list, u::Vector, hash::table)
  # gitter ... ein Gitter
  # optional:    u    ... Koeffizientenvektor
  #              hash ... eine passende Hashtabelle
  if nargs = 1 then
    return plot([seq(phi_li(gitterpunkt[1],gitterpunkt[2])(x),
                     gitterpunkt in gitter)], x=0..1, poptions);
  else
    return plot([seq(u[hash[gitterpunkt]]*phi_li(gitterpunkt[1],gitterpunkt[2])(x),
                     gitterpunkt in gitter)], x=0..1, poptions);
  end if;
end proc:

Und 2D:

>    plot_basisfunktionen_2D := proc(gitter::list, u::Vector, hash::table)
  # gitter ... ein Gitter
  # optional:    u    ... Koeffizientenvektor
  #              hash ... eine passende Hashtabelle
  if nargs = 1 then
    return plot3d([seq(phi_li(gitterpunkt[1],gitterpunkt[2])(x1)
                       *phi_li(gitterpunkt[3],gitterpunkt[4])(x2),
                       gitterpunkt in gitter)], x1=0..1, x2=0..1, poptions2d);
  else
    return plot3d([seq(u[hash[gitterpunkt]]
                       *phi_li(gitterpunkt[1],gitterpunkt[2])(x1)
                       *phi_li(gitterpunkt[3],gitterpunkt[4])(x2),
                       gitterpunkt in gitter)], x1=0..1, x2=0..1, poptions2d);
  end if;
end proc:

Hashtabelle, die uns passend zu einem Eintrag [l, i]  einen Index im Vektor u  zurückliefert.

>    hash_tabelle := proc(gitter::list)
  local i, hash;
  # Hashtabelle erstellen
  hash := table();
  # Hashtabelle füllen:
  for i from 1 to nops(gitter) do
    hash[gitter[i]] := i;
  end do:
  return hash;
end proc:

Prozedur, die uns zu gegebenem Level n  das Gitter zurückliefert. Dies geht direkt nach der Definition der Basis:

>    gitter_hierarch_1D := proc(n::integer)
  return [seq(
              seq([l,2*i-1],
                  i=1..2^(l-1)),
              l=1..n)];
end proc:

Hierarchisieren

>    hierarchisieren_1D := proc(u::evaln(Vector), n::integer, hash::table)
  # u    ... Koeffizientenvektor
  # n    ... max. Level
  # hash ... Hashtabelle [l,i] -> Index-in-u
  local hierarch_rec;

  hierarch_rec := proc(fl, fr, l, i)

    if l < n then
      hierarch_rec(fl, u[hash[[l,i]]], l+1, 2*i-1);
      hierarch_rec(u[hash[[l,i]]], fr, l+1, 2*i+1);
    end if;

    u[hash[[l,i]]] := u[hash[[l,i]]] - 1/2*(fl+fr);

  end proc:

  hierarch_rec(0, 0, 1, 1);

end proc:

2D-Fall zu gegebenem n, Hierarchisierung (leicht modifiziert):

>    hierarchisieren_2D := proc(u::evaln(Vector), n::integer, hash::table)
  # u    ... Koeffizientenvektor
  # n    ... max. Level
  # hash ... Hashtabelle [l1,i1,l2,i2] -> Index-in-u
  local hierarch_rec, gitter_1D, gitterpunkt;

  hierarch_rec := proc(praefix::list, suffix::list, fl, fr, l, i, lmax)
    local gitterpunkt;
    gitterpunkt := [op(praefix),l,i,op(suffix)];

    # hier machen wir eine leichte Modifikation, damit wir beliebige
    # (auch adaptive) Gitter hierarchisieren können. Statt:
    #if l < lmax then
    # als Abfrage prüfen wir nun, ob die Kinder-Basisfunktionen existieren:
    if assigned(hash[[op(praefix),l+1,2*i-1,op(suffix)]]) then
      hierarch_rec(praefix, suffix, fl, u[hash[gitterpunkt]], l+1, 2*i-1, lmax);
    end if;
    if assigned(hash[[op(praefix),l+1,2*i+1,op(suffix)]]) then
      hierarch_rec(praefix, suffix, u[hash[gitterpunkt]], fr, l+1, 2*i+1, lmax);
    end if;

    u[hash[gitterpunkt]] := u[hash[gitterpunkt]] - 1/2*(fl+fr);

  end proc:

  # Erstellen Gitter 1D
  gitter_1D := gitter_hierarch_1D(n);

  # Hierarchisieren in x1-Richtung, d.h. die x2-Komponente wird fixiert:
  for gitterpunkt in gitter_1D do
    hierarch_rec([], gitterpunkt, 0, 0, 1, 1, n-gitterpunkt[1]+1);
  end do;
  # Hierarchisieren in x2-Richtung, d.h. die x1-Komponente wird fixiert:
  for gitterpunkt in gitter_1D do
    hierarch_rec(gitterpunkt, [], 0, 0, 1, 1, n-gitterpunkt[1]+1);
  end do;

end proc:

>   

Gittererzeugung

Kartes. Produkt zweier Listen.

>    kartprod := proc(l1::list, l2::list)
  return map(n -> op(map(m -> [op(m), op(n)], l1)), l2);
end proc:

Volles Gitter:

>    gitter_dD := proc(dim::integer, level::integer)
  if dim = 1 then
    return gitter_1D(level);
  else
    return kartprod(gitter_1D(level),
                    gitter_dD(dim-1, level));
  end if;
end proc:

>   

Teilraumschema

Jeder der Teilräume W[l]  wird festgelegt durch den Vektor l , den wir im Folgenden als Liste relisieren.
Für einen bestimmten Teilraum können die Basisfunktionen ähnlich wie oben erzeugt werden:

>    gitter_W := proc(l::list)
  if nops(l) = 1 then   # nur 1D übrig
    return [seq([l[1],2*i-1], i=1..2^(l[1]-1))];
  else
    return kartprod([seq([l[1],2*i-1], i=1..2^(l[1]-1))],
                    gitter_W(l[2..nops(l)]));
  end if;
end proc:

Die alte Dünngitter-Prozedur:
Alle Teilräume ablaufen. Bei den dünnen Gittern bezüglich der
L[2] -Norm sind dies alle W[l]  mit sum(l[i],i = 1 .. d) <= n+d-1 .
Wir beginnen mit l=[1,...,1]

>    gitter_hierarch_dD := proc(dim::integer, n::integer)
  local gitter_rec;

  gitter_rec := proc(l::list, d::integer)
    local gitter_d_rec, gitter_l_rec, d_rec;

    # Rekursion im Level in aktueller Dimension d?
    if (l[d] < n) and (add(l[i],i=1..dim) < n+dim-1) then
      gitter_l_rec := gitter_rec([seq(l[i],i=1..d-1),
                                  l[d]+1,
                                  seq(l[i],i=d+1..dim)], d);
    else
      gitter_l_rec := [];
    end if;

    # Rekursion in der Dimension?
    gitter_d_rec := [];
    if (d < dim) and (add(l[i],i=1..dim) < n+dim-1) then
      for d_rec from d+1 to dim do
        gitter_d_rec := [op(gitter_d_rec),
                         op(gitter_rec([seq(l[i],i=1..d_rec-1),
                                        l[d_rec]+1,
                                        seq(l[i],i=d_rec+1..dim)], d_rec))];
      end do;
    end if;

    # Gitterpunkte aktueller Teilraum + Ergebnis aus Rekursionen
    return [op(gitter_W(l)), op(gitter_l_rec), op(gitter_d_rec)];
  end proc:

  return gitter_rec([1$dim], 1);
end proc:

>   

>   

Und, noch neu, eine Prozedur zum plotten von 1-, 2- bzw, 3-D-Gittern:

>    plot_grid := proc(l::list)
  description "plot eines Gitters, gegeben als Liste von Gitterpunkten in 1D, 2D or 3D";
  local d;
  if nops(l) > 0 then d := nops(l[1])/2; else d := 0; end if;
  if d = 1 then
    return display([seq(point([l[i][2]/2^l[i][1],0], symbol=circle), i=1..nops(l))],
                   view=[0..1,-1..1], axes=normal);
  elif d = 2 then
    return display([seq(point([seq(l[i][2*j]/2^l[i][2*j-1], j=1..d)], symbol=circle), i=1..nops(l))],
                   view=[seq(0..1, j=1..d)], poptions, axes=boxed);
  elif d = 3 then
    return display([seq(point([seq(l[i][2*j]/2^l[i][2*j-1], j=1..d)], symbol=circle), i=1..nops(l))],
                   view=[seq(0..1, j=1..d)], poptions2d, orientation=[-80, 80]);
  else
    error("Can't plot grid with dimensionality", d);
  end if;
end proc:

Für den Teilraum W[3,5]  sieht das Gitter dann so aus:

>    plot_grid(gitter_W([3,5]));

[Maple Plot]

>   

>   

>   

>   

>   

Normen und volles Gitter:

Zuerst wollen wir unsere drei Normen aufstellen, wir beginnen mit der Maximumsnorm.
Eigentlich würden wir etwas wie folgt machen (da bei uns der Rand 0 ist, genügt das):

>    Norm_max := proc(dim::integer, f::procedure)
  local x, sols;
  # 1 Gleichung pro dim: diff(f,x[i])=0 und diese lösen
  sols := solve({seq(diff(g(seq(x[j],j=1..dim)),x[i])=0,i=1..dim)},
                {seq(x[i],i=1..dim)});
  # Maximum von |f| bestimmen, ausgewertet an den berechneten Stellen
  return eval(max(seq(subs(l,abs(g(seq(x[j],j=1..dim)))),l in sols)));
end proc:

Praktisch funktioniert das für stückweise lineare Funktionen natürlich nicht.
Wir müssen f an verschiedenen Stellen auswerten (ja, da würde der Fluch der
Dimensionalität wieder zuschlagen, für unsere Zwecke wird es aber reichen...),
und hoffen, dass wir alle Maxima damit erwischen:

>    Norm_max := proc(dim::integer, f::procedure, level::integer)
  local punkte, j;
  # Maximum von |f| bestimmen, ausgewertet an den Gitterpunkten des
  # vollen Gitters zum nächsthöheren Level
  punkte := [[]]; # alle Auswertungspunkte
  for j from 1 to dim do
    punkte := kartprod(punkte, [seq([i*2^(-level-1)],i=1..2^(level+1)-1)]);
  end do;
map(x->f(op(x)),punkte);
  return eval(max(op(map(x->f(op(x)),punkte))));
end proc:

Nun die L[2] -Norm:

>    Norm_L2_alt := proc(dim::integer, f::procedure)
  local x;
  return sqrt(evalf(Int(f(seq(x[i],i=1..dim))^2,[seq(x[i]=0..1,i=1..dim)])));
end proc:

So wäre das simpel mit Maples Integrationsmethode (siehe ?int[numeric] ). Allerdings kommt Maple beim höherdimensionalen Integrieren mit den eigenen piecewise -Funktionen nicht mehr klar. Daher verwenden wir für 2 <= d  die Monte-Carlo Integration von Maple:

>    Norm_L2 := proc(dim::integer, f::procedure)
  local x;
  return sqrt(Int_SimpsonSumme(dim, n, x->f(x)^2));
end proc:
Norm_L2 := proc(dim::integer, f::procedure)
  local x;
if dim = 1 then
  return sqrt(evalf(Int(f(seq(x[i],i=1..dim))^2,
                        [seq(x[i]=0..1,i=1..dim)])));
else
  return sqrt(evalf(Int(f(seq(x[i],i=1..dim))^2,
                    [seq(x[i]=0..1,i=1..dim)],
                     method=_MonteCarlo, epsilon = 0.5e-2)));
end if;
end proc:

Und die Energienorm:

>    Norm_E_alt := proc(dim::integer, f::procedure)
  local x;
  return sqrt(evalf(Int(add(diff(f(seq(x[i],i=1..dim)),x[j])^2, j=1..dim),
                        [seq(x[i]=0..1,i=1..dim)])));
end proc:

Und wieder mit Monte-Carlo:

>    Norm_E := proc(dim::integer, f::procedure)
  local x;
  if dim = 1 then
    return sqrt(evalf(Int(add(diff(f(seq(x[i],i=1..dim)),x[j])^2, j=1..dim),
                          [seq(x[i]=0..1,i=1..dim)])));
  else
    return sqrt(evalf(Int(add(diff(f(seq(x[i],i=1..dim)),x[j])^2, j=1..dim),
                      [seq(x[i]=0..1,i=1..dim)],
                       method=_MonteCarlo, epsilon = 0.5e-2)));
  end if;
end proc:

>   

>   

Testen wir das gleich für unsere Basisfunktionen im 1D und vergleichen das Ergebnis mit den Ergebnissen aus der Vorlesung:

>    seq(simplify(Norm_max(1,phi_li(l,1),l))
    = 1,l=1..4);

>    seq(simplify(Norm_L2(1,phi_li(l,1)))
    = sqrt(evalf(2*2^(-l)/3)),l=1..4);

>    seq(simplify(Norm_E(1,phi_li(l,1)))
    = sqrt(evalf(2/2^(-l))),l=1..4);

1 = 1, 1 = 1, 1 = 1, 1 = 1

.5773502692 = .5773502692, .4082482905 = .4082482905, .2886751346 = .2886751346, .2041241452 = .2041241452
.5773502692 = .5773502692, .4082482905 = .4082482905, .2886751346 = .2886751346, .2041241452 = .2041241452

2. = 2.000000000, 2.828427125 = 2.828427125, 4. = 4.000000000, 5.656854249 = 5.656854249

Um nicht nur die L[2] -Norm-optimalen Gitter, sondern Gitter erzeugen zu können,
die eine bestimmte Auswahl an Teilräumen umfassen, verändern wir obige Prozedur
und machen das Abbruchkriterium variabel. Der neue Parameter
abbr_krit  bekommt
als Parameter
l , das Array der Level, übergeben.

>    gitter_hierarch_dD := proc(dim::integer, n::integer, abbr_krit::procedure)
  local gitter_rec;

  gitter_rec := proc(l::list, d::integer)
    local gitter_d_rec, gitter_l_rec, d_rec;

    if not abbr_krit(l) then
      # Rekursion im Level in aktueller Dimension d?
      if (l[d] < n) then
        gitter_l_rec := gitter_rec([seq(l[i],i=1..d-1),
                                    l[d]+1,
                                    seq(l[i],i=d+1..dim)], d);
      else
        gitter_l_rec := [];
      end if;

      # Rekursion in der Dimension?
      gitter_d_rec := [];
      if (d < dim) then
        for d_rec from d+1 to dim do
          gitter_d_rec := [op(gitter_d_rec),
                           op(gitter_rec([seq(l[i],i=1..d_rec-1),
                                          l[d_rec]+1,
                                          seq(l[i],i=d_rec+1..dim)], d_rec))];
        end do;
      end if;

      # Gitterpunkte aktueller Teilraum + Ergebnis aus Rekursionen
      return [op(gitter_W(l)), op(gitter_l_rec), op(gitter_d_rec)];
    else
      # Aktueller Teilraum soll nicht mehr dazugenommen werden
      return [];
    end if;
  end proc:

  return gitter_rec([1$dim], 1);

end proc:

>   

Volle Gitter erhalten wir damit mit:

>    gitter_V_inf := proc(dim::integer, n::integer)
  return gitter_hierarch_dD(dim, n, l->max(seq(l[i],i=1..dim))>n);
end proc:

Das volle Gitter 2D, Level 3 sieht wie zu erwarten so aus:

>    plot_grid(gitter_V_inf(2,3));

[Maple Plot]

>   

>   

Interpolation volles Gitter, Fehlernorm

Nun wollen wir, zunächst für den 1D-Fall, für gegebenes n  den Interpolanten bestimmen.
Dazu gehen wir vor, wie auf dem letzten Blatt: Funktionsauswertungen am Gitterpunkt und Hierarchisieren.

>    interpol_V_inf_1D := proc(n::integer, f::procedure)
  local gitter, hash, u, i;
  # Gitter und Hashtabelle erstellen
  gitter := gitter_V_inf(1,n);
  hash := hash_tabelle(gitter);
  # Koeffizientenvektor u
  u := Vector(nops(gitter)):
  # mit Funktionswerten füllen
  for i from 1 to nops(gitter) do
    u[hash[gitter[i]]] := f( gitter[i][2] * 2^(-gitter[i][1]));
  end do:
  # und hierarchisieren
  hierarchisieren_1D(u, n, hash):
  # Interpolanten explizit aufstellen und zurückgeben
  return unapply(simplify(add(u[hash[gitterpunkt]] * phi_li(gitterpunkt[1], gitterpunkt[2])(x),
                          gitterpunkt in gitter)),x);
end proc:

Nun stellen wir den Interpolanten für unsere 1D-Funktion f  auf:

>    f := x->4*x*(1-x);
f_N := interpol_V_inf_1D(3,f);
plot([f(x),f_N(x)],x=0..1, poptions);

f := proc (x) options operator, arrow; 4*x*(1-x) end proc

f_N := proc (x) options operator, arrow; piecewise(x <= 0,0,x <= 1/8,7/2*x,x <= 1/4,5/2*x+1/8,x <= 3/8,3/2*x+3/8,x <= 1/2,1/2*x+3/4,x <= 5/8,5/4-1/2*x,x <= 3/4,15/8-3/2*x,x <= 7/8,21/8-5/2*x,x <= 1,7/2...
f_N := proc (x) options operator, arrow; piecewise(x <= 0,0,x <= 1/8,7/2*x,x <= 1/4,5/2*x+1/8,x <= 3/8,3/2*x+3/8,x <= 1/2,1/2*x+3/4,x <= 5/8,5/4-1/2*x,x <= 3/4,15/8-3/2*x,x <= 7/8,21/8-5/2*x,x <= 1,7/2...

[Maple Plot]

Nun haben wir alles, um für n=1,...,8 im 1D-Fall den Interpolationsfehler in allen drei Normen zu messen und zu plotten:

>    err_max := []:
err_L2 := []:
err_E := []:
for n from 1 to 8 do
  print(n);
  f_N := interpol_V_inf_1D(n,f);
  f(x)-f_N(x);
  f_err := unapply(f(x)-f_N(x),x);
  print(plot([f(x), f_N(x), f_err(x)],x=0..1,color=[red,green,blue]));
  err_max := [op(err_max), [n, Norm_max(1, f_err,n)]];
  err_L2 := [op(err_L2), [n, Norm_L2(1, f_err)]];
  err_E := [op(err_E), [n, Norm_E(1, f_err)]];
end do:

1

[Maple Plot]

2

[Maple Plot]

3

[Maple Plot]

4

[Maple Plot]

5

[Maple Plot]

6

[Maple Plot]

7

[Maple Plot]

8

[Maple Plot]

Man sieht schon, dass die Berechnung eine Weile dauert. Quadratisch ist doch nicht ganz so gut...
Die Fehler sind dann:

>    err_max;
err_L2;
err_E;

[[1, 1/4], [2, 1/16], [3, 1/64], [4, 1/256], [5, 1/1024], [6, 1/4096], [7, 1/16384], [8, 1/65536]]

[[1, .1825741858], [2, .4564354646e-1], [3, .1141088661e-1], [4, .2852721654e-2], [5, .7131804134e-3], [6, .1782951034e-3], [7, .4457377583e-4], [8, .1114344396e-4]]
[[1, .1825741858], [2, .4564354646e-1], [3, .1141088661e-1], [4, .2852721654e-2], [5, .7131804134e-3], [6, .1782951034e-3], [7, .4457377583e-4], [8, .1114344396e-4]]
[[1, .1825741858], [2, .4564354646e-1], [3, .1141088661e-1], [4, .2852721654e-2], [5, .7131804134e-3], [6, .1782951034e-3], [7, .4457377583e-4], [8, .1114344396e-4]]

[[1, 1.154700538], [2, .5773502692], [3, .2886751346], [4, .1443375673], [5, .7216878365e-1], [6, .3608439182e-1], [7, .1804219591e-1], [8, .9021097956e-2]]
[[1, 1.154700538], [2, .5773502692], [3, .2886751346], [4, .1443375673], [5, .7216878365e-1], [6, .3608439182e-1], [7, .1804219591e-1], [8, .9021097956e-2]]

Um welchen Faktor ändert sich der Fehler bei einer Halbierung der Maschenweite?

>    seq(err_max[i+1][2]/err_max[i][2], i=1..nops(err_max)-1);
seq(err_L2[i+1][2]/err_L2[i][2], i=1..nops(err_L2)-1);
seq(err_E[i+1][2]/err_E[i][2], i=1..nops(err_E)-1);  

1/4, 1/4, 1/4, 1/4, 1/4, 1/4, 1/4

.2500000001, .2499999999, .2500000001, .2500000000, .2500000001, .2499999999, .2500000001
.2500000001, .2499999999, .2500000001, .2500000000, .2500000001, .2499999999, .2500000001

.5000000002, .5000000000, .5000000000, .5000000000, .4999999999, .5000000000, .5000000001
.5000000002, .5000000000, .5000000000, .5000000000, .4999999999, .5000000000, .5000000001

Bei der Maximumsnorm sieht man sehr schön, dass der Fehler O(h[n]^2)  ist, ebenso bei der L[2] -Norm,
bei der Energienorm jedoch
O(h[n]) . Im Plot sieht das dann so aus:

>    logplot([err_max, err_L2, err_E], poptions);

[Maple Plot]

>   

>   

Und wie sieht das für den 2D-Fall aus?

>    interpol_V_inf_2D := proc(n::integer, f::procedure)
  local gitter, hash, u, i;
  # Gitter und Hashtabelle erstellen
  gitter := gitter_V_inf(2,n);
  hash := hash_tabelle(gitter);
  # Koeffizientenvektor u
  u := Vector(nops(gitter)):
  # mit Funktionswerten füllen
  for i from 1 to nops(gitter) do
    u[hash[gitter[i]]] := f( gitter[i][2] * 2^(-gitter[i][1]),
                             gitter[i][4] * 2^(-gitter[i][3]));
  end do:
  # und hierarchisieren
  hierarchisieren_2D(u, n, hash):
  # Interpolanten explizit aufstellen und zurückgeben
  return unapply(simplify(add(u[hash[gitterpunkt]] * phi_li(gitterpunkt[1], gitterpunkt[2])(x1)
                                                   * phi_li(gitterpunkt[3], gitterpunkt[4])(x2),
                              gitterpunkt in gitter)),x1,x2);
end proc:

Messen wir wieder die Interpolationsfehler, dieses Mal gleich für weniger n:

>    f := (x1,x2) -> 16*x1*(x1-1)*x2*(x2-1);

>    err_max := []:
err_L2 := []:
err_E := []:
for n from 1 to 6 do
  print(n);
  f_N := interpol_V_inf_2D(n,f):
  f_err := unapply(f(x1,x2)-f_N(x1,x2),x1,x2):
  #print(plot3d(f_err(x1,x2),x1=0..1,x2=0..1, poptions2d));
  err_max := [op(err_max), [n, Norm_max(2, f_err, n)]];
  err_L2 := [op(err_L2), [n, Norm_L2(2, f_err)]];
  err_E := [op(err_E), [n, Norm_E(2, f_err)]];
end do:

f := proc (x1, x2) options operator, arrow; 16*x1*(x1-1)*x2*(x2-1) end proc

1

2

3

4

5

6

Die Fehler sind dann:

>    err_max;
err_L2;
err_E;

[[1, 5/16], [2, 29/256], [3, 125/4096], [4, 509/65536], [5, 2045/1048576], [6, 8189/16777216]]

[[1, .2200218753], [2, .5983629170e-1], [3, .1525282425e-1], [4, .3827264460e-2], [5, .9586309662e-3], [6, .2398117505e-3]]
[[1, .2200218753], [2, .5983629170e-1], [3, .1525282425e-1], [4, .3827264460e-2], [5, .9586309662e-3], [6, .2398117505e-3]]

[[1, 1.299591926], [2, .6146050169], [3, .3003693594], [4, .1492194929], [5, .7455250007e-1], [6, .3720646731e-1]]
[[1, 1.299591926], [2, .6146050169], [3, .3003693594], [4, .1492194929], [5, .7455250007e-1], [6, .3720646731e-1]]

Um welchen Faktor ändert sich der Fehler bei einer Halbierung der Maschenweite?

>    seq(err_max[i+1][2]/err_max[i][2], i=1..nops(err_max)-1);
seq(err_L2[i+1][2]/err_L2[i][2], i=1..nops(err_L2)-1);
seq(err_E[i+1][2]/err_E[i][2], i=1..nops(err_E)-1);  

29/80, 125/464, 509/2000, 2045/8144, 8189/32720

.2719561026, .2549092502, .2509216914, .2504741902, .2501606551

.4729215414, .4887193419, .4967866669, .4996163612, .4990639787

Wieder sieht man, sehr schön, dass der Fehler bei der Maximums- und der L[2] -norm O(h[n]^2)  ist,
bei der Energienorm jedoch
O(h[n]) . Im Plot sieht das dann so aus:

>    logplot([err_max, err_L2, err_E], poptions);

[Maple Plot]

>   

>   

>   

>   

>   

>   

Kosten/Nutzen

Betrachten wir nun die Kosten pro Teilraum.

>    c_1D := l -> 2^(l-1);
c_2D := (l1,l2) -> 2^(l1+l2-2);

c_1D := proc (l) options operator, arrow; 2^(l-1) end proc

c_2D := proc (l1, l2) options operator, arrow; 2^(l1+l2-2) end proc

>   

Das war noch unabhängig von der Norm. Der Nutzen muss jedoch für alle drei Normen einzeln aufgestellt werden:

>    b_1D_max :=       l -> 2^(-2*l-1);

>    b_1D_L2  :=       l -> 2^(-2*l)/3;

>    b_1D_E   :=       l -> 2^(-l-1);

>    b_2D_max := (l1,l2) -> 2^(-2*l1-2*l2-2);

>    b_2D_L2  := (l1,l2) -> 2^(-2*l1-2*l2)/9;

>    b_2D_E   := (l1,l2) -> sqrt(1/48 * (2^(-2*l1-4*l2)+2^(-4*l1-2*l2)));

b_1D_max := proc (l) options operator, arrow; 2^(-2*l-1) end proc

b_1D_L2 := proc (l) options operator, arrow; 1/3*2^(-2*l) end proc

b_1D_E := proc (l) options operator, arrow; 2^(-l-1) end proc

b_2D_max := proc (l1, l2) options operator, arrow; 2^(-2*l1-2*l2-2) end proc

b_2D_L2 := proc (l1, l2) options operator, arrow; 1/9*2^(-2*l1-2*l2) end proc

b_2D_E := proc (l1, l2) options operator, arrow; sqrt(1/48*2^(-2*l1-4*l2)+1/48*2^(-4*l1-2*l2)) end proc

>   

Im Zweidimensionalen wollen wir uns näher anschauen, für welche Teilräume W[l]  das Kosten/Nutzen-Verhältnis konstant ist.
Für die Maximumsnorm gilt:

>    cbr_max := simplify(c_2D(l1,l2)/b_2D_max(l1,l2));

cbr_max := 8^(l1+l2)

Im Schaubild der "kontinuierlichen" (l1,l2) ergibt dies folgende Äquipotentiallinien:

>    contourplot(cbr_max, l1=0..10,l2=0..10, contours=10, poptions);

[Maple Plot]

Entsprechend verhält es sich für die L[2] -Norm:

>    cbr_L2 := simplify(c_2D(l1,l2)/b_2D_L2(l1,l2));

cbr_L2 := 9/4*8^(l1+l2)

>    contourplot(cbr_L2, l1=0..10,l2=0..10, contours=10, poptions);

[Maple Plot]

Neu ist das Verhalten für die Energienorm:

>    cbr_E := simplify(c_2D(l1,l2)/b_2D_E(l1,l2));

cbr_E := 3*2^(l1+l2)/(3*4^(-l1)*16^(-l2)+3*16^(-l1)*4^(-l2))^(1/2)

Hier sieht das Verhältnis grundlegend anders aus. Wir erwarten daher auch andersgeartete Äquipotientiallinien:

>    contourplot(cbr_E, l1=0..10,l2=0..10, grid=[100,100], contours=10, poptions);

[Maple Plot]

Noch etwas deutlicher wird der Unterschied nach der Umformung aus der Übung (und alles außer n auf die Linke Seite gebracht):

Für die Max.- bzw. L[2] -Norm:

>    contourplot(l1+l2-1, l1=0..10,l2=0..10, grid=[100,100], contours=10, poptions);

Für die Energienorm:

>    contourplot(((l1+l2)-log[2](2^(-2*l1)+2^(-2*l2))-5)/4, l1=0..10,l2=0..10, grid=[100,100], contours=10, poptions);

[Maple Plot]

[Maple Plot]

Anmerkung: die rechte Seite wird plausibel, wenn man den Term für verschiedene l[2] , l[1] = 1 , auswertet:

>    ls := 4*(l1+l2)-log[2](2^(-2*l1)+2^(-2*l2)):
evalf(seq(subs({l1=j,l2=1},ls),j=1..10));

9.000000000, 13.67807190, 17.91253716, 21.97763219, 25.99437545, 29.99859180, 33.99964782, 37.99991195, 41.99997799, 45.99999450
9.000000000, 13.67807190, 17.91253716, 21.97763219, 25.99437545, 29.99859180, 33.99964782, 37.99991195, 41.99997799, 45.99999450

>   

>   

>   

>   

Dünne Gitter

Das reguläre dünne Gitter nach der Max- bzw. L[2] -Norm:

>    gitter_V_1 := proc(dim::integer, n::integer)
  return gitter_hierarch_dD(dim, n, l->add(l[i],i=1..dim)>n+dim-1);
end proc:

Das Gitter Level 6 sieht dann so aus:

>    plot_grid(gitter_V_1(2,6));

[Maple Plot]

>   

Und für die Energienorm, hier nur für 2D. Hier müssen wir ein wenig tricksen und abrunden, da unsere Schranke sonst zu ungenau ist:

>    gitter_V_E := proc(n::integer)
  return gitter_hierarch_dD(2, n, l->floor(4*(l[1]+l[2])-log[2](2^(-2*l[1])+2^(-2*l[2])))>4*n+5);
end proc:

Das Gitter Level 6 sieht dann so aus:

>    plot_grid(gitter_V_E(6));

[Maple Plot]

>   

>   

>   

Approximationsordnung, 2D

Nun wollen wir für die dünnen Gitter noch die Approximationsordnungen für gegebenes n  ermitteln.
Dazu gehen wir vor wie beim vollen Gitter oben. Copy'n'paste macht hier das Leben einfacher.
Für den 1D-Fall sieht es aus wie oben, daher gleich den 2D-Fall.

Hier sehen die Gitter für die unterschiedlichen Normen allerdings nicht mehr gleich aus:

>    interpol_V_1 := proc(n::integer, f::procedure)
  local gitter, hash, u, i;
  # Gitter und Hashtabelle erstellen
  gitter := gitter_V_1(2,n);
  hash := hash_tabelle(gitter);
  # Koeffizientenvektor u
  u := Vector(nops(gitter)):
  # mit Funktionswerten füllen
  for i from 1 to nops(gitter) do
    u[hash[gitter[i]]] := f( gitter[i][2] * 2^(-gitter[i][1]),
                             gitter[i][4] * 2^(-gitter[i][3]));
  end do:
  # und hierarchisieren
  hierarchisieren_2D(u, n, hash):
  # Interpolanten explizit aufstellen und zurückgeben
  return unapply(simplify(add(u[hash[gitterpunkt]] * phi_li(gitterpunkt[1], gitterpunkt[2])(x1)
                                                   * phi_li(gitterpunkt[3], gitterpunkt[4])(x2),
                              gitterpunkt in gitter)),x1,x2);
end proc:

>    interpol_V_E := proc(n::integer, f::procedure)
  local gitter, hash, u, i;
  # Gitter und Hashtabelle erstellen
  gitter := gitter_V_E(n);
  hash := hash_tabelle(gitter);
  # Koeffizientenvektor u
  u := Vector(nops(gitter)):
  # mit Funktionswerten füllen
  for i from 1 to nops(gitter) do
    u[hash[gitter[i]]] := f( gitter[i][2] * 2^(-gitter[i][1]),
                             gitter[i][4] * 2^(-gitter[i][3]));
  end do:
  # und hierarchisieren
  hierarchisieren_2D(u, n, hash):
  # Interpolanten explizit aufstellen und zurückgeben
  return unapply(simplify(add(u[hash[gitterpunkt]] * phi_li(gitterpunkt[1], gitterpunkt[2])(x1)
                                                   * phi_li(gitterpunkt[3], gitterpunkt[4])(x2),
                              gitterpunkt in gitter)),x1,x2);
end proc:

Und gleich testweise einen Interpolanten:

>    f := (x1,x2) -> 16*x1*(x1-1)*x2*(x2-1):
f_N := interpol_V_E(3,f):
plot3d([f_N(x1,x2)],x1=0..1,x2=0..1, poptions2d);

[Maple Plot]

>   

Messen wir wieder die Interpolationsfehler für einige n :

>    err_max := []:
err_L2 := []:
err_E := []:
for n from 1 to 6 do
  print(n);
  f_N := interpol_V_1(n,f):
  f_err := unapply(f(x1,x2)-f_N(x1,x2),x1,x2):
  #print(plot3d(f_err(x1,x2),x1=0..1,x2=0..1, poptions2d));
  err_max := [op(err_max), [n, Norm_max(2, f_err, n)]];
  err_L2 := [op(err_L2), [n, Norm_L2(2, f_err)]];

  f_N := interpol_V_E(n,f):
  f_err := unapply(f(x1,x2)-f_N(x1,x2),x1,x2):
  #print(plot3d(f_err(x1,x2),x1=0..1,x2=0..1, poptions2d));
  err_E := [op(err_E), [n, Norm_E(2, f_err)]];
end do:

1

2

3

4

5

6

Die Fehler sind dann:

>    err_max;
err_L2;
err_E;

[[1, 5/16], [2, 33/256], [3, 161/4096], [4, 817/65536], [5, 3825/1048576], [6, 17841/16777216]]

[[1, .2197808050], [2, .7618018596e-1], [3, .2428778048e-1], [4, .7382911644e-2], [5, .2174130709e-2], [6, .6251736673e-3]]
[[1, .2197808050], [2, .7618018596e-1], [3, .2428778048e-1], [4, .7382911644e-2], [5, .2174130709e-2], [6, .6251736673e-3]]

[[1, 1.301329560], [2, .6610283665], [3, .4166823234], [4, .2170640648], [5, .1098948399], [6, .5507686593e-1]]
[[1, 1.301329560], [2, .6610283665], [3, .4166823234], [4, .2170640648], [5, .1098948399], [6, .5507686593e-1]]

Um welchen Faktor ändert sich der Fehler bei einer Halbierung der Maschenweite?

>    evalf(seq(err_max[i+1][2]/err_max[i][2], i=1..nops(err_max)-1));
seq(err_L2[i+1][2]/err_L2[i][2], i=1..nops(err_L2)-1);
seq(err_E[i+1][2]/err_E[i][2], i=1..nops(err_E)-1);  

.4125000000, .3049242424, .3171583851, .2926101591, .2915196078

.3466189232, .3188201784, .3039763823, .2944814748, .2875510956

.5079638447, .6303546784, .5209341808, .5062783653, .5011779077

Jetzt kann man sehen, dass die Konvergenzordnung des Fehlers bei der Maximums- und der L[2] -Norm schlechter als O(h[n]^2)  ist,
während es bei der Energienorm noch so gut aussieht wie beim vollen Gitter, nämlich
O(h[n]) .
Im Plot sieht das dann so aus:

>    logplot([err_max, err_L2, err_E], poptions);

[Maple Plot]

>   

>   

>   

>   

>   

>   

>   

>   

>   

>   

>   

>   

>   

>   

>