Algorithmen des wissenschaftlichen Rechnens II
Übungsblatt 4, 19.12.2006
Hierarchische Basis, dünne Gitter und Interpolation

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

Warning, the name changecoords has been redefined

Ein paar generelle Optionen zum Plotten:

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

poptions := thickness = 2

poptions2d := thickness = 2, axes = boxed

>   

Basisfunktionen 1D

Knotenbasis

Zuerst basteln wir uns 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

Dann bekommen wir die Knotenbasis  zur Maschenweite h[n] := 2^(-n)  mit n = 3 :

>    Phi := [seq(phi_li(3,i)(x), i=1..2^3-1)];

Phi := [PIECEWISE([0, x <= 0],[8*x, x <= 1/8],[2-8*x, x <= 1/4],[0, otherwise]), PIECEWISE([0, x <= 1/8],[8*x-1, x <= 1/4],[3-8*x, x <= 3/8],[0, otherwise]), PIECEWISE([0, x <= 1/4],[8*x-2, x <= 3/8],[...
Phi := [PIECEWISE([0, x <= 0],[8*x, x <= 1/8],[2-8*x, x <= 1/4],[0, otherwise]), PIECEWISE([0, x <= 1/8],[8*x-1, x <= 1/4],[3-8*x, x <= 3/8],[0, otherwise]), PIECEWISE([0, x <= 1/4],[8*x-2, x <= 3/8],[...
Phi := [PIECEWISE([0, x <= 0],[8*x, x <= 1/8],[2-8*x, x <= 1/4],[0, otherwise]), PIECEWISE([0, x <= 1/8],[8*x-1, x <= 1/4],[3-8*x, x <= 3/8],[0, otherwise]), PIECEWISE([0, x <= 1/4],[8*x-2, x <= 3/8],[...

Und die können wir auch gleich zeichnen:

>    plot(Phi, x=0..1, poptions);

[Maple Plot]

>   

Nun erstellen wir ein Gitter, d.h eine Liste mit Gitterpunkten, spezifiziert durch eine Liste aus Level und Index

>    gitter_1D := proc(level::integer)
  return [seq([level,i], i=1..2^level-1)];
end proc:

Für Level 3 gibt das:

>    gitter := gitter_1D(3);

gitter := [[3, 1], [3, 2], [3, 3], [3, 4], [3, 5], [3, 6], [3, 7]]

Für ein gegebenes Gitter können wir uns die Basisfunktionen gleich mit Hilfe einer Prozedur 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 gleich ausprobieren:

>    plot_basisfunktionen(gitter_1D(4));

[Maple Plot]

Nun erstellen wir den Koeffizientenvektor u :

>    u := Vector(nops(gitter));

u := Vector(%id = 23528524)

Diesen wollen wir mit Koeffizienten passend zu den Gitterpunkten in der Liste der Gitterpunkte füllen. Wir hätten im Hinblick auf die Verwendung im Höherdimensionalen gerne wahlfreien Zugriff und erstellen uns deshalb eine Hahstabelle, die uns passend zu einem Eintrag [l, i]  einen Index im Vektor u  zurückliefert. (Für den 1D-Fall lohnt dies nicht wirklich.)

>    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:

Für das Gitter zum Level 3 erhalten wir:

>    hash := hash_tabelle(gitter);
print(hash);

hash := hash

TABLE([[3, 4] = 4, [3, 7] = 7, [3, 2] = 2, [3, 5] = 5, [3, 6] = 6, [3, 3] = 3, [3, 1] = 1])

>   

Nun wenden wir uns unserer Funktion f(x)  zu

>    f := x -> 4*x*(1-x);
plot(f(x), x=0..1, poptions);

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

[Maple Plot]

... und füllen den Vektor u mit den Koeffizienten der Basisfunktionen für den Interpolanten f[N] , was im Fall der Knotenbasis einfach nur die Funktionswerte von f  am jeweilgen Gitterpunkt sind:

>    for i from 1 to nops(gitter) do
  u[hash[gitter[i]]] := f( gitter[i][2] * 2^(-gitter[i][1]));
end do;
u;

u[1] := 7/16

u[2] := 3/4

u[3] := 15/16

u[4] := 1

u[5] := 15/16

u[6] := 3/4

u[7] := 7/16

Vector(%id = 23528524)

Und zeichnen das gleich:

>    plot_basisfunktionen(gitter, u, hash);

[Maple Plot]

Nun können wir unseren Interpolanten f[N]  explizit aufstellen.
Hier zeigt sich bereits der Vorteil der Hashtabelle: wir benötigen keinerlei Indexrechnungen:

>    f_N := add(u[hash[gitterpunkt]]
           * phi_li(gitterpunkt[1], gitterpunkt[2])(x),
           gitterpunkt in gitter);

f_N := 7/16*PIECEWISE([0, x <= 0],[8*x, x <= 1/8],[2-8*x, x <= 1/4],[0, otherwise])+3/4*PIECEWISE([0, x <= 1/8],[8*x-1, x <= 1/4],[3-8*x, x <= 3/8],[0, otherwise])+15/16*PIECEWISE([0, x <= 1/4],[8*x-2,...
f_N := 7/16*PIECEWISE([0, x <= 0],[8*x, x <= 1/8],[2-8*x, x <= 1/4],[0, otherwise])+3/4*PIECEWISE([0, x <= 1/8],[8*x-1, x <= 1/4],[3-8*x, x <= 3/8],[0, otherwise])+15/16*PIECEWISE([0, x <= 1/4],[8*x-2,...
f_N := 7/16*PIECEWISE([0, x <= 0],[8*x, x <= 1/8],[2-8*x, x <= 1/4],[0, otherwise])+3/4*PIECEWISE([0, x <= 1/8],[8*x-1, x <= 1/4],[3-8*x, x <= 3/8],[0, otherwise])+15/16*PIECEWISE([0, x <= 1/4],[8*x-2,...

Dass diese Funkion wirklich stückweise linear ist mit Maschenweite h[3] = 1/8  sehen wir direkt, wenn wir Maple vereinfachen lassen:

>    simplify(f_N);

PIECEWISE([0, x <= 0],[7/2*x, x <= 1/8],[1/8+5/2*x, x <= 1/4],[3/8+3/2*x, x <= 3/8],[3/4+1/2*x, x <= 1/2],[5/4-1/2*x, x <= 5/8],[15/8-3/2*x, x <= 3/4],[21/8-5/2*x, x <= 7/8],[7/2-7/2*x, x <= 1],[0, 1 <...

Oder natürlich auch, wenn wir Funktion und Interpolanten plotten:

>    plot([f_N, f(x)], x=0..1, poptions);

[Maple Plot]

>   

>   

>   

Hierarchische Basis

Nun wenden wir uns der Darstellung des Interpolanten in der hierarchischen Basis zu. Einen großen Teil der bisherigen Funktionen ( phi_li , plot_basisfunktionen , hash_tabelle ) können wir direkt wiederverwenden.
Was wir neu benötigen, ist eine 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:

Für Level 3 lassen wir uns gleich zur Kontrolle die Basisfunktionen darstellen:

>    gitter := gitter_hierarch_1D(3);
plot_basisfunktionen(gitter);

gitter := [[1, 1], [2, 1], [2, 3], [3, 1], [3, 3], [3, 5], [3, 7]]

[Maple Plot]

Die passende Hashtabelle dazu:

>    hash := hash_tabelle(gitter);
print(hash);

hash := hash

TABLE([[2, 3] = 3, [1, 1] = 1, [3, 7] = 7, [3, 5] = 6, [3, 3] = 5, [3, 1] = 4, [2, 1] = 2])

Wieder füllen wir den Koeffizientenvektor u  mit den Funktionswerten, die nun u.A. in anderer Reihenfolge im Vektor stehen:

>    u := Vector(nops(gitter)):
for i from 1 to nops(gitter) do
  u[hash[gitter[i]]] := f( gitter[i][2] * 2^(-gitter[i][1]));
end do:
u;

Vector(%id = 23194904)

Nun müssen wir von der nicht-hierarchischen in die hierarchische Darstellung umrechnen, also Hierarchisieren.
Man könnte jetzt jeden hierarchischen Überschuss direkt nach der Formel
u[l,i] := f((a+b)/2)-(f(a)+f(b))/2  berechnen - dies bedeutet allerdings, dass manche Funktionswerte mehrfach berechnet werden. Sind Funktionsauswertungen teuer, so ist dies sicherlich unerwünscht, und wir hätten uns auch nicht die Mühe machen müssen, den Vektor u  bereits zu füllen. Es geht auch ohne erneute Funktionsauswertungen -- wir haben bereits alle benötigten Funktionswerte in u  stehen.

Das geht am Besten rekursiv. Dabei ist fl  der Funktionswert am linken, fr  am rechten Rand des Trägers der Basisfunktion:

>    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:

Nun testen wir, was dies für unser Gitter und unsere Basisfunktion gibt

>    u;
hierarchisieren_1D(u, 3, hash):
u;

Vector(%id = 23194904)

Vector(%id = 23194904)

Schauen wir uns unsere gewichteten Basisfunktionen an:

>    plot_basisfunktionen(gitter, u, hash);

[Maple Plot]

Das sieht schon sehr gut aus.
Nun stellen wir wieder unseren Interpolanten
f[N]  auf.

>    f_N := add(u[hash[gitterpunkt]]
           * phi_li(gitterpunkt[1], gitterpunkt[2])(x),
           gitterpunkt in gitter);

f_N := PIECEWISE([0, x <= 0],[2*x, x <= 1/2],[2-2*x, x <= 1],[0, otherwise])+1/4*PIECEWISE([0, x <= 0],[4*x, x <= 1/4],[2-4*x, x <= 1/2],[0, otherwise])+1/4*PIECEWISE([0, x <= 1/2],[4*x-2, x <= 3/4],[4...
f_N := PIECEWISE([0, x <= 0],[2*x, x <= 1/2],[2-2*x, x <= 1],[0, otherwise])+1/4*PIECEWISE([0, x <= 0],[4*x, x <= 1/4],[2-4*x, x <= 1/2],[0, otherwise])+1/4*PIECEWISE([0, x <= 1/2],[4*x-2, x <= 3/4],[4...
f_N := PIECEWISE([0, x <= 0],[2*x, x <= 1/2],[2-2*x, x <= 1],[0, otherwise])+1/4*PIECEWISE([0, x <= 0],[4*x, x <= 1/4],[2-4*x, x <= 1/2],[0, otherwise])+1/4*PIECEWISE([0, x <= 1/2],[4*x-2, x <= 3/4],[4...

Und plotten ihn zusammen mit f :

>    plot([f_N, f(x)], x=0..1, poptions);

[Maple Plot]

Und wie gewünscht erhalten wir denselben Interpolanten.

>   

>   

>   

Basisfunktionen 2D

Unsere neue Funktion f :

>    f := (x1,x2) -> 16*x1*(1-x1)*x2*(1-x2);
plot3d(f(x1,x2), x1=0..1, x2=0..1, poptions2d);

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

[Maple Plot]

Analog zum 1D-Fall erstellen wir eine Prozedur, die uns für ein gegebenes Gitter die Basisfunktionen malt:

>    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:

>   

Gittererzeugung

Die Gittererzeugung wollen wir uns gleich d-dimensional betrachten, der wesentliche Schritt erfolgt von einer auf zwei Dimensionen.
Für ein reguläres, d-dimensionales Gitter ist dies relativ einfach. Wir haben die 1D-Version zum Level n gegeben, benötigen die in jeder Dimension und bilden das d-dimensionale kartesische Produkt daraus.

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

>    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:

Gleich der Test, ob wir auch etwas Sinnvolles tun für d = 2  und n = 3 :

>    gitter_dD(2,3);

[[3, 1, 3, 1], [3, 2, 3, 1], [3, 3, 3, 1], [3, 4, 3, 1], [3, 5, 3, 1], [3, 6, 3, 1], [3, 7, 3, 1], [3, 1, 3, 2], [3, 2, 3, 2], [3, 3, 3, 2], [3, 4, 3, 2], [3, 5, 3, 2], [3, 6, 3, 2], [3, 7, 3, 2], [3, ...
[[3, 1, 3, 1], [3, 2, 3, 1], [3, 3, 3, 1], [3, 4, 3, 1], [3, 5, 3, 1], [3, 6, 3, 1], [3, 7, 3, 1], [3, 1, 3, 2], [3, 2, 3, 2], [3, 3, 3, 2], [3, 4, 3, 2], [3, 5, 3, 2], [3, 6, 3, 2], [3, 7, 3, 2], [3, ...
[[3, 1, 3, 1], [3, 2, 3, 1], [3, 3, 3, 1], [3, 4, 3, 1], [3, 5, 3, 1], [3, 6, 3, 1], [3, 7, 3, 1], [3, 1, 3, 2], [3, 2, 3, 2], [3, 3, 3, 2], [3, 4, 3, 2], [3, 5, 3, 2], [3, 6, 3, 2], [3, 7, 3, 2], [3, ...
[[3, 1, 3, 1], [3, 2, 3, 1], [3, 3, 3, 1], [3, 4, 3, 1], [3, 5, 3, 1], [3, 6, 3, 1], [3, 7, 3, 1], [3, 1, 3, 2], [3, 2, 3, 2], [3, 3, 3, 2], [3, 4, 3, 2], [3, 5, 3, 2], [3, 6, 3, 2], [3, 7, 3, 2], [3, ...
[[3, 1, 3, 1], [3, 2, 3, 1], [3, 3, 3, 1], [3, 4, 3, 1], [3, 5, 3, 1], [3, 6, 3, 1], [3, 7, 3, 1], [3, 1, 3, 2], [3, 2, 3, 2], [3, 3, 3, 2], [3, 4, 3, 2], [3, 5, 3, 2], [3, 6, 3, 2], [3, 7, 3, 2], [3, ...
[[3, 1, 3, 1], [3, 2, 3, 1], [3, 3, 3, 1], [3, 4, 3, 1], [3, 5, 3, 1], [3, 6, 3, 1], [3, 7, 3, 1], [3, 1, 3, 2], [3, 2, 3, 2], [3, 3, 3, 2], [3, 4, 3, 2], [3, 5, 3, 2], [3, 6, 3, 2], [3, 7, 3, 2], [3, ...
[[3, 1, 3, 1], [3, 2, 3, 1], [3, 3, 3, 1], [3, 4, 3, 1], [3, 5, 3, 1], [3, 6, 3, 1], [3, 7, 3, 1], [3, 1, 3, 2], [3, 2, 3, 2], [3, 3, 3, 2], [3, 4, 3, 2], [3, 5, 3, 2], [3, 6, 3, 2], [3, 7, 3, 2], [3, ...

>    plot_basisfunktionen_2D(gitter_dD(2,3));

[Maple Plot]

Und noch der Test für d=3 und n=2:

>    gitter_dD(3,2);

[[2, 1, 2, 1, 2, 1], [2, 2, 2, 1, 2, 1], [2, 3, 2, 1, 2, 1], [2, 1, 2, 2, 2, 1], [2, 2, 2, 2, 2, 1], [2, 3, 2, 2, 2, 1], [2, 1, 2, 3, 2, 1], [2, 2, 2, 3, 2, 1], [2, 3, 2, 3, 2, 1], [2, 1, 2, 1, 2, 2], ...
[[2, 1, 2, 1, 2, 1], [2, 2, 2, 1, 2, 1], [2, 3, 2, 1, 2, 1], [2, 1, 2, 2, 2, 1], [2, 2, 2, 2, 2, 1], [2, 3, 2, 2, 2, 1], [2, 1, 2, 3, 2, 1], [2, 2, 2, 3, 2, 1], [2, 3, 2, 3, 2, 1], [2, 1, 2, 1, 2, 2], ...
[[2, 1, 2, 1, 2, 1], [2, 2, 2, 1, 2, 1], [2, 3, 2, 1, 2, 1], [2, 1, 2, 2, 2, 1], [2, 2, 2, 2, 2, 1], [2, 3, 2, 2, 2, 1], [2, 1, 2, 3, 2, 1], [2, 2, 2, 3, 2, 1], [2, 3, 2, 3, 2, 1], [2, 1, 2, 1, 2, 2], ...
[[2, 1, 2, 1, 2, 1], [2, 2, 2, 1, 2, 1], [2, 3, 2, 1, 2, 1], [2, 1, 2, 2, 2, 1], [2, 2, 2, 2, 2, 1], [2, 3, 2, 2, 2, 1], [2, 1, 2, 3, 2, 1], [2, 2, 2, 3, 2, 1], [2, 3, 2, 3, 2, 1], [2, 1, 2, 1, 2, 2], ...
[[2, 1, 2, 1, 2, 1], [2, 2, 2, 1, 2, 1], [2, 3, 2, 1, 2, 1], [2, 1, 2, 2, 2, 1], [2, 2, 2, 2, 2, 1], [2, 3, 2, 2, 2, 1], [2, 1, 2, 3, 2, 1], [2, 2, 2, 3, 2, 1], [2, 3, 2, 3, 2, 1], [2, 1, 2, 1, 2, 2], ...
[[2, 1, 2, 1, 2, 1], [2, 2, 2, 1, 2, 1], [2, 3, 2, 1, 2, 1], [2, 1, 2, 2, 2, 1], [2, 2, 2, 2, 2, 1], [2, 3, 2, 2, 2, 1], [2, 1, 2, 3, 2, 1], [2, 2, 2, 3, 2, 1], [2, 3, 2, 3, 2, 1], [2, 1, 2, 1, 2, 2], ...

Wir erwarten (2^2-1)^3 = 27  Gitterpunkte:

>    nops(%);

27

>   

Zur Erstellung eines dünnen Gitters muss man sich schon mehr einfallen lassen.
Es gibt mehrere Möglichkeiten, eine Liste der Gitterpunkte eines dünnen Gitters zu erstellen. Wir wollen uns zwei davon näher betrachten

>   

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:

Testen wir es für ein paar W[l] :

>    gitter_W([1,1]);

>    gitter_W([1,3]);

>    gitter_W([3,2]);

[[1, 1, 1, 1]]

[[1, 1, 3, 1], [1, 1, 3, 3], [1, 1, 3, 5], [1, 1, 3, 7]]

[[3, 1, 2, 1], [3, 3, 2, 1], [3, 5, 2, 1], [3, 7, 2, 1], [3, 1, 2, 3], [3, 3, 2, 3], [3, 5, 2, 3], [3, 7, 2, 3]]
[[3, 1, 2, 1], [3, 3, 2, 1], [3, 5, 2, 1], [3, 7, 2, 1], [3, 1, 2, 3], [3, 3, 2, 3], [3, 5, 2, 3], [3, 7, 2, 3]]

Und graphisch:

>    plot_basisfunktionen_2D(gitter_W([1,1]));

>    plot_basisfunktionen_2D(gitter_W([1,3]));

>    plot_basisfunktionen_2D(gitter_W([3,2]));

[Maple Plot]

[Maple Plot]

[Maple Plot]

Nun müssen wir nur noch 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 gleich wieder ein paar Tests:

>    gitter_hierarch_dD(1,1);

[[1, 1]]

>    gitter_hierarch_1D(3);
gitter_hierarch_dD(1,3);

[[1, 1], [2, 1], [2, 3], [3, 1], [3, 3], [3, 5], [3, 7]]

[[1, 1], [2, 1], [2, 3], [3, 1], [3, 3], [3, 5], [3, 7]]

>    gitter_hierarch_dD(3,1);

[[1, 1, 1, 1, 1, 1]]

>    gitter_hierarch_dD(2,2);

[[1, 1, 1, 1], [2, 1, 1, 1], [2, 3, 1, 1], [1, 1, 2, 1], [1, 1, 2, 3]]

>   

>   

>   

Rekursiver Gitteraufbau

Alternativ können wir die rekursive Struktur dünner Gitter ausnutzen, also die Rekursionsformel a[n,d] := a[n,d-1]+2*a[n-1,d]  aus der Vorlesung.
Hierzu benötigen wir im Wesentlichen eine Prozedur, die uns, "zentriert" an einem bestimmten Gitterpunkt, ein dünnes Gitter zum Level
n  zurückliefert.
Abgebrochen wird, wenn die Rekursion in der Dimension die letzte verbliebene Dimension erreicht. In diesem Fall können wir in der letzten Komponente ein dünnes Gitter zum verbliebenen Level mit
gitter_hierarch_1D  erzeugen, und die vorherigen, bereits betrachteten Komponenten mittels kartprod  davor stellen.
Also beispielsweise in
x[2] -Richtung:

>    gitter_hierarch_1D(2);

Mit der Gitterkomponente [2, 1]  in x[1] -Richtung

>    kartprod([[2,1]], gitter_hierarch_1D(2));

[[1, 1], [2, 1], [2, 3]]

[[2, 1, 1, 1], [2, 1, 2, 1], [2, 1, 2, 3]]

>   

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

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

    if d=1 then
      # return 1D-Gitter
      return gitter_hierarch_1D(l);
    else
      # rekursiver Aufruf
      # Rekursion in der Dimension "a_{n,d-1}"
      gitter_d_rec := kartprod([current],
                               gitter_rec([1,1], d-1, l));
      # Rekursion im Level "a_{n-1,d}"
      if l > 1 then
        gitter_l_rec :=
      [op(gitter_rec([current[1]+1,2*current[2]-1], d, l-1)),
       op(gitter_rec([current[1]+1,2*current[2]+1], d, l-1))];
      else
        gitter_l_rec := [];
      end if;
      # return Gitter
      return [op(gitter_d_rec), op(gitter_l_rec)];
    end if;
  end proc;

  # Start implizit mit dem Gitterpunkt [1,1, 1,1, ..., 1,1] und damit
  # mit noch (n+d-1)-d verfügbaren Level für den Abstieg
  return gitter_rec([1,1], dim, n);
end proc:

Und gleich wieder ein paar Tests:

>    gitter_hierarch_rec_dD(1,1);

[[1, 1]]

>    gitter_hierarch_1D(3);
gitter_hierarch_rec_dD(1,3);

[[1, 1], [2, 1], [2, 3], [3, 1], [3, 3], [3, 5], [3, 7]]

[[1, 1], [2, 1], [2, 3], [3, 1], [3, 3], [3, 5], [3, 7]]

>    gitter_hierarch_rec_dD(3,1);

[[1, 1, 1, 1, 1, 1]]

>    gitter_hierarch_rec_dD(2,2);

[[1, 1, 1, 1], [1, 1, 2, 1], [1, 1, 2, 3], [2, 1, 1, 1], [2, 3, 1, 1]]

>   

>   

Zur Sicherheit prüfen wir, ob beide Varianten für höheres Level und höhere Dimension dasselbe erzeugen.
Dazu lassen wir uns beide Listen erzeugen, wandeln sie in Sets um, damit die Reihenfolge keine Rolle spielt und lassen sie vergleichen:

>    evalb(  convert(gitter_hierarch_dD(5,5), set)
      = convert(gitter_hierarch_rec_dD(5,5), set));

true

>   

>   

>   

>   

>   

>   

Hierarchisierung, 2D

Wir können uns nun der Hierarchisierung zuwenden, und das gleich ganz anschaulich zum V[3]  im 2D:

>    gitter := gitter_hierarch_dD(2,3);

gitter := [[1, 1, 1, 1], [2, 1, 1, 1], [2, 3, 1, 1], [3, 1, 1, 1], [3, 3, 1, 1], [3, 5, 1, 1], [3, 7, 1, 1], [2, 1, 2, 1], [2, 3, 2, 1], [2, 1, 2, 3], [2, 3, 2, 3], [1, 1, 2, 1], [1, 1, 2, 3], [1, 1, 3...
gitter := [[1, 1, 1, 1], [2, 1, 1, 1], [2, 3, 1, 1], [3, 1, 1, 1], [3, 3, 1, 1], [3, 5, 1, 1], [3, 7, 1, 1], [2, 1, 2, 1], [2, 3, 2, 1], [2, 1, 2, 3], [2, 3, 2, 3], [1, 1, 2, 1], [1, 1, 2, 3], [1, 1, 3...
gitter := [[1, 1, 1, 1], [2, 1, 1, 1], [2, 3, 1, 1], [3, 1, 1, 1], [3, 3, 1, 1], [3, 5, 1, 1], [3, 7, 1, 1], [2, 1, 2, 1], [2, 3, 2, 1], [2, 1, 2, 3], [2, 3, 2, 3], [1, 1, 2, 1], [1, 1, 2, 3], [1, 1, 3...

Wir erstellen wieder unsere Hashtabelle für den wahlfreien Indexzugriff

>    hash := hash_tabelle(gitter);
print(hash);

hash := hash

TABLE([[2, 3, 2, 1] = 9, [2, 1, 1, 1] = 2, [1, 1, 1, 1] = 1, [3, 5, 1, 1] = 6, [1, 1, 3, 1] = 14, [3, 3, 1, 1] = 5, [2, 3, 1, 1] = 3, [3, 1, 1, 1] = 4, [2, 1, 2, 3] = 10, [1, 1, 2, 3] = 13, [1, 1, 3, 5...
TABLE([[2, 3, 2, 1] = 9, [2, 1, 1, 1] = 2, [1, 1, 1, 1] = 1, [3, 5, 1, 1] = 6, [1, 1, 3, 1] = 14, [3, 3, 1, 1] = 5, [2, 3, 1, 1] = 3, [3, 1, 1, 1] = 4, [2, 1, 2, 3] = 10, [1, 1, 2, 3] = 13, [1, 1, 3, 5...
TABLE([[2, 3, 2, 1] = 9, [2, 1, 1, 1] = 2, [1, 1, 1, 1] = 1, [3, 5, 1, 1] = 6, [1, 1, 3, 1] = 14, [3, 3, 1, 1] = 5, [2, 3, 1, 1] = 3, [3, 1, 1, 1] = 4, [2, 1, 2, 3] = 10, [1, 1, 2, 3] = 13, [1, 1, 3, 5...
TABLE([[2, 3, 2, 1] = 9, [2, 1, 1, 1] = 2, [1, 1, 1, 1] = 1, [3, 5, 1, 1] = 6, [1, 1, 3, 1] = 14, [3, 3, 1, 1] = 5, [2, 3, 1, 1] = 3, [3, 1, 1, 1] = 4, [2, 1, 2, 3] = 10, [1, 1, 2, 3] = 13, [1, 1, 3, 5...
TABLE([[2, 3, 2, 1] = 9, [2, 1, 1, 1] = 2, [1, 1, 1, 1] = 1, [3, 5, 1, 1] = 6, [1, 1, 3, 1] = 14, [3, 3, 1, 1] = 5, [2, 3, 1, 1] = 3, [3, 1, 1, 1] = 4, [2, 1, 2, 3] = 10, [1, 1, 2, 3] = 13, [1, 1, 3, 5...
TABLE([[2, 3, 2, 1] = 9, [2, 1, 1, 1] = 2, [1, 1, 1, 1] = 1, [3, 5, 1, 1] = 6, [1, 1, 3, 1] = 14, [3, 3, 1, 1] = 5, [2, 3, 1, 1] = 3, [3, 1, 1, 1] = 4, [2, 1, 2, 3] = 10, [1, 1, 2, 3] = 13, [1, 1, 3, 5...

und füllen den Vektor u mit den Funktionswerten von f :

>    u := Vector(nops(gitter));
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;
u;

u := Vector(%id = 21200756)

u[1] := 1

u[2] := 3/4

u[3] := 3/4

u[4] := 7/16

u[5] := 15/16

u[6] := 15/16

u[7] := 7/16

u[8] := 9/16

u[9] := 9/16

u[10] := 9/16

u[11] := 9/16

u[12] := 3/4

u[13] := 3/4

u[14] := 7/16

u[15] := 15/16

u[16] := 15/16

u[17] := 7/16

Vector(%id = 21200756)

>   

Nun müssen wir uns "nur" noch um das 2D-Hierarchisieren kümmern.
Wir müssen für die Berechnung des hierarchischen Überschusses einer Basisfunktion die Beiträge aller gewichteten Basisfunktionen, deren Träger am Gitterpunkt nicht null ist, abziehen. Das können nur hierarchisch höhere Basisfunktionen (kleineres Level) sein.
Da unser Interpolant
f[N]  stückweise bilinear ist genügt es wieder, die Werte an den beiden Nachbar-Gitterpunkten zu betrachten.
Als (hoffentlich) anschauliches Beispiel betrachten wir den Gitterpunkt
[2, 1, 2, 1] , sowie die hierarchisch höheren Basisfunktionen [1, 1, 2, 1] , [2, 1, 1, 1]  und [1, 1, 1, 1] .

>    gitter := [[1,1,1,1],[1,1,2,1],[2,1,1,1],[2,1,2,1]];
u := Vector(nops(gitter)):
hash := hash_tabelle(gitter):
u[hash[[1,1,1,1]]] := f(0.5,0.5);
u[hash[[1,1,2,1]]] := f(0.5,0.25);
u[hash[[2,1,1,1]]] := f(0.25,0.5);
u[hash[[2,1,2,1]]] := f(0.25,0.25);

gitter := [[1, 1, 1, 1], [1, 1, 2, 1], [2, 1, 1, 1], [2, 1, 2, 1]]

u[1] := 1.0000

u[2] := .750000

u[3] := .750000

u[4] := .56250000

>    poptions2d_old := poptions2d:
poptions2d := poptions2d, style=WIREFRAME, orientation=[-173,79], color=[grey,green,blue,red]:
plot_basisfunktionen_2D(gitter, u, hash);
poptions2d := poptions2d_old:

[Maple Plot]

Wenn wir zunächst nur die x[1] -Richtung betrachten, so können wir für die Gitterpunkte [2, 1, 1, 1]  und [2, 1, 2, 1]  die Einträge in u  neu berechnen; für [1, 1, 1, 1]  und [1, 1, 2, 1]  gibt es, bei Randwerten 0, nichts zu tun:

>    u[hash[[2,1,1,1]]] := u[hash[[2,1,1,1]]]
                      - ((0 + u[hash[[1,1,1,1]]]) / 2);

>    u[hash[[2,1,2,1]]] := u[hash[[2,1,2,1]]]
                      - ((0 + u[hash[[1,1,2,1]]]) / 2);

u[3] := .2500000000

u[4] := .1875000000

Dann sieht das ganze so aus:

>    poptions2d_old := poptions2d:
poptions2d := poptions2d, style=WIREFRAME, orientation=[-137,79], color=[grey,green,blue,red]:
plot_basisfunktionen_2D(gitter, u, hash);
poptions2d := poptions2d_old:

[Maple Plot]

In x[1] -Richtung sieht es schon gut aus; in -Richtung stimmt es allerdings noch überhaupt nicht. Hier müssen wir einfach ebenfalls eindimensional Hierarchisieren:

>    u[hash[[1,1,2,1]]] := u[hash[[1,1,2,1]]]
                      - ((0 + u[hash[[1,1,1,1]]]) / 2);

>    u[hash[[2,1,2,1]]] := u[hash[[2,1,2,1]]]
                      - ((0 + u[hash[[2,1,1,1]]]) / 2);

u[2] := .2500000000

u[4] := .625000000e-1

Und wir erhalten:

>    poptions2d_old := poptions2d:
poptions2d := poptions2d, style=WIREFRAME, orientation=[-136,89], color=[grey,green,blue,red]:
plot_basisfunktionen_2D(gitter, u, hash);
poptions2d := poptions2d_old:

[Maple Plot]

Der Interpolant f[N]  ist dann:

>    f_N := add(u[hash[gitterpunkt]]
           *phi_li(gitterpunkt[1],gitterpunkt[2])(x1)
           *phi_li(gitterpunkt[3],gitterpunkt[4])(x2),
           gitterpunkt in gitter):
p1 := plot3d(f_N, x1=0..1, x2=0..1, poptions2d):
p2 := plot3d(f(x1,x2), x1=0..1, x2=0..1, poptions2d, style=WIREFRAME, orientation=[-137,79]):
display([p2,p1]);

[Maple Plot]

>   

Nun wollen wir für den 2D-Fall zu gegebenem n die Hierarchisierung durchführen:

>    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)];

    if l < lmax then
      hierarch_rec(praefix, suffix, fl, u[hash[gitterpunkt]], l+1, 2*i-1, lmax);
      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:

Jetzt alles zusammen:

>    dim := 2: n := 3:
gitter := gitter_hierarch_dD(dim,n):
hash := hash_tabelle(gitter):
u := Vector(nops(gitter)):
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:
hierarchisieren_2D(u,n,hash):
f_N := add(u[hash[gitterpunkt]]
           *phi_li(gitterpunkt[1],gitterpunkt[2])(x1)
           *phi_li(gitterpunkt[3],gitterpunkt[4])(x2),
           gitterpunkt in gitter):
p1 := plot3d(f_N, x1=0..1, x2=0..1, poptions2d):
p2 := plot3d(f(x1,x2), x1=0..1, x2=0..1, poptions2d, style=WIREFRAME, orientation=[-137,79], color=grey):
display([p2,p1]);

[Maple Plot]

Die Überschüsse für n=3 sind:

>    for gitterpunkt in gitter do
  print(gitterpunkt, u[hash[gitterpunkt]]);
end do:

[1, 1, 1, 1], 1

[2, 1, 1, 1], 1/4

[2, 3, 1, 1], 1/4

[3, 1, 1, 1], 1/16

[3, 3, 1, 1], 1/16

[3, 5, 1, 1], 1/16

[3, 7, 1, 1], 1/16

[2, 1, 2, 1], 1/16

[2, 3, 2, 1], 1/16

[2, 1, 2, 3], 1/16

[2, 3, 2, 3], 1/16

[1, 1, 2, 1], 1/4

[1, 1, 2, 3], 1/4

[1, 1, 3, 1], 1/16

[1, 1, 3, 3], 1/16

[1, 1, 3, 5], 1/16

[1, 1, 3, 7], 1/16

>   

>   

>   

>   

>   

Auswertung Interpolant 2D:

Hier können wir ausnützen, dass die Basisfunktionen in jedem der Teilräume W[l]  paarweise disjunkte Träger besitzen und somit für einen Auswertepunkt X  nur genau eine Basisfunktion einen Funktionswert ungleich null besitzen kann. Welche das ist, lässt sich einfach mit Hilfe einer Indexrechnung feststellen:

>    dim := 4;
l := [2,4,3,1];
X := Vector([0.8, 0.48, 0.33, 0.78]);
for d from 1 to dim do
  print(x[d], [l[d], ceil(X[d]*2^(l[d]-1))*2-1]);
end do:

dim := 4

l := [2, 4, 3, 1]

X := Vector(%id = 718680)

x[1], [2, 3]

x[2], [4, 7]

x[3], [3, 3]

x[4], [1, 1]

D.h. wie oben müssen wir nur noch alle Teilräume ablaufen.
Wir beginnen mit l=[1,...,1]

>    eval_gitter_hierarch_dD := proc(dim::integer, n::integer,
                               X::Vector, u::Vector, hash::table)
  local gitter_rec, res;

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

    # 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 := 0;
    end if;

    # Rekursion in der Dimension?
    gitter_d_rec := 0;
    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 := gitter_d_rec +
                        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;

    # Ergebnis aktueller Teilraum:
    gitterpunkt :=
          [seq(op([l[i], ceil(X[i]*2^(l[i]-1))*2-1]), i=1..dim)];
    res := u[hash[gitterpunkt]]*
         mul(phi_li(gitterpunkt[2*i-1],gitterpunkt[2*i]) (X[i]),
             i=1..dim);

    #  aktueller Teilraum + Ergebnis aus Rekursionen
    return res + gitter_l_rec + gitter_d_rec;
  end proc:

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

Und testen für eine Funktion im 2D:

>    dim := 2: n := 3:
gitter := gitter_hierarch_dD(dim,n):
hash := hash_tabelle(gitter):
u := Vector(nops(gitter)):
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:
hierarchisieren_2D(u,n,hash):
f_N := unapply(add(u[hash[gitterpunkt]]
               *phi_li(gitterpunkt[1],gitterpunkt[2])(x1)
               *phi_li(gitterpunkt[3],gitterpunkt[4])(x2),
               gitterpunkt in gitter), [x1,x2]):
for i from 1 to 10 do
  X := Vector([(rand(2^10)()+1)/2^10, (rand(2^10)()+1)/2^10]);
  print(eval_gitter_hierarch_dD(dim, n, X, u, hash)
        = f_N(X[1],X[2]));
end do:

469643/524288 = 469643/524288

79709/131072 = 79709/131072

707233/1048576 = 707233/1048576

247645/524288 = 247645/524288

66885/262144 = 66885/262144

187/16384 = 187/16384

157781/524288 = 157781/524288

223933/524288 = 223933/524288

364741/524288 = 364741/524288

63369/131072 = 63369/131072

>   

>   

>   

>   

>   

>