numerische Differentiation

[rayman]

Ich hab mal aus Langeweile ein klitzekleines (überhaput noch nicht fertiges) Programm zur Approximation reelwertiger Funktionen mittels TAYLOR-Polynomen geschrieben.

Dafür muss man ja die Funktion an der entsprechenden Stelle mehrmals differenzieren (je nachdem, wie genau man das haben will). Ich hab dazu zuerst eine rekursive Function foo_ gemacht, die einfach einen Differenzenquotienten bildet:

Delphi-Quellcode:

function foo_(x: Double; degree: integer): Double;
begin
if degree = 0
  then Result := foo(x)
  else Result := (foo_(x + PRECISION, degree - 1) - foo_(x - PRECISION, degree - 1)) / (2*PRECISION);
end;

Nun hab ich gemerkt, dass die für eine 10. Ableitung schon ziemlich lange braucht (logisch: rekursiv --> ca. 2^10 Aufrufe). Also wollte ich eine iterative machen, die ein array of Double benutzt um Schicht für Schicht die entsprechende Ableitung zu berechnen:

Delphi-Quellcode:

function foo_(x: Double; degree: integer): Double;
var
  diff: array of Double; // für die Differenzenquotienten
  n, m: integer;
begin
SetLength(diff, degree + 1);

for n := 0 to degree do
  diff[n] := foo(x + (-degree/2 + n) * PRECISION);

for n := 0 to degree do
  for m := 0 to n - 1 do
    diff[m] := (diff[m+1] - diff[m]) / PRECISION;

Result := diff[0];
end;

Für die n-te Ableitung suche ich mir n+1 gleichmäßig um die Stelle x verteilte Funktionswerte. Aus zwei benachbarten Stellen lässt sich jeweils ein Differenzenquotient berechnen, der dann wieder in das array kommt. Für den nächsten Schritt hat man also einen Wert weniger. Nach n solchen Durchläufen hat man dann die geünschte n-te Ableitung an x. Insgesamt brauch man nur etwa 1/2 * n² Differenzenquotienten zu machen - geht also sehr viel schneller.

So. Problem:
Leider erfüllt diese 2. Function ihre Aufgabe nur für 1.Ableitungen. Bei der 2. kommt immer ein ganz komischer Wert raus. Einfach mal angucken. Ich hab alles dreiviermal geprüft, und weiß nicht, wo der Fehler liegt.

Kann das irgendwie an dem Zahlenformat (Double) liegen? Die rekursive nutzt doch im Prinzip die gleichen Werte, nur dass sie dort mehrfach berechnet werden.
Bin ich mit meiner Methode auf dem Holzweg?
Oda was?

Wäre echt gut, wenn ihr mir da helfen könntet, weil ich jetz momentan nich weiterweiß.