Ableitung einer Funktion numerisch bestimmen

Hallo, ich implementiere gerade das Newton-Verfahren zur Nullstellenbestimmung. Das funktioniert auch schon bestens, allerdings habe ich bisher in meinem Sourcecode die Ableitung einer Funktion f per Hand reingeschrieben.

Nun möchte ich das aber ändern und die Ableitung zur Laufzeit berechnen. Das Sekantenverfahren könnte zwar auch Nullstellen ohne die Ableitung bestimmen, aber dies nicht von Gleichungssystemen mit mehreren Unbekannten. Ich suche im Prinzip ein Verfahren um das Ergebnis einer Ableitung einer Funktion zur Laufzeit zu bestimmen.

double f(double x)
{
  return sin(x)*2*x*x;
}

double f_ableitung()
{
  return cos(x)*4*x; // diesen Schritt möchte ich dynamisch machen
}

double newton()
{
  //Implementierung von Newton
}

Die Lösung ist wahrscheinlich sehr einfach - ich komm aber einfach nicht drauf.

Bloops

Hm, das geht entweder

symbolisch (was jedoch sehr aufwendig werden dürfte -> CAS) oder
numerisch (z. B. über den Differentialquotienten...)

Bloops schrieb:

numerisch (z. B. über den Differentialquotienten...)

Wie würde das dann genau gehen? Würde das auch mit partiellen Ableitungen klappen? Hast Du Namen oder Links für die Verfahren?

Danke!

Walli

z.B. für zentrale Differenzen kann man eine (partielle) Ableitung so annähern: df(x,y,...)/dx = (f(x+h,y,...)-f(x-h,y,...))/(2*h). Das kannst du im Prinzip so runterprogrammieren. h wählst du einfach passend.

DrGreenthumb

btw, die ableitung von sin(x)*2*x*x ist aber nicht cos(x)*4*x

Taurin

Wenn du ganz naiv rangehst, geht folgenes: Im eindimensionalen Fall schnappst du dir einfach den Differenzenquotient: $f'(x) = \lim \frac{f(x) - f(x\_0)}{x - x\_0}$ . Weil du den Grenzwert numerisch aber nicht analytisch exakt ausrechnen willst, nimmst du einfach ein leicht verändertes $x_0 = x + \epsilon$ mit irgendeinem kleinen Wert für $\epsilon$ .

(Geometrisch legst du statt einer Tangente also nur eine Sekante an)

Im mehrdimensionalen Fall muss diesen "Trick" einfach auf jede Komponente des Gradienten/der Jakobi-Matrix.

Hab letztes Semester sowas mal in Matlab implementiert:

function J = jacobi(fun, x, varargin)
    h = 1e-7;

    y = feval(fun, x, varargin{:});

    DIM = length(x);

    J = zeros(length(y), DIM);
    for i = 1:DIM
        xDisturbed = x;
        if abs(x(i)) <= 1e-10 
          delta = h;
        else
          delta = h * x(i);
        xDisturbed(i) = x(i) + delta;  
        J(:,i) = (feval(fun, xDisturbed, varargin{:}) - y) / delta;
    end
end

b7f7

http://www.library.cornell.edu/nr/cbookcpdf.html
in der 5.7.