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ß.

[rayman]

Darf ich das mal pushen? Gestern abend hat das wohl keiner mitbekommen.
Bitte helft mir!

[Der_Unwissende]

Hi,
ich bin mir nicht ganz sicher, aber kannst du nicht deine Rekursive Formel sehr leicht in eine lineare umwandeln?

Zitat von rayman:

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;

Weiß natürlich nicht was foo macht, werd vielleicht mal in deinen Code (war doch welcher angehangen?) schauen. Aber die foo_ (gute Benennung), kannst du denke ich auch so schreiben:

Delphi-Quellcode:

function foo_(x : Double; degree: Integer): Double;
begin
  if degree = 0 then
    begin
      result := foo(x);
    end
  else 
    begin
      result := (foo(x + ((degree - 1) * PRECISION)) / degree * PRECISION) - (foo(x - ((degree - 1) * PRECISION)) / degree * PRECISION); 
    end; 
end;

Hoffe Klammerung stimmt, sorry hab's mir nicht genauer angeschaut.

Gruß Der Unwissende

[rayman]

Ich kanns jetz zwar nicht überprüfen, aber das dürfte imho so nicht hinhauen.

Du vergrößerst mit (degree * PRECISION) ja einfach das Anstiegsdreieck für den Quotienten. Dann hat man jedesmal eine erste Ableitung. Nur das die mit größerem degree eben ungenauer wird.

[rayman]

*push* *heul*

[tigerman33]

Ich schmeiß hier jetzt mal vollkommen ungetestet wilde Behauptungen in den Raum:

Zitat:

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;

Erstens: In der zweiten for-Schleife kannst du dir eigentlich den Durchlauf für n=0 sparen. Da liegt aber sicher nicht der Fehler.
Aber: In der inneren Schleife (for m := ...) ist IMHO die Obergrenze falsch berechnet. Es müsste heißen
for m := 0 to degree - (n-1) do...

[rayman]

na klar! Dass ich das nich sehe, bzw wie ich auf son müll komme...

Sorum klappt das natürlich.
dickes danke.