Worum geht es?
Die Animationen zeigen eine konfigurierte kinematische Szene in einem inertialen Bezugsystem S sowie dann die entsprechende lorentztransformierte Szene in S'. Das ist soweit richtig, nur: Die kinematischen Szenen sind sozusagen im galileischen Sinne konfiguriert und nicht im einsteinschen Sinne. Bei Galilei sind z.B. beide Panzer in S gleich lang, bei Einstein aber nicht.
Die neue Version der Software, die ScænaSRT heißt, berücksichtigt nun diesen Umstand. Die kinematische Szene wird für ein `ruhendes' galileisches inertiales Bezugsystem G konfiguriert. Falls die Szene auch Lichtausbreitung beschreibt, ist das das eine absolute Bezugsystem, in dem sich das Licht ausbreitet. Die Software erstellt dann die Animation mit folgenden vier Darstellungen:
`Ruhendes' galileisches Inertialsystem G
Zeigt die konfigurierte kinematische Szene.
`Bewegtes' galileisches Inertialsystem G'
Zeigt dieselbe Szene aus G in einem relativ bewegten inertialen Bezugssystem G',
indem die Darstellung in G' aus der Darstellung in G per Galileitransformation gewonnen wird.
`Ruhendes' einsteinsches Inertialsystem E
Zeigt dieselbe Szene aus G wiederum in `demselben' `ruhenden' Bezugssystem, aber mit `Einstein-Semantik'. D.h. bewegte Objekte erscheinen in Bewegungsrichtung lorentzkontrahiert. Um die Darstellung im `ruhenden' Bezugssystem E aus der Darstellung im `ruhenden' Bezugssystem G zu gewinnen, werden die Trajektorien eines jeden bewegten Vertex (der kein Licht repräsentiert) zu jedem Zeitpunkt t in je eigener Weise transformiert. Der Vertex wird zunächst per Galileitransformation in ein inertiales Bezugssystem überführt, in dem er ruht, dann wird zwecks Kontraktion der Gammafaktor angewendet und dann die inverse Galileitransformation. (Die kinematischen Szenen basieren inzwischen ausschließlich auf einer Sammlung differenzierbarer Funktionen, so dass die Geschwindigkeiten effizient und hochgenau zur Verfügung stehen.)
`Bewegtes' einsteinsches Inertialsystem E'
Zeigt die Szene aus E in einem relativ bewegten inertialen Bezugssystem E', indem die Darstellung in E' aus der Darstellung in E per vollständiger Lorentztransformation gewonnen wird.
Auf dieser Seite habe ich die mathematischen Zusammenhänge genauer beschrieben.
Das folgende einfache Beispiel demonstriert, worum es geht:
- Im `ruhenden' System G sind drei baugleiche Panzer zu sehen, die sich in verschiedene Richtungen bewegen.
- Im `bewegten' System G' die galieitransformierte Szene dazu.
- Im `ruhenden' System E dieselbe Szene, aber mit lorentzkontrahierten Objekten.
- Im `bewegten' System E' die lorentztransformierte Szene zu E.
Für die meisten bisher gezeigten Bilder ergeben sich nur unbedeutende Änderungen, die das Wesentliche des Ergebnisses nicht betreffen. Ein wesentlicher Unterschied ergibt sich aber für das rotierende Zahnrad:
Erklärungen zu dem Bild habe ich im Strang Zum `Ehrenfest-Paradoxon' (Zur Zahnradbahn) formuliert.
Die weiteren Animationen kann man sich hier anschauen:
Zahnrad mit Zahnschienen
Zum Phipps-Versuch [1]
Zum Phipps-Versuch [2]
Ausbreitung dreier Lichtblitze
Licht vs. Ballon
Reisende Zwillinge
Zwillinge mit Uhren
Zum symmetrischen Panzerparadoxon
Gruß
Faber
P.S.: Ich werde die API und Beispiele (Java 1.6 eclipse Projekt, nothing left to buy) voraussichtlich veröffentlichen, wenn ich die Sache gänzlich dokumentiert und weitere Vereinfachungen implementiert habe. Falls Interesse an Einsicht in die Funktionsweise des Kerns der API besteht, hier Auszüge aus den wichtigsten Quelltexten:
com.mahag.neufor.faber.srt.trajectory.Trajectory.java
- Code: Alles auswählen
package com.mahag.neufor.faber.srt.trajectory;
import java.util.HashMap;
import java.util.Map;
/**
* Abstrakte Basisklasse für transformierbare Trajektorien.
* <p>
* Es wird empfohlen, möglichst {@link GeneralTrajectory} zu verwenden.
* </p>
*
* @author Faber@www.mahag.com/neufor
*/
public abstract class Trajectory
{
/** Transformationskontext. */
protected final TransformationContext tc;
/** Startzeitpunkt der Trajektorie in S. */
protected final double t0;
/** Endzeitpunkt der Trajektorie in S. */
protected final double t1;
/** Vorberechnete Werte. */
protected Map<Double, Point3D> valueMap = null;
/** Transformierte Trajektorie. */
protected Trajectory einsteinRestTrajectory = null;
/** Transformierte Trajektorie. */
protected Trajectory galileiMovedTrajectory = null;
/** Transformierte Trajektorie. */
protected Trajectory einsteinMovedTrajectory = null;
/**
* Liefert neue Instanz.
*
* @param tc
* Transformationskontext.
* @throws IllegalTrajectoryException
* falls Trajektorie nicht gültig ist (schneller als Lichtgeschwindigkeit).
*/
protected Trajectory(final TransformationContext tc) throws IllegalTrajectoryException
{
this(tc, tc.t0, tc.t1);
}
/**
* Liefert neue Instanz.
*
* @param tc
* Transformationskontext.
* @param t0
* Startzeitpunkt der Trajektorie.
* @param t1
* Endzeitpunkt der Trajektorie.
* @throws IllegalTrajectoryException
* falls Trajektorie nicht gültig ist (schneller als Lichtgeschwindigkeit).
*/
protected Trajectory(final TransformationContext tc, final double t0, final double t1)
throws IllegalTrajectoryException
{
this.tc = tc;
this.t0 = t0;
this.t1 = t1;
}
/**
* Liefert den Transformationskontext.
*/
public TransformationContext getTransformationContext()
{
return tc;
}
/**
* Zeigt an, ob {@code t} im Bereich der Lebensdauer der Trajektorie liegt.
*/
public boolean isAliveAt(final double t)
{
return t0 - tc.eps / tc.c <= t && t <= t1 + tc.eps / tc.c;
}
/**
* Liefert {@code [ x(t), y(t), z(t) ]}.
*
* @throws IllegalTrajectoryException
* falls Trajektorie nicht gültig ist.
*/
public Point3D getVertex(final double t) throws IllegalTrajectoryException
{
if (valueMap == null)
{
valueMap = new HashMap<Double, Point3D>();
for (int k = 0, N = tc.getNoOfPictures(); k < N; k++)
{
final double tau = tc.getT(k);
valueMap.put(tau, calculateVertex(tau));
}
}
Point3D vertex = valueMap.get(t);
if (vertex == null)
{
vertex = calculateVertex(t);
}
return vertex;
}
/**
* Berechnet und liefert {@code [ x(t), y(t), z(t) ]}.
* <p>
* Wird für kein {@code t} mehr als einmal aufgerufen.
* </p>
*
* @param t
* Zeitpunkt t.
*/
protected abstract Point3D calculateVertex(final double t);
/**
* Liefert die zeitliche Ableitung zum Zeitpunkt {@code t}.
*/
public abstract Point3D getDerivative(final double t);
/**
* Liefert die Schnelligkeit der Trajektorie zum Zeitpunkt {@code t}.
*/
public double getSpeed(final double t)
{
return getDerivative(t).getDistance();
}
/**
* Liefert den Parameter {@code t}.
*
* @param tMoved
* Parameter {@code t'}.
* @return Parameter {@code t}.
*/
public double getTRest(final double tMoved)
{
final double t = tMoved;
return t;
}
/**
* Liefert den Parameter {@code tStrich}.
*
* @param tRest
* Parameter {@code t}.
* @return Parameter {@code t'}.
*/
public double getTMoved(final double tRest)
{
final double tStrich = tRest;
return tStrich;
}
/**
* Liefert neue transformierte Trajektorie für die angegebene Darstellung.
*
* @param systemType
* Art der Darstellung.
* @return Transformierte Trajektorie für die angegebene Darstellung (kann {@code this} sein).
*/
public Trajectory to(final SystemType systemType)
{
switch (systemType)
{
default :
case GALILEI_REST :
return this;
case GALILEI_MOVED :
return toGalileiMoved();
case EINSTEIN_REST :
return toEinsteinRest();
case EINSTEIN_MOVED :
return toEinsteinRest().toEinsteinMoved();
}
}
/**
* Liefert neue transformierte Trajektorie für das Einstein-Ruhesystem.
*/
private Trajectory toEinsteinRest()
{
class EinsteinRestTrajectory extends Trajectory
{
public EinsteinRestTrajectory(final TransformationContext ctx)
{
super(ctx, Trajectory.this.t0, Trajectory.this.t1);
}
/*
* @see com.mahag.neufor.faber.srt.trajectory.Trajectory#calculateVertex(double)
*/
@Override
public Point3D calculateVertex(final double tRest)
{
final double t = tRest;
final double c = tc.c;
final Point3D p = Trajectory.this.getVertex(t);
if (p.z != 0)
{
throw new UnsupportedOperationException("3D not yet implemented");
}
final Point3D dp = Trajectory.this.getDerivative(t);
final double w = dp.getDistance();
if (Math.abs(w - c) < tc.eps / c)
{
return p;
}
final double wx = dp.x;
final double wy = dp.y;
final double phi = Math.atan2(wy, wx);
final double cosPhi = Math.cos(phi);
final double sinPhi = Math.sin(phi);
final double scale = Math.sqrt(Math.abs(1 - w * w / (c * c)));
final double x0 = p.x;
final double y0 = p.y;
final double x1 = x0 * cosPhi + y0 * sinPhi;
final double y1 = -x0 * sinPhi + y0 * cosPhi;
final double x2 = scale * (x1 - w * t) + w * t;
final double y2 = y1;
final double x3 = x2 * cosPhi - y2 * sinPhi;
final double y3 = x2 * sinPhi + y2 * cosPhi;
return new Point3D(x3, y3);
}
/*
* @see com.mahag.neufor.faber.srt.trajectory.Trajectory#getDerivative(double)
*/
@Override
public Point3D getDerivative(double tRest)
{
return getNumericDerivative(tRest);
}
}
if (einsteinRestTrajectory == null)
{
einsteinRestTrajectory = new EinsteinRestTrajectory(tc);
}
return einsteinRestTrajectory;
}
/**
* Liefert neue transformierte Trajektorie für das bewegte Galilei-System.
*/
private Trajectory toGalileiMoved()
{
class GalileiMovedTrajectory extends Trajectory
{
public GalileiMovedTrajectory(final TransformationContext ctx)
{
super(ctx, Trajectory.this.t0, Trajectory.this.t1);
}
/*
* @see com.mahag.neufor.faber.srt.trajectory.Trajectory#calculateVertex(double)
*/
@Override
public Point3D calculateVertex(final double tMoved)
{
final double t = getTRest(tMoved);
final Point3D p = Trajectory.this.getVertex(t);
final double x = p.x - tc.v * t;
final double y = p.y;
final double z = p.z;
return new Point3D(x, y, z);
}
/*
* @see com.mahag.neufor.faber.srt.trajectory.Trajectory#getDerivative(double)
*/
@Override
public Point3D getDerivative(double tMoved)
{
return getNumericDerivative(tMoved);
}
}
if (galileiMovedTrajectory == null)
{
galileiMovedTrajectory = new GalileiMovedTrajectory(tc);
}
return galileiMovedTrajectory;
}
/**
* Liefert neue transformierte Trajektorie für das bewegte Einstein-System.
*/
private Trajectory toEinsteinMoved()
{
class EinsteinMovedTrajectory extends Trajectory
{
private double previousTRest = Double.NaN;
protected final Map<Double, Double> tRestOfTMoved = new HashMap<Double, Double>();
public EinsteinMovedTrajectory(final TransformationContext ctx, final double t0,
final double t1)
{
super(ctx, t0, t1);
}
/*
* @see com.mahag.neufor.faber.srt.trajectory.Trajectory#calculateVertex(double)
*/
@Override
public Point3D calculateVertex(final double tMoved)
{
final Point3D vertex = valueMap.get(tMoved);
if (vertex != null)
{
return vertex;
}
final double t = getTRest(tMoved);
final Point3D p = Trajectory.this.getVertex(t);
final double x = tc.gamma * (p.x - tc.v * t);
final double y = p.y;
final double z = p.z;
return new Point3D(x, y, z);
}
/*
* @see com.mahag.neufor.faber.srt.trajectory.Trajectory#getTRest(double)
*/
@Override
public double getTRest(final double tMoved)
{
if (tMoved <= tc.t0)
{
previousTRest = tc.t0;
}
else if (tMoved > tc.t1)
{
previousTRest = tc.t1;
}
else
{
final double previousTMoved = tc.floor(tMoved);
if (previousTMoved != tMoved)
{
final Double precalculatedTRest = tRestOfTMoved.get(previousTMoved);
if (precalculatedTRest != null)
{
previousTRest = precalculatedTRest;
}
}
else
{ // valueMap wird gerade berechnet
if (Double.isNaN(previousTRest))
{
previousTRest = Trajectory.this.t0;
}
}
}
final Double precalculatedT = tRestOfTMoved.get(tMoved);
final double t = precalculatedT != null ? precalculatedT
: (previousTRest = new LorentzTimeMapper(Trajectory.this).getT(tMoved,
previousTRest));
tRestOfTMoved.put(tMoved, t);
return t;
}
/*
* @see com.mahag.neufor.faber.srt.trajectory.Trajectory#getTMoved(double)
*/
@Override
public double getTMoved(final double tRest)
{
return tc.gamma * (tRest - tc.v / (tc.c * tc.c) * Trajectory.this.getVertex(tRest).x);
}
/*
* @see com.mahag.neufor.faber.srt.trajectory.Trajectory#getDerivative(double)
*/
@Override
public Point3D getDerivative(double tMoved)
{
final double t0 = tc.floor(tMoved);
final double t1 = tc.ceiling(tMoved);
final Point3D p0 = getVertex(t0);
final Point3D p1 = getVertex(t1);
final double dx = (p1.x - p0.x) / (t1 - t0);
final double dy = (p1.y - p0.y) / (t1 - t0);
final double dz = (p1.z - p0.z) / (t1 - t0);
return new Point3D(dx, dy, dz);
}
}
if (einsteinMovedTrajectory == null)
{
final double t0Moved = tc.gamma * t0;
final double t1Moved = tc.gamma * t1;
einsteinMovedTrajectory = new EinsteinMovedTrajectory(tc, t0Moved, t1Moved);
}
return einsteinMovedTrajectory;
}
/**
* Liefert numerisch ermittelte zeitliche Ableitung zum Zeitpunkt {@code t}.
* <p>
* Wird nur verwendet, wo nichts besseres zur Verfügung steht.
* </p>
*/
protected Point3D getNumericDerivative(final double t)
{
final double t0 = t - 0.5 * tc.eps * tc.c;
final double t1 = t + 0.5 * tc.eps * tc.c;
final Point3D p0 = getVertex(t0);
final Point3D p1 = getVertex(t1);
final double dx = (p1.x - p0.x) / (t1 - t0);
final double dy = (p1.y - p0.y) / (t1 - t0);
final double dz = (p1.z - p0.z) / (t1 - t0);
return new Point3D(dx, dy, dz);
}
}
com.mahag.neufor.faber.srt.trajectory.GeneralTrajectory.java
- Code: Alles auswählen
package com.mahag.neufor.faber.srt.trajectory;
import com.mahag.neufor.faber.srt.function.Add;
import com.mahag.neufor.faber.srt.function.Function;
import com.mahag.neufor.faber.srt.function.Const;
import com.mahag.neufor.faber.srt.function.Mult;
/**
* Trajektorie, die durch drei ableitbare, skalare Funktionen der Zeit definiert ist.
* <p>
* Es wird empfohlen, {@link Function}en zu verwenden, die (möglichst) exakte Ableitungen liefern.
* </p>
*
* @author Faber@www.mahag.com/neufor
*/
public class GeneralTrajectory extends Trajectory
{
/** x-Koordinate. */
private final Function x;
/** y-Koordinate. */
private final Function y;
/** z-Koordinate. */
private final Function z;
/**
* Liefert neue Instanz.
*
* @param ctx
* Transformationskontext.
*/
public GeneralTrajectory(final TransformationContext ctx)
{
this(ctx, 0, 0);
}
/**
* Liefert neue Instanz.
*
* @param ctx
* Transformationskontext.
* @param x
* x-Koordinate.
* @param y
* y-Koordinate
*/
public GeneralTrajectory(final TransformationContext ctx, final double x, final double y)
{
this(ctx, ctx.t0, ctx.t1, x, y);
}
/**
* Liefert neue Instanz.
*
* @param ctx
* Transformationskontext.
* @param T0
* Startzeitpunkt der Trajektorie.
* @param T1
* Endzeitpunkt der Trajektorie.
* @param x
* x-Koordinate.
* @param y
* y-Koordinate
*/
public GeneralTrajectory(final TransformationContext ctx, final double T0, final double T1,
final double x, final double y)
{
this(ctx, T0, T1, new Const(x), new Const(y));
}
/**
* Liefert neue Instanz.
*
* @param ctx
* Transformationskontext.
* @param x
* x-Koordinate.
* @param y
* y-Koordinate
*/
public GeneralTrajectory(final TransformationContext ctx, final Function x, final Function y)
{
this(ctx, ctx.t0, ctx.t1, x, y);
}
/**
* Liefert neue Instanz.
*
* @param ctx
* Transformationskontext.
* @param T0
* Startzeitpunkt der Trajektorie.
* @param T1
* Endzeitpunkt der Trajektorie.
* @param x
* x-Koordinate.
* @param y
* y-Koordinate
*/
public GeneralTrajectory(final TransformationContext ctx, final double T0, final double T1,
final Function x, final Function y)
{
super(ctx, T0, T1);
this.x = x;
this.y = y;
this.z = new Const(0);
}
/**
* Liefert neue Instanz.
*
* @param ctx
* Transformationskontext.
* @param T0
* Startzeitpunkt der Trajektorie.
* @param T1
* Endzeitpunkt der Trajektorie.
* @param x
* x-Koordinate.
* @param y
* y-Koordinate
* @param z
* z-Koordinate
*/
public GeneralTrajectory(final TransformationContext ctx, final double T0, final double T1,
final Function x, final Function y, final Function z)
{
super(ctx, T0, T1);
this.x = x;
this.y = y;
this.z = z;
}
/**
* Liefert neue Instanz {@code [ x + vx * t; y + vx * t; 0 ]}.
*
* @param ctx
* Transformationskontext.
* @param T0
* Startzeitpunkt der Trajektorie.
* @param T1
* Endzeitpunkt der Trajektorie.
* @param x
* Erster x-Koeffizient.
* @param y
* Erster y-Koeffizient.
* @param vx
* Zweiter x-Koeffizient.
* @param vy
* Zweiter y-Koeffizient.
*/
public GeneralTrajectory(final TransformationContext ctx, final double T0, final double T1,
final double x, final double y, final double vx, final double vy)
{
this(ctx, T0, T1, new Const(x), new Const(y), new Const(vx), new Const(vy));
}
/**
* Liefert neue Instanz {@code [ x + vx * t; y + vx * t; 0 ]}.
*
* @param ctx
* Transformationskontext.
* @param T0
* Startzeitpunkt der Trajektorie.
* @param T1
* Endzeitpunkt der Trajektorie.
* @param x
* Erster x-Koeffizient.
* @param y
* Erster y-Koeffizient.
* @param vx
* Zweiter x-Koeffizient.
* @param vy
* Zweiter y-Koeffizient.
*/
public GeneralTrajectory(final TransformationContext ctx, final double T0, final double T1,
final Function x, final Function y, final Function vx, final Function vy)
{
this(ctx, T0, T1, new Add(x, new Mult(vx, Function.T)), new Add(y, new Mult(vy, Function.T)),
Function.NIL);
}
/*
* @see com.mahag.neufor.faber.srt.trajectory.Trajectory#calculateVertex(double)
*/
@Override
protected Point3D calculateVertex(final double t)
{
return new Point3D(x.value(t), y.value(t), z.value(t));
}
/*
* @see com.mahag.neufor.faber.srt.trajectory.Trajectory#getDerivative(double)
*/
@Override
public Point3D getDerivative(final double t)
{
final double x = this.x.derivative().value(t);
final double y = this.y.derivative().value(t);
final double z = this.z.derivative().value(t);
return new Point3D(x, y, z);
}
}
com.mahag.neufor.faber.srt.trajectory.LorentzTimeMapper.java
- Code: Alles auswählen
package com.mahag.neufor.faber.srt.trajectory;
import com.mahag.neufor.faber.srt.function.Function;
/**
* Bildet für eine gegebene Trajektorie die t'-Achse von S' auf die t-Achse von S ab, wobei
* Mehrdeutigkeiten mithilfe der Forderung aufgelöst werden, dass t(t') stetig ist.
*
* @author Faber@www.mahag.com/neufor
*/
class LorentzTimeMapper
{
/** Ausgabe von Debug-Infos ein-/ausschalten. */
private static final boolean DEBUG = false;
/** Lorentz-Kontext. */
private final TransformationContext tc;
/** Trajektorie. */
private final Trajectory trajectory;
/**
* Liefert neue Instanz.
*
* @param trajectory
* Trajektorie.
*/
public LorentzTimeMapper(final Trajectory trajectory)
{
this.tc = trajectory.getTransformationContext();
this.trajectory = trajectory;
}
/**
* Berechnet und liefert {@code t(t')} mit {@code t'(t) = gamma * (t - v / (c * c) * x(t))}.
* <p>
* Falls t'(t) nicht eindeutig ist, wird der t-Wert geliefert, für den {@code t(t')} stetig ist.
* Diese Methode sollte sukkzessive von t0Strich an für tStrich-Werte im Intervall
* [t0Strich;t1Strich] aufgerufen werden.
* </p>
*
* @param tStrich
* Zeitwert {@code t'}.
* @param previousT
* Der zuletzt errechnete {@code t}-Wert oder irgendein bekannter {@code t}-Wert
* möglichst in der Nähe von dem gesuchten {@code t}-Wert. Ist er weit weg, kann es
* dauern, bis die Methode zurückkehrt.
* @return {@code t(t')}
* @throws NumericalException
* wenn der t-Wert nicht berechnet werden konnte, was auf eine entartete Trajektorie
* hinweist.
*/
public double getT(final double tStrich, double previousT) throws NumericalException
{
final Function tStrichMinusTStrichOfT = new Function()
{
@Override
public double value(double t)
{
return tStrich - tc.gamma * (t - tc.v / (tc.c * tc.c) * trajectory.getVertex(t).x);
}
@Override
public Function derivative()
{
throw new IllegalStateException("unreached");
}
};
final Function tStrichOfT = new Function()
{
@Override
public double value(double t)
{
return tc.gamma * (t - tc.v / (tc.c * tc.c) * trajectory.getVertex(t).x);
}
@Override
public Function derivative()
{
throw new IllegalStateException("unreached");
}
};
// TODO: die 20 irgendwie global angebbar machen
final double dt = 20 / tc.zoom * tc.gamma * (tc.t1 - tc.t0);
final double aMin = tc.t0 - dt;
final double bMax = tc.t1 + dt;
final double T = 0.1 * tc.dt;
double a = previousT;
double b = a + T;
double tStrichA = tStrichOfT.value(a);
double tStrichB = tStrichOfT.value(b);
if (DEBUG)
{
System.out.format("continuous bracketing:\n\n");
System.out.format("tStrichA = %15.12f [%15.12f]\n", tStrichA, a);
System.out.format("tStrich = %15.12f\n", tStrich);
System.out.format("tStrichB = %15.12f [%15.12f]\n", tStrichB, b);
System.out.format("\n");
}
while (true)
{
if (a < aMin)
{
throw new NumericalException("continuous bracketing failed [a < aMin]");
}
if (b > bMax)
{
throw new NumericalException("continuous bracketing failed [b > bMax]");
}
if (tStrichA - tc.eps <= tStrich && tStrich <= tStrichB + tc.eps)
{
break;
}
if (tStrichB - tc.eps <= tStrich && tStrich <= tStrichA + tc.eps)
{
break;
}
final boolean moveUp = tStrich > tStrichB;
if (moveUp)
{
a = b;
b += T;
tStrichA = tStrichB;
tStrichB = tStrichOfT.value(b);
}
else
{
b = a;
a -= T;
tStrichB = tStrichA;
tStrichA = tStrichOfT.value(a);
}
}
final double t = solve_VanWijngaardenDekkerBrent(tStrichMinusTStrichOfT, a, b);
if (DEBUG)
{
System.out.format("tStrichA = %15.12f [%15.12f]\n", tStrichA, a);
System.out.format("tStrich = %15.12f [%15.12f]\n", tStrich, t);
System.out.format("tStrichB = %15.12f [%15.12f]\n", tStrichB, b);
System.out.format("\n");
}
return t;
}
/**
* Liefert Nullstelle in einem Intervall, falls die Intervallgrenzen ungleiches Vorzeichen haben,
* oder {@code Double.NaN} sonst.
* <p>
* Siehe William H. Press <i>et. al.</i>:
* "Numerical Recipes in C : the art of scientific computing" 2nd edition, Cambridge University
* Press 1992, Abschnitt 9.3, Seite359.
* </p>
*
* @param f
* Funktion.
* @param a
* Untere Intervallgrenze.
* @param b
* Obere Intervallgrenze.
* @return Nullstelle oder {@code NaN}, falls die Intervallgrenzen das gleiche Vorzeichen haben.
* @throws NumericalException
* Falls das Intervall ungültig ist.
*/
public double solve_VanWijngaardenDekkerBrent(final Function f, double a, double b)
throws NumericalException
{
// TODO: Toleranzen korrigieren. Steigt wegen ITMAX aus. Das betrifft die Effizienz.
if (DEBUG)
{
System.out.format("VanWijngaardenDekkerBrent:\n", a);
System.out.format(" a = %12.8f\n", a);
System.out.format(" b = %12.8f\n", b);
System.out.format("\n");
}
final double epsF = tc.eps;
final double epsT = epsF / tc.c;
final int ITMAX = 100;
final double EPS = 3E-8;
double c = b;
double d = 0;
double e = 0;
double min1;
double min2;
double fa = f.value(a);
double fb = f.value(b);
double fc = fb;
double p;
double q;
double r;
double s;
if (Math.abs(fa) < epsF)
{
return a;
}
if (Math.abs(fb) < epsF)
{
return b;
}
if (a >= b)
{
throw new NumericalException("illegal interval");
}
if ((fa > 0) == (fb > 0))
{
return Double.NaN;
}
for (int iter = 1; iter <= ITMAX; iter++)
{
if ((fb > 0 && fc > 0) || (fb < 0 && fc < 0))
{
c = a;
fc = fa;
e = d = b - a;
}
if (Math.abs(fc) < Math.abs(fb))
{
a = b;
b = c;
c = a;
fa = fb;
fb = fc;
fc = fa;
}
final double tol1 = 2 * EPS * Math.abs(b) + 0.5 * epsT;
final double xm = 0.5 * (c - b);
if (fb == 0 || (Math.abs(xm) <= tol1 && Math.abs(fb) <= epsF))
{
return b;
}
if (Math.abs(e) >= tol1 && Math.abs(fa) > Math.abs(fb))
{
s = fb / fa;
if (a == c)
{
p = 2 * xm * s;
q = 1 - s;
}
else
{
q = fa / fc;
r = fb / fc;
p = s * (2 * xm * q * (q - r) - (b - a) * (r - 1));
q = (q - 1) * (r - 1) * (s - 1);
}
if (p > 0)
{
q = -q;
}
p = Math.abs(p);
min1 = 3 * xm * q - Math.abs(tol1 * q);
min2 = Math.abs(e * q);
if (2 * p < (min1 < min2 ? min1 : min2))
{
e = d;
d = p / q;
}
else
{
d = xm;
e = d;
}
}
else
{
d = xm;
e = d;
}
a = b;
fa = fb;
if (Math.abs(d) > tol1)
{
b += d;
}
else
{
b += xm > 0 ? Math.abs(tol1) : -Math.abs(tol1);
}
fb = f.value(b);
}
if (DEBUG)
{
System.out.format("VanWijngaardenDekkerBrent [iterations exceeded]:\n", a);
System.out.format(" a = %12.8f\n", a);
System.out.format(" b = %12.8f\n", b);
}
return 0.5 * (a + b);
}
}