Zwei Ingenieurbüros rechnen denselben Träger: links steckt sein Ende im kühlen 20-°C-Auflager, rechts heizt der Brandherd auf 300 °C. Beide liefern einen sauberen, professionell aussehenden Bericht ab — ein Farbbild des Temperaturverlaufs, glatt von Blau nach Rot, dazu die Kennzahl, auf die es dem Prüfstatiker ankommt: die Wärmeleistung, die der Träger vom heißen zum kühlen Ende durchreicht. Büro A schreibt 140 W. Büro B schreibt 280 W — doppelt so viel. Beide Bilder sehen gleich vertrauenswürdig aus.
Eines der beiden ist um den Faktor zwei falsch. Vom bloßen Ansehen: nicht zu entscheiden. Und — das ist der eigentliche Schrecken — die beiden Farbbilder des Temperaturverlaufs sind sogar identisch. Der Fehler versteckt sich in einer Zahl, die das schöne Bild gar nicht zeigt.
Abbildung 12.1: Das Duell der zwei Berichte. Beide zeigen dasselbe Temperaturbild des Trägers — dieselbe Farbskala, dieselbe Gerade von 20 °C (blau, links) nach 300 °C (rot, rechts). Nur die Kennzahl darunter unterscheidet sich: 140 W gegen 280 W. Welcher Bericht ist falsch? Am Bild ist es nicht zu sehen — das ist die Lektion des Kapitels. Die Auflösung steht am Kapitelende.
Was man hier sieht: zwei Bilder, die man nicht auseinanderhalten kann, und zwei Zahlen, die sich um den Faktor zwei unterscheiden. Ein FEM-Programm rechnet immer ein hübsches Bild aus — auch wenn es sich verrechnet hat. Der Merksatz dieses Kapitels lautet darum:
Bunt ist kein Beweis.
Dieses Kapitel gibt dir die Werkzeuge, in fünf Minuten zu entscheiden, welches der beiden Büros recht hat — und allgemeiner: jeder Simulation die drei Fragen zu stellen, die ein Profi ihr stellt, bevor er ihr glaubt.
Lernziele
Nach diesem Kapitel kannst du …
… Verifikation („löse ich die Gleichungen richtig?“) und Validierung („löse ich die richtigen Gleichungen?”) auseinanderhalten,
… drei analytische Prüffälle des Trägers einsetzen — die Gerade als Nulltest, die Parabel mit Wärmequelle und die transiente Abkühlkurve — und aus einer Abweichung auf ihre Ursache schließen,
… eine Konvergenzstudie durchführen: den Fehler bei \(h\), \(h/2\), \(h/4\) messen und am Viertelungs-Rhythmus erkennen, dass die Rechnung sauber konvergiert,
… die Energiebilanz als Alltagstest nutzen (rein \(-\) raus \(=\) gespeichert, auf 1 % genau?),
… die fünf häufigsten Anfängerfehler an ihren Symptomen erkennen und mit der Simulations-Checkliste systematisch ausschließen.
WarnungNaheliegende Vermutung
Vermutung:„Je feiner das Netz, desto richtiger — irgendwann stimmt es dann von selbst.”
Warum sie naheliegt: In Kapitel 1 haben wir es ja selbst gesehen: mehr Knoten bilden die Wirklichkeit genauer ab. Und tatsächlich schrumpft der Fehler mit jedem Verfeinern.
Was stattdessen stimmt: Verfeinern heilt nur den einen Fehler, der vom groben Netz kommt — den Diskretisierungsfehler. Eine vergessene Randbedingung, ein Einheitenfehler, ein vertauschtes Vorzeichen oder ein Programmierfehler bleiben auf jedem Netz falsch. Schlimmer noch: Auf einem feinen Netz sieht das falsche Ergebnis glatter und vertrauenswürdiger aus als auf einem groben. Ein feines Netz macht ein falsches Bild nicht richtig, sondern nur schöner. Genau deshalb braucht es Prüffälle mit bekannter Antwort.
12.1 Die drei Fragen — und zwei Wörter, die man nicht verwechseln darf
Bevor ein Profi einer Simulation glaubt, stellt er ihr drei Fragen:
Stimmt das Programm? Rechnet es die Gleichungen, die es zu rechnen vorgibt, überhaupt richtig — ohne Programmier- oder Einheitenfehler?
Reicht das Netz? Ist der Diskretisierungsfehler klein genug für den Zweck?
Passen die Randbedingungen zur Wirklichkeit? Bildet das Modell die echte Situation ab?
Die ersten beiden Fragen prüft man rechnerisch, ganz ohne Labor — das ist die Verifikation, das Thema dieses Kapitels. Die dritte führt irgendwann ins Labor, zu echten Messungen — das ist die Validierung.
TippVerifikation und Validierung — zwei Fragen, nicht eine
Zwei Wörter, die im Alltag oft durcheinandergehen und die man als Simulationsingenieur sauber trennt (Roache 1998):
Verifikation:„Löse ich die Gleichungen richtig?“ — Rechnet das Programm die gewählte Mathematik korrekt? Werkzeuge: analytische Prüffälle, Netzkonvergenz, Bilanzen. Braucht kein Labor, nur bekannte Vergleichslösungen. Das kann dieses Buch vollständig.
Validierung:„Löse ich die richtigen Gleichungen?“ — Bildet das Modell die Wirklichkeit ab? Werkzeug: der Vergleich mit Messungen an echten Bauteilen. Deckt auf, was keine noch so saubere Rechnung findet — z. B. dass wir die Wärmeabgabe an die Luft weggelassen haben, oder dass Stahl bei 300 °C anders leitet als bei 20 °C.
Merkregel: Erst verifizieren, dann validieren. Ein Programm mit einem Rechenfehler gegen eine Messung zu halten, verwirrt nur — man weiß dann nicht, ob das Modell oder der Code danebenliegt. Dieses Kapitel verifiziert; die Validierung bleibt sein ehrlicher Ausblick (siehe Das Kleingedruckte).
12.2 Prüffall 1 — die Gerade: der Nulltest
Der erste Prüffall ist der schärfste, und wir haben ihn längst in der Hand. In Kapitel 6 haben wir bewiesen: Ohne Wärmequelle ist das eingeschwungene Temperaturprofil des homogenen Trägers exakt eine Gerade — von den 20 °C links geradlinig auf die 300 °C rechts. Die Formel dazu kennt jeder Schüler:
Das Besondere: Die linearen Elemente dieses Buches stellen eine Gerade exakt dar (ihre Hütchen interpolieren geradlinig, Kapitel 7). Zwischen der exakten Lösung und dem, was die FEM darstellen kann, klafft hier also keine Lücke — kein Diskretisierungsfehler, auf keinem Netz. Die FEM muss die Gerade an jedem Knoten und an jedem Punkt dazwischen punktgenau treffen. Das macht die Gerade zum Nulltest:
TippDer Nulltest — das schärfste Werkzeug
Ein Nulltest ist ein Prüffall, dessen exakte Lösung die FEM ohne jeden Diskretisierungsfehler treffen muss. Läuft er durch: schön. Zeigt er auch nur ein hundertstel Grad Abweichung, ist das garantiert ein Programm- oder Randbedingungsfehler — es gibt keine gutmütige Ausrede „das Netz war zu grob”. Deshalb ist der Nulltest das erste, was man rechnet, wenn ein neues Programm läuft.
Wir lassen die FEM aus Kapitel 7 den quellfreien Träger auf drei Netzen rechnen und vergleichen jeden Knoten mit der Geraden:
Code
def gerade(ort, links_grad, rechts_grad, laenge):"""Die exakte Gerade zwischen den beiden Randtemperaturen."""return links_grad + ort / laenge * (rechts_grad - links_grad)def groesste_abweichung(knotentemperaturen, links_grad, rechts_grad, laenge):"""Die groesste Abweichung der FEM-Knoten von der exakten Geraden (Grad).""" knotenzahl =len(knotentemperaturen) elementlaenge = laenge / (knotenzahl -1) schlimmste =0.0for knoten inrange(knotenzahl): ort = knoten * elementlaenge abweichung = knotentemperaturen[knoten] - gerade( ort, links_grad, rechts_grad, laenge)if abweichung <0.0: abweichung =-abweichungif abweichung > schlimmste: schlimmste = abweichungreturn schlimmste# Die FEM-Knotenwerte des quellfreien Traegers (aus dem Kapitel-Programm) fuer# 5, 9 und 17 Knoten -- hier zur Anschauung schon eingesetzt:fem_5 = [20.0, 90.0, 160.0, 230.0, 300.0]print("Nulltest -- FEM gegen die exakte Gerade:")print(" 5 Knoten: groesste Abweichung =","%.2e Grad"% groesste_abweichung(fem_5, 20.0, 300.0, 1.0))
Nulltest -- FEM gegen die exakte Gerade:
5 Knoten: groesste Abweichung = 0.00e+00 Grad
Interpretation: null (bis auf die letzte Stelle, die der Computer beim Rechnen mit Kommazahlen ohnehin verwackelt). Auf 9 und 17 Knoten sieht es genauso aus — der Nulltest ist auf jedem Netz exakt. Büro A und Büro B aus dem Aufhänger würden diesen Test übrigens beide bestehen; merke dir das für die Auflösung.
12.3 Prüffall 2 — die Parabel: der erste echte Diskretisierungsfehler
Jetzt bauen wir einen Fehler ein, den die FEM nicht mehr wegzaubern kann. Wir stecken den Heizdraht aus Kapitel 6 in den Träger: eine gleichmäßig über die Länge verteilte Wärmequelle, insgesamt \(Q = 160\,\mathrm{W}\) (das ist die Aufgabe aus Übung 7.3). Die Bilanz geht dann nicht mehr mit einer Geraden auf; die exakte Lösung ist eine Parabel über der Geraden:
In der Mitte wölbt sich die Parabel also genau 40 °C über die Gerade — auf 200 °C statt 160 °C. Das ist dieselbe Parabel, die der Kasten in Kapitel 6 angekündigt und die Übung 7.3 ausgerechnet hat. Jetzt ist sie unsere geprüfte Referenz.
Und nun kommt das Verblüffende. Wir lassen die FEM die Parabel auf 5, 9 und 17 Knoten rechnen und messen den Fehler auf zwei Arten: einmal nur an den Knoten, einmal zwischen den Knoten (dort, wo der FEM-Streckenzug die gebogene Kurve geradlinig abkürzt).
Code
import os, syssys.path.insert(0, os.path.join("..", "programme", "kap12"))sys.path.insert(0, os.path.join("..", "programme", "gemeinsam"))from kap12_pruefstation import (loese_stationaer, exakte_parabel, groesster_zwischenfehler)import matplotlib.pyplot as pltknotenliste = [5, 9, 17, 33, 65]fehler = []for kn in knotenliste: fehler.append(groesster_zwischenfehler(loese_stationaer(kn, 160.0), exakte_parabel))fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12.0, 4.4))ax1.loglog(knotenliste, fehler, "o-", color="tab:blue", ms=8)for i inrange(len(knotenliste)): ax1.annotate("%.3f °C"% fehler[i], (knotenliste[i], fehler[i]), textcoords="offset points", xytext=(6, 6), fontsize=8)ax1.set_xlabel("Knotenzahl (feiner →)")ax1.set_ylabel("größter Fehler (°C)")ax1.set_title("Viertelung: jede Netzhalbierung ¼-tel den Fehler")ax1.grid(True, which="both", ls=":", alpha=0.5)fem = loese_stationaer(5, 160.0)xs = [i /200.0for i inrange(201)]fem_kurve = []exakt_kurve = []for x in xs: el =min(int(x /0.25), 3) a = (x - el *0.25) /0.25 fem_kurve.append(fem[el] + a * (fem[el +1] - fem[el])) exakt_kurve.append(exakte_parabel(x))ax2.plot(xs, exakt_kurve, "-", color="black", lw=2.0, label="exakte Parabel")ax2.plot(xs, fem_kurve, "--", color="tab:red", lw=1.8, label="FEM 5 Knoten")ax2.plot([i *0.25for i inrange(5)], fem, "o", color="tab:red", ms=9, label="Knoten (exakt getroffen)")ax2.fill_between(xs, fem_kurve, exakt_kurve, color="tab:red", alpha=0.15)ax2.set_xlabel("Ort (m)")ax2.set_ylabel("Temperatur (°C)")ax2.set_title("Fehlerlandkarte: die Lücke liegt ZWISCHEN den Knoten")ax2.legend(fontsize=8)plt.tight_layout()plt.show()
Abbildung 12.2: Links der Konvergenzplot: der größte Zwischen-Knoten-Fehler der Parabel über der Knotenzahl (beide Achsen logarithmisch — gleiche Faktoren werden zu gleichen Abständen). Von Netz zu Netz viertelt sich der Fehler: 2,50 → 0,625 → 0,156 → 0,039 °C. Rechts die Fehlerlandkarte für 5 Knoten: Die FEM (rot gestrichelt) trifft jeden Knoten (rote Punkte) exakt auf der schwarzen Parabel, schneidet aber zwischen den Knoten die Kurve ab — die rote Fläche ist der ganze Fehler.
Interpretation: Zwei Überraschungen. Erstens: An den Knoten ist die FEM sogar exakt — die roten Punkte liegen punktgenau auf der schwarzen Parabel (der Knotenfehler ist so klein wie beim Nulltest). Das ist eine Besonderheit der 1D-Rechnung. Zweitens, und das ist die eigentliche Lektion: Der Fehler steckt ganz zwischen den Knoten, wo der gerade FEM-Streckenzug die Kurve der Parabel abkürzt. Genau dort greift ein buntes Bild, wenn es zwischen den Knoten interpoliert. Und dieser Zwischen-Knoten-Fehler viertelt sich bei jeder Netzhalbierung: \(2{,}50 \to 0{,}625 \to 0{,}156
\to 0{,}039\,°\mathrm{C}\).
TippDer Viertelungs-Rhythmus (2. Ordnung), ohne Angst vor dem Wort
Halbiert man die Maschenweite \(h\), so viertelt sich der Fehler. Warum viertelt, nicht halbiert? Weil der Fehler eines geraden Stücks über einer Kurve mit dem Quadrat der Stücklänge wächst (\(h^2\)) — halb so lang heißt \((\tfrac{1}{2})^2 = \tfrac{1}{4}\) des Fehlers. Das nennt man zweite Ordnung oder Konvergenzordnung 2. Der Viertelungs-Rhythmus ist der praktische Test dafür: Rechne dieselbe Aufgabe auf \(h\), \(h/2\), \(h/4\) und schau, ob sich der Fehler brav viertelt. Tut er das nicht (viertelt zu wenig oder gar nicht), stimmt etwas nicht — dann ist der Fehler kein Netzfehler, sondern ein echter Bug. Der Log-Plot links macht den Rhythmus sichtbar: Weil beide Achsen gleiche Faktoren zu gleichen Abständen stauchen, liegen die Punkte auf einer Geraden — je steiler, desto höher die Ordnung.
12.4 Dasselbe in der Fläche — und eine Falle beim Testen
In Kapitel 10 und Kapitel 11 haben wir zwei Elementfamilien gebaut, Dreiecke und Rechtecke. Welche ist genauer? Der naheliegende Plan: dieselbe Parabel in 2D rechnen und die Fehler vergleichen. Wir legen den Träger als Fläche an (oben und unten isoliert, also nur in \(x\)-Richtung ein Verlauf) und lassen beide Familien auf immer feineren Netzen los.
Interpretation: null und wieder null — beide Familien treffen die Parabel an jedem Knoten exakt, auf jedem Netz. Heißt das, Dreiecke und Rechtecke sind gleich gut? Nein — es heißt, dass dieser Test sie gar nicht auseinanderhalten kann. Die Parabel krümmt sich nur in einer Richtung; in dieser einen Richtung sind beide Elemente geradlinig und darum an den Knoten exakt. Ein Prüffall testet nur, was er auch beansprucht.
WichtigVorhersage-Punkt
Bevor du weiterliest: Wir bauen jetzt einen Prüffall, der in beide Richtungen krümmt — einen glatten Wärmehügel. Jetzt muss jedes Element die Krümmung mit geraden bzw. sattelförmigen Flächen annähern; ein echter Fehler entsteht. Was glaubst du: Viertelt er sich auch hier (zweite Ordnung)? Und gewinnt das Dreieck oder das Rechteck? Leg dich fest.
Der neue Prüffall ist ein quadratischer Fleck Blech, an allen vier Rändern auf 20 °C gehalten, mit einer Quelle, die in der Mitte am stärksten heizt. Seine exakte Lösung ist eine glatte Kuppel, \(20 + 40\sin(\pi x)\sin(\pi y)\) — ein „Parabel-Hügel”, der in \(x\)und\(y\) krümmt.
Abbildung 12.3: Links der Wärmehügel als Farbbild (exakte Lösung, Peak 60 °C in der Mitte). Rechts die Konvergenz beider Familien: Dreiecke (blau) und Rechtecke (orange) vierteln beide ihren Fehler bei jeder Verfeinerung (zweite Ordnung), aber mit verschiedenem Vorfaktor — hier liegen die Dreiecke etwas günstiger. Ein echter 2D-Test trennt die Familien, die entartete Parabel konnte es nicht.
Interpretation: Jetzt entsteht ein echter Fehler, und beide Familien zähmen ihn im Viertel-Takt (zweite Ordnung — die Log-Geraden sind gleich steil). Sie unterscheiden sich nur im Vorfaktor: Auf demselben Netz sind die Dreiecke hier rund dreimal genauer als die Rechtecke (\(2{,}1\) gegen \(6{,}7\,°\mathrm{C}\) beim gröbsten Netz). Welche Familie „gewinnt”, hängt von der Form der Lösung ab — Kapitel 11 hat schon angedeutet, dass es keinen generellen Sieger gibt. Wichtiger als der Sieger ist die Lehre über das Testen: Ein Prüffall, den zwei verschiedene Verfahren identisch bestehen, vergleicht sie nicht — er ist zu einfach. Man muss den Test so wählen, dass er das beansprucht, was man wissen will.
Die Verfeinerung lässt sich auch als Film ansehen: Netz und Fehler schmelzen gemeinsam im Viertel-Takt.
Code
import os, sys, mathsys.path.insert(0, os.path.join("..", "programme", "kap12"))sys.path.insert(0, os.path.join("..", "programme", "gemeinsam"))from dreiecke import (strukturiertes_dreiecksnetz, assembliere_2d, reduziere_system, setze_ein)from loeser import gauss_eliminationimport numpy as npimport matplotlib.pyplot as pltfrom matplotlib.animation import FuncAnimationfrom IPython.display import HTMLdef exakt(x, y):return20.0+40.0* math.sin(math.pi * x) * math.sin(math.pi * y)def loese(nx): knoten, dreiecke = strukturiertes_dreiecksnetz(nx, nx, 1.0, 1.0) matrix = assembliere_2d(knoten, dreiecke, 50.0, 1.0) seite = [0.0] *len(knoten)for tri in dreiecke: (x1, y1), (x2, y2), (x3, y3) = (knoten[tri[0]], knoten[tri[1]], knoten[tri[2]]) flaeche =abs((x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1)) /2.0for i in tri: xi, yi = knoten[i] q =50.0*2* math.pi**2*40.0* math.sin(math.pi*xi) * math.sin(math.pi*yi) seite[i] += q * flaeche /3.0 randwerte = [None] *len(knoten)for k inrange(len(knoten)): x, y = knoten[k]if x <1e-9or x >1-1e-9or y <1e-9or y >1-1e-9: randwerte[k] = exakt(x, y) km, ks, frei = reduziere_system(matrix, seite, randwerte) voll = setze_ein(gauss_elimination(km, ks), frei, randwerte)return knoten, dreiecke, vollframes = [4, 8, 16, 32]loesungen = [loese(nx) for nx in frames]fig, (axn, axf) = plt.subplots(1, 2, figsize=(11.0, 4.4))def zeichne(idx): axn.clear(); axf.clear() knoten, dreiecke, voll = loesungen[idx] xs = [k[0] for k in knoten]; ys = [k[1] for k in knoten] tris = [list(t) for t in dreiecke] axn.triplot(xs, ys, tris, color="0.4", lw=0.5) axn.set_title("Netz: %d×%d"% (frames[idx], frames[idx]), fontsize=10) axn.set_aspect("equal"); axn.set_xlabel("x (m)") fehler = [abs(voll[k] - exakt(xs[k], ys[k])) for k inrange(len(knoten))] bild = axf.tripcolor(xs, ys, tris, fehler, cmap="Reds", vmin=0, vmax=2.2) axf.set_title("Fehlerbetrag (°C), Peak %.2f"%max(fehler), fontsize=10) axf.set_aspect("equal"); axf.set_xlabel("x (m)")return []ani = FuncAnimation(fig, zeichne, frames=len(frames), interval=900, blit=False)plt.close(fig)HTML(ani.to_jshtml())
Abbildung 12.4: Netzverfeinerung des Wärmehügels (Dreiecke). Von Frame zu Frame verdoppelt sich die Knotenzahl je Richtung; links das Netz, rechts das Fehlerbild (Betrag der Abweichung von der exakten Kuppel, gleiche Farbskala über alle Frames). Der Fehler schmilzt sichtbar im Viertel-Takt — die Skala bleibt fest, damit das Schrumpfen ehrlich sichtbar ist.
12.5 Prüffall 3 — die Abkühlkurve gegen eine exakte Reihe
Die ersten beiden Prüffälle waren stationär. Jetzt bekommt die Zeit ihre Prüfung. Der Träger startet kalt (überall 20 °C), das rechte Ende springt auf die Temperatur des Brandherds, und wir sehen der Wärme beim Hereinwandern zu (Kapitel 9). Gibt es dafür eine exakte Vergleichslösung? Ja — und wir sind ihr längst begegnet.
HinweisWoher kam die Fühler-Kurve aus Kapitel 5?
In Kapitel 5 hatten wir eine glatte Funktion temperatur(ort) als „Messfühler” benutzt — die 30-Minuten-Momentaufnahme des Trägers, an der wir Steigung und Fläche geübt haben. Damals stand dort: „Woher die Formel kommt, entwickeln die Kapitel 6 und 12.” Hier ist die Einlösung. Jene Kurve war die exakte Lösung der Wärmeleitungsgleichung — die stationäre Gerade plus eine abklingende Summe von Sinuswellen (eine Fourier-Reihe):
Man muss diese Formel nicht herleiten — sie ist ein geprüftes Werkzeug, kein Lehrstoff dieses Buches. Zwei Dinge genügen zum Verständnis: Jede Sinuswelle klingt mit der Zeit ab (der Faktor \(e^{-\dots t}\)), und die langsamste Welle (\(n=1\)) bestimmt, wie lange es dauert — ihre Abklingzeit ist \(\tau_1 = \tau/\pi^2 \approx 7320\,\mathrm{s}\), mit der Zeitkonstante \(\tau = L^2\rho c/\lambda \approx 72\,220\,\mathrm{s}\) (rund 20 Stunden) aus Kapitel 9. Genau diese Reihe ist ab jetzt unsere Referenz. Wir prüfen, ob die transiente FEM gegen sie konvergiert.
Wir lassen die implizite FEM aus Kapitel 9 bis \(t = 1800\,\mathrm{s}\) (30 Minuten) rechnen und vergleichen mit der Reihe — einmal bei Netzverfeinerung, einmal bei Zeitschritt-Verfeinerung:
Code
import os, syssys.path.insert(0, os.path.join("..", "programme", "kap12"))from kap12_pruefstation import transienter_fehlerprint("Transient gegen die Reihe (t = 1800 s, rechtes Ende 500 °C):")print()print("Netz verfeinern (dt = 1 s fest):")for kn in (5, 9, 17):print(" %2d Knoten: max Fehler = %.3f °C"% (kn, transienter_fehler(kn, 1.0, 1800.0, 500.0)))print()print("Zeitschritt verfeinern (33 Knoten fest):")vorher =Nonefor dt in (200.0, 100.0, 50.0, 25.0): f = transienter_fehler(33, dt, 1800.0, 500.0) verh =""if vorher isNoneelse" (%.2f-fach kleiner)"% (vorher / f)print(" dt = %6.1f s: max Fehler = %5.3f °C%s"% (dt, f, verh)) vorher = f
Transient gegen die Reihe (t = 1800 s, rechtes Ende 500 °C):
Netz verfeinern (dt = 1 s fest):
5 Knoten: max Fehler = 11.328 °C
9 Knoten: max Fehler = 3.865 °C
17 Knoten: max Fehler = 1.022 °C
Zeitschritt verfeinern (33 Knoten fest):
dt = 200.0 s: max Fehler = 7.414 °C
dt = 100.0 s: max Fehler = 3.751 °C (1.98-fach kleiner)
dt = 50.0 s: max Fehler = 1.917 °C (1.96-fach kleiner)
dt = 25.0 s: max Fehler = 1.001 °C (1.91-fach kleiner)
Interpretation: Die FEM konvergiert sauber gegen die exakte Reihe — die „Fühler-Kurve” aus Kapitel 5 ist also wirklich die richtige Antwort, und die FEM trifft sie im Grenzfall. Aber die beiden Verfeinerungen verhalten sich unterschiedlich: Das Netz strebt die Viertelung an (Raum, zweite Ordnung) — und die Tabelle ist dabei ehrlicher als jede Faustregel: erst 2,9-fach, dann 3,8-fach kleiner. Auf groben Netzen ist der Viertelungs-Rhythmus noch nicht ganz eingerastet, mit jeder Verfeinerung rückt er näher an die 4 heran (genau das meint „asymptotisch”). Der Zeitschritt aber halbiert den Fehler nur (jede Halbierung von \(\Delta t\) halbiert ihn — die zweite Tabelle zeigt es fast aufs Komma).
TippZeit ist gröber als Raum
Der implizite Euler ist nur erster Ordnung: Halbierst du \(\Delta t\), halbierst du den Zeitfehler — nicht mehr. Der Raum ist bei linearen Elementen zweiter Ordnung (viertelt), die Zeit nur erster (halbiert). Wer also die Zeit genauso genau haben will wie den Raum, muss viel kleinere Zeitschritte machen, als das Netz an Feinheit nahelegt. Es gibt bessere Zeitschrittverfahren — das Crank-Nicolson-Verfahren etwa mittelt alten und neuen Zeitpunkt und ist zweiter Ordnung — aber das ist ein anderes Kapitel. Für uns zählt die ehrliche Faustregel: Zeit ist gröber als Raum, also Zeitschritt nicht vergessen, wenn man verfeinert (Übung 12.2).
12.6 Die Energiebilanz — der Wachhund, der immer mitläuft
Prüffälle mit bekannter Lösung sind Gold, aber es gibt sie nicht für jede echte Aufgabe. Ein Test läuft dagegen immer mit, ganz ohne Vergleichslösung: die Energiebilanz. Die Physik erlaubt keine Wärme aus dem Nichts. Über jeden Zeitraum muss darum aufgehen:
Das rechnet man im Programm einfach mit: an jedem Rand den Wärmestrom mal Zeitschritt aufsummieren, am Ende die gespeicherte Wärme aus den Temperaturanstiegen bilden. Geht die Bilanz nicht auf 1 % genau auf, ist etwas faul — eine Wärmequelle, die man vergessen hat, ein Rand, über den heimlich Wärme entweicht, ein Vorzeichen verkehrt herum.
Interpretation: Die Bilanz schließt praktisch perfekt — jedes Joule, das rechts hereinkommt, findet sich als gespeicherte oder links abgeflossene Wärme wieder. Das ist kein Beweis, dass alles stimmt (die Bilanz merkt manche Fehler nicht, siehe Aufhänger-Auflösung), aber ein Ausschlusskriterium: Schließt sie nicht, hast du garantiert einen Fehler. Als Wachhund gehört sie in jedes ernsthafte Programm.
12.7 Die Fehler-Sprechstunde: die fünf Klassiker an ihren Symptomen
Die häufigsten Anfängerfehler hinterlassen typische Symptombilder. Wer sie kennt, diagnostiziert in Sekunden. Hier die Galerie der fünf Klassiker — je mit Symptom und Diagnose-Rezept. Klapp jede Diagnose erst auf, wenn du selbst geraten hast.
Code
import os, syssys.path.insert(0, os.path.join("..", "programme", "kap12"))from kap12_pruefstation import (loese_stationaer, exakte_parabel, sabotiere, simuliere_implizit, exakte_reihe)import matplotlib.pyplot as pltfig, achsen = plt.subplots(1, 5, figsize=(15.0, 3.2))xs = [i /100.0for i inrange(101)]exakt_parabel_kurve = [exakte_parabel(x) for x in xs]def zeichne_profil(ax, temperaturen, titel): kn =len(temperaturen) orte = [i * (1.0/ (kn -1)) for i inrange(kn)] ax.plot(xs, exakt_parabel_kurve, "-", color="black", lw=1.5) ax.plot(orte, temperaturen, "o-", color="tab:red", ms=4) ax.set_title(titel, fontsize=9) ax.set_xlabel("Ort (m)")zeichne_profil(achsen[0], sabotiere(9, 1)[0], "1: Rand vergessen")zeichne_profil(achsen[1], sabotiere(9, 2)[0], "2: Netz zu grob")zeichne_profil(achsen[2], sabotiere(9, 3)[0], "3: Einheiten-Salat")zeichne_profil(achsen[3], sabotiere(9, 4)[0], "4: Vorzeichen")# 5: zu grosser Zeitschritt (transient, verschmierte Front gegen feine Referenz)grob = simuliere_implizit(21, 750.0, 1500.0, 300.0)orte21 = [i * (1.0/20.0) for i inrange(21)]referenz = [exakte_reihe(o, 1500.0, 300.0) for o in orte21]achsen[4].plot(orte21, referenz, "-", color="black", lw=1.5)achsen[4].plot(orte21, grob, "o-", color="tab:red", ms=4)achsen[4].set_title("5: Zeitschritt zu groß", fontsize=9)achsen[4].set_xlabel("Ort (m)")plt.tight_layout()plt.show()
Abbildung 12.5: Die Symptom-Galerie der fünf Klassiker, jeweils als Temperaturprofil des Trägers (schwarz die korrekte Parabel-Lösung, rot der fehlerhafte Lauf). Von links: vergessene Randbedingung (rechtes Ende läuft weg), zu grobes Netz (Knoten stimmen, Form zu kantig), Einheiten-Salat (Höcker zehnfach), Vorzeichenfehler (Kurve hängt durch statt aufzuwölben), zu großer Zeitschritt (Front verschmiert, transient). Jedes Bild ist bunt und plausibel — und trotzdem falsch.
HinweisKlassiker 1 — die vergessene Randbedingung
Symptom: Ein Rand zeigt nicht den Wert, den er zeigen soll — das rechte Ende „läuft weg” statt bei 300 °C zu stehen; oft läuft die halbe Rechnung auf einen unsinnigen Wert. Ursache: Über einen Rand wurde nichts gesagt — und in der FEM heißt „nichts sagen” isoliert (Fluss null), nicht „offen” (Kapitel 8). Diagnose-Rezept: Punkt 1 der Checkliste — jede Randkante durchgehen: Hat sie eine gewollte Bedingung? Der Nulltest deckt das sofort auf, weil die Gerade dann nicht mehr stimmt.
HinweisKlassiker 2 — das zu grobe Netz
Symptom: der heimtückischste. In 1D sehen die Knoten sogar richtig aus (die Mitte zeigt korrekt 200 °C!) — nur zwischen den Knoten ist die Form zu kantig. Ein Farbbild, das nur die Knoten zeigt, verrät nichts. Ursache: Der Diskretisierungsfehler ist noch zu groß für den Zweck. Diagnose-Rezept: die Konvergenzprobe. Rechne dieselbe Aufgabe doppelt so fein. Ändert sich das Ergebnis noch merklich, war das Netz zu grob; viertelt sich die Änderung brav, kannst du abschätzen, wie fein es sein muss. Niemals einem einzigen Netz trauen.
HinweisKlassiker 3 — der Einheiten-Salat
Symptom: die Größenordnung stimmt nicht — der Höcker ist zehnmal zu hoch (560 °C statt 200 °C in der Mitte), oder die Temperaturen sind absurd. Ursache: eine Materialzahl in der falschen Einheit — \(\lambda\) in W/(cm·K) statt W/(m·K), Länge in mm statt m, Grad statt Kelvin, wo das Programm es anders erwartet. Diagnose-Rezept: Checkliste Punkt 5. Rechne eine Kennzahl von Hand grob nach (den Höcker \(q L^2/(8\lambda A)\), den Wärmestrom \(\lambda A\,\Delta T/L\)) und vergleiche die Größenordnung. Einheiten konsequent mitschreiben — der teuerste Fehler der Raumfahrt war ein Einheiten-Salat.
HinweisKlassiker 4 — das falsche Vorzeichen
Symptom: die Kurve krümmt sich falsch herum — sie hängt durch statt sich aufzuwölben (Mitte 120 °C statt 200 °C). Ursache: ein vertauschtes Vorzeichen — der Heizdraht kühlt statt zu heizen, oder ein Wärmestrom fließt rein, wo er raus soll. Diagnose-Rezept: Checkliste Punkt 5 (Vorzeichen des Stroms) und die Energiebilanz. Frag dich: In welche Richtung muss die Wärme fließen? Eine Quelle wölbt die Kurve über die Gerade, eine Senke unter sie.
HinweisKlassiker 5 — der zu große Zeitschritt
Symptom (transient): die Wärmefront ist verschmiert — sie ist weiter ins Kalte gewandert und flacher, als sie sein dürfte. (Beim expliziten Euler wäre das Symptom drastischer: wachsende Zickzack-Ausschläge, die Explosion aus Kapitel 9.) Ursache:\(\Delta t\) zu groß — der implizite Euler bleibt zwar stabil, aber ungenau (erste Ordnung). Diagnose-Rezept: Zeitschritt halbieren und schauen, ob sich die Lösung noch merklich ändert — dieselbe Konvergenzprobe wie beim Netz, nur in der Zeit (Übung 12.2). Zeit ist gröber als Raum.
Diese Erfahrung gießen wir — wie schon die Randbedingungen in Kapitel 8 — in eine feste Prüfroutine. Sie erweitert die dortige Randbedingungs-Checkliste zur Simulations-Checkliste des ganzen Buches; sie kehrt im Finale (Kapitel 15) neben den Reglern wieder.
TippDie Simulations-Checkliste des Buches
Vor jeder Simulation — und vor jedem Glauben an ein buntes Bild:
Randbedingungen (die Checkliste aus Kapitel 8): Hat jeder Rand eine gewollte Bedingung? Mindestens eine Temperatur (Dirichlet) vorgegeben, sonst schwimmt das System? Passt der Typ zur Physik (Thermostat → Dirichlet, Brandherd → Neumann)? Vorzeichen und Einheiten der Randwerte richtig?
Nulltest: Gibt es einen Fall mit exakt bekannter Lösung, den das Programm punktgenau treffen muss (die Gerade, freie Dehnung)? Läuft er durch?
Konvergenz: Auf zwei Netzen (\(h\) und \(h/2\)) gerechnet — viertelt sich die Änderung? Falls transient: auch \(\Delta t\) halbiert?
Energiebilanz: Geht rein \(=\) gespeichert \(+\) raus \((+\) erzeugt\()\) auf 1 % genau auf?
Plausibilität von Hand: Stimmen Größenordnung, Vorzeichen und Richtung einer nachgerechneten Kennzahl (Höcker, Wärmestrom, Zeitkonstante)?
Erst wenn alle fünf Haken sitzen, ist das bunte Bild ein Beweis — vorher ist es nur bunt.
12.8 Die interaktive Einheit: die Prüfstation
Jetzt bedienst du die Prüfstation selbst. Du wählst einen Prüffall, ein Netz und einen Elementtyp; die Station rechnet und zeigt dir Fehlertabelle und Bilanzprotokoll. Und sie hat einen Sabotage-Modus: Stell SABOTAGE auf 1 bis 5, dann baut sie heimlich einen Fehler ein — die vier stationären Klassiker (1 = Rand vergessen, 2 = Netz zu grob, 3 = Einheiten-Salat, 4 = Vorzeichen) und als fünften den gemeinen Aufhänger-Fehler. Du spielst Detektiv: Nutze die Checkliste, sag vorher, welcher Fehler es ist, und klapp dann die Auflösung auf. (Den fünften Galerie-Klassiker, den zu großen Zeitschritt, jagst du in Übung 12.2 — er zeigt sich nur im transienten Fall.)
WichtigVorhersage-Punkt (das Detektivspiel)
Bevor du ausführst: Die Station startet mit SABOTAGE = 5 — dem gemeinsten Fall. Die Temperaturgerade wird völlig normal aussehen, der Nulltest wird bestehen. Trotzdem stimmt etwas nicht. Was für ein Fehler versteckt sich so gut, dass nur eine der fünf Prüfungen der Checkliste ihn fängt — und welche Prüfung ist es? Leg dich fest, dann führe aus.
Abbildung 12.6: Vorgerenderte Fassung der Prüfstation mit dem Startfall SABOTAGE = 5: Die rote FEM-Kurve liegt genau auf der Geraden 20 → 300 (die Parabel-Referenz gehört zum Quellfall und dient hier nur als Maßstab). Am Temperaturbild ist nichts zu sehen — der Faktor-2-Fehler im λ verrät sich erst am Wärmestrom (280 W statt 140 W). Im Browser ersetzt dein eigener Lauf dieses Bild.
12.9 So machen es die Großen
Diese Kultur — nichts glauben, was nicht gegen eine bekannte Lösung geprüft ist — ist keine Schrulle dieses Buches, sondern der Alltag ernsthafter Simulationssoftware. Das quelloffene Programm OpenGeoSys, mit dem der Playground dieses Buches arbeitet (Kolditz u. a. 2012), pflegt eine große Sammlung von Benchmarks, die bei jeder Codeänderung automatisch nachgerechnet werden. Einer davon ist der Thiem-Brunnentest: Pumpt man aus einem Brunnen im Grundwasser, senkt sich der Wasserspiegel zu einem Trichter, dessen genaue Form Günther Thiem schon 1906 als exakte Formel angab (Thiem 1906). Die Software rechnet denselben Trichter numerisch und legt ihn über die analytische Kurve — genau unser Nulltest-Gedanke, nur mit Wasser statt Wärme (im Playground unter sim_thiem/). Wer so prüft, darf seinem Programm auch dort trauen, wo es keine Formel mehr gibt.
12.10 Die Auflösung des Aufhängers
Zurück zu den zwei Berichten. Büro A sagt 140 W, Büro B sagt 280 W — und die Temperaturbilder sind identisch. Jetzt haben wir die Werkzeuge. Der Fehler: Büro B hat sich bei der Wärmeleitfähigkeit um den Faktor zwei vertan — ein klassischer Einheiten-Salat, \(\lambda = 100\) statt \(50\ \mathrm{W/(m\,K)}\).
Warum ist das am Bild nicht zu sehen? Weil im quellfreien Träger die Temperaturgerade gar nicht von \(\lambda\) abhängt — die Wärmeleitfähigkeit kürzt sich aus der stationären Gleichung heraus. Ob \(\lambda\) nun 50 oder 100 ist: Die Lösung ist dieselbe Gerade 20 → 300. Darum bestehen beide Büros den Nulltest mit Bravour — das schärfste Werkzeug ist hier blind.
Code
import os, syssys.path.insert(0, os.path.join("..", "programme", "kap12"))from kap12_pruefstation import loese_stationaer, waermestrom_durch_traegerdef als_text(temperaturen):"""Rundet eine Temperaturliste zu lesbarem Text (ohne Comprehension).""" stuecke = []for wert in temperaturen: stuecke.append("%.1f"% wert)return" ".join(stuecke)richtig = loese_stationaer(5, 0.0, leitfaehigkeit=50.0)falsch = loese_stationaer(5, 0.0, leitfaehigkeit=100.0)print("Temperaturgerade bei lambda = 50:", als_text(richtig))print("Temperaturgerade bei lambda = 100:", als_text(falsch))print("-> identisch: der Nulltest schweigt.")print()print("Waermestrom lambda = 50 W/(m K):", "%.0f W"% waermestrom_durch_traeger(50.0))print("Waermestrom lambda = 100 W/(m K):", "%.0f W"% waermestrom_durch_traeger(100.0))
Temperaturgerade bei lambda = 50: 20.0 90.0 160.0 230.0 300.0
Temperaturgerade bei lambda = 100: 20.0 90.0 160.0 230.0 300.0
-> identisch: der Nulltest schweigt.
Waermestrom lambda = 50 W/(m K): 140 W
Waermestrom lambda = 100 W/(m K): 280 W
Interpretation: Der Fehler steckt nicht in der Temperatur, sondern im Wärmestrom — und der hängt von \(\lambda\) ab: \(Q = \lambda A\,\Delta T/L
= 50\cdot 0{,}01\cdot 280/1 = 140\,\mathrm{W}\) mit dem richtigen \(\lambda\), und glatt das Doppelte, 280 W, mit dem falschen. Die eine Kennzahl, die das schöne Bild nicht zeigt, ist genau die, in der der Fehler wohnt. Büro A hat recht. Und weil man den wahren Wärmestrom aus dem Brandherd unabhängig kennt (oder von Hand nachrechnet), fällt Büro B beim Punkt 5 der Checkliste sofort durch — die Plausibilität von Hand.
Das ist die Quintessenz des Kapitels: Ein FEM-Programm liefert immer ein plausibles, buntes Bild. Ob es stimmt, entscheidet nicht das Bild, sondern die Prüfung — Nulltest, Konvergenz, Bilanz, Wärmestrom, Handrechnung. Bunt ist kein Beweis.
TippMerkkasten
Verifikation (rechne ich richtig?) prüft man rechnerisch — mit Prüffällen, Konvergenz und Bilanzen; Validierung (rechne ich das Richtige?) braucht Messungen. Erst das eine, dann das andere.
Der Nulltest (die Gerade) muss die FEM exakt treffen — jede Abweichung ist ein Bug, kein Netzfehler. Sein blinder Fleck: Fehler, die die exakte Lösung nicht verändern.
Bei der Parabel viertelt sich der Fehler mit jeder Netzhalbierung (zweite Ordnung); in 1D sind die Knoten sogar exakt, der Fehler lebt dazwischen.
Zeit ist gröber als Raum: der implizite Euler halbiert den Fehler mit \(\Delta t\) nur (erste Ordnung) — Zeitschritt beim Verfeinern nicht vergessen.
Die Energiebilanz läuft immer mit: rein \(=\) gespeichert \(+\) raus, auf 1 % genau.
Bunt ist kein Beweis. Erst wenn die fünf Haken der Simulations-Checkliste sitzen, glaubt man dem Bild.
Roter Faden
Wo kam das schon vor, wo kommt es wieder? — Dieses Kapitel hat die Versprechen eingelöst, die der ganze Teil III und IV gestreut hat: die Gerade aus Kapitel 6 wird zum Nulltest, die Parabel aus Kapitel 6 und Übung 7.3 zum Diskretisierungs-Prüffall, die Fühler-Kurve aus Kapitel 5 entpuppt sich als die exakte Reihe, und die Randbedingungs-Checkliste aus Kapitel 8 wächst zur Simulations-Checkliste. Die Elementfamilien aus Kapitel 10 und Kapitel 11 treten gegeneinander an. — Nach vorn: Der Nulltest eröffnet Kapitel 13, wo freie Wärmedehnung Spannung null ergeben muss (Übung 12.3 entwirft ihn schon); und die Simulations-Checkliste hängt im Finale (Kapitel 15) neben den Reglern, wenn du selbst am vollständigen FEM-Programm drehst.
Übungen
Ü 12.1 (Verstehen). Sieh dir den Konvergenzplot (Abbildung 12.2) an. Der Zwischen-Knoten-Fehler der Parabel beträgt bei 17 Knoten 0,156 °C, bei 33 Knoten 0,039 °C. Du brauchst für deine Anwendung eine Genauigkeit von 0,1 °C. Wie viele Knoten reichen — und wie viele sind es ungefähr, wenn du den Viertelungs-Rhythmus zum Abschätzen nutzt?
HinweisMusterlösung zu Ü 12.1
17 Knoten (0,156 °C) reichen nicht ganz, 33 Knoten (0,039 °C) reichen locker — also liegt die Grenze dazwischen. Genauer: Der Fehler wächst mit \(h^2\), also \(\text{Fehler} \approx 40\cdot h^2\) (in °C, mit \(h\) in Metern). Für 0,1 °C: \(h^2 = 0{,}1/40 = 0{,}0025\), also \(h = 0{,}05\,\mathrm{m}\). Das sind \(L/h = 20\) Elemente, mithin 21 Knoten. Wer nur die tabellierten Netze 5/9/17/33 zur Hand hat, nimmt sicherheitshalber die 33 — lieber etwas zu fein als zu grob. Die Lehre: Man liest die nötige Feinheit direkt aus der Konvergenzkurve ab, statt zu raten.
Ü 12.2 (Verändern). In der Prüfstation (bzw. im Skript loesungen/kap12_ue2.py) spielst du Zeitschritt und Maschenweite gegeneinander aus. Verfeinere beim transienten Fall nur das Netz und lass \(\Delta t = 200\,\mathrm{s}\) fest. Was passiert mit dem Fehler ab einer gewissen Knotenzahl — und warum? Wann lohnt feineres Netz nichts mehr?
HinweisMusterlösung zu Ü 12.2
Der Gesamtfehler ist ungefähr die Summe aus einem Raum-Anteil (viertelt sich mit \(h\)) und einem Zeit-Anteil (halbiert sich mit \(\Delta t\)). Verfeinert man nur das Netz, fällt der Fehler zunächst — bis er auf dem Zeit-Anteil aufsitzt und sich durch noch mehr Knoten nicht weiter drücken lässt: Bei \(\Delta t = 200\,\mathrm{s}\) bleibt er ab etwa 17 Knoten bei rund 7,4 °C stehen, ganz gleich, wie fein das Netz noch wird (Skript kap12_ue2.py, Teil a). Erst ein kleineres \(\Delta t\) senkt ihn weiter (Teil b, brave Halbierung). Merke: Feineres Netz und kleineres \(\Delta t\) müssen zusammen gehen — wer nur eines verfeinert, läuft in den Boden des anderen. Zeit ist gröber als Raum.
Ü 12.3 (Übertragen). Entwirf einen Nulltest für die Thermomechanik (Vorgriff auf Kapitel 13). Ein frei aufgehängter Träger (nirgends festgeklemmt) wird gleichmäßig erwärmt. Er dehnt sich aus — aber welche Spannung muss dabei im Inneren herrschen? Beschreibe das Prüf-Protokoll: Was ist der bekannte Sollwert, und was schließt du aus einer Abweichung?
HinweisMusterlösung zu Ü 12.3
Sollwert: Spannung exakt null, überall. Ein frei beweglicher Körper, der sich gleichmäßig erwärmt, dehnt sich einfach ungehindert aus — er darf länger werden, ohne dass irgendwo eine Kraft entsteht. Spannung baut sich erst auf, wenn die Ausdehnung behindert wird (eingeklemmte Enden). Protokoll: (1) Träger frei lagern (keine Einspannung, aber gerade so viele Lager, dass er nicht davonschwebt — statisch bestimmt). (2) Gleichmäßige Temperaturerhöhung \(\Delta T\) aufbringen. (3) Die berechnete Spannung an allen Knoten ablesen. (4) Erwartung: null (bis auf Rechenungenauigkeit). Jede von null verschiedene Spannung ist ein Bug — meist eine versehentlich zu starre Lagerung, ein falscher Ausdehnungskoeffizient oder ein Vorzeichenfehler in der Kopplung Temperatur → Dehnung. Das ist derselbe Nulltest-Gedanke wie bei der Geraden: eine Aufgabe mit exakt bekannter, im Ansatzraum liegender Lösung, an der ein Fehler garantiert auffällt. Kapitel 13 baut ihn wirklich.
Das Kleingedruckte
Dieses Kapitel verifiziert, es validiert nicht. Wir haben die FEM gegen exakte Lösungen derselben Gleichung gehalten — aber nie gefragt, ob diese Gleichung die Wirklichkeit trifft. Der echte Stahlträger gibt Wärme an die Luft ab (das haben wir weggelassen), sein \(\lambda\) ändert sich mit der Temperatur (wir halten es konstant), und bei mehreren hundert Grad spielen Effekte hinein — allen voran der Festigkeitsverlust des Stahls —, die unser Modell gar nicht kennt. Solche Fragen beantwortet nur die Validierung — der Vergleich mit einer echten Messung. Sie ist der ehrliche nächste Schritt, den dieses Buch nicht mehr geht.
Zwei weitere Grenzen: Wir haben keine Fehlerabschätzungs-Theorie getrieben (keine Normen, keine Beweise — nur den handfesten Viertelungs-Rhythmus), und wir haben die Netze von Hand verfeinert. Professionelle Programme verfeinern adaptiv: Sie schätzen selbst, wo der Fehler groß ist (etwa an der wandernden Front), und setzen genau dort mehr Knoten. Das ist eine eigene Kunst — für uns bleibt es ein Namedrop am Wegesrand.
Kolditz, Olaf, Sebastian Bauer, Lars Bilke, Norbert Böttcher, Jens-Olaf Delfs, u. a. 2012. „OpenGeoSys: an open-source initiative for numerical simulation of thermo-hydro-mechanical/chemical (THM/C) processes in porous media“. Environmental Earth Sciences 67 (2): 589–99. https://doi.org/10.1007/s12665-012-1546-x.
Roache, Patrick J. 1998. Verification and Validation in Computational Science and Engineering. Hermosa Publishers.