1.2.2. Utilizzo dei coefficienti della formula del binomio di Newton per le derivate n-esime

// Calcolo la derivata n-esima di una funzione di classe C^(n-1)

// utilizzando come criterio di convergenza

// il confronto dei differenziali di sinistra con i differenziali di destra.

//

// 09/10/2019 saurischio@gmail.com

//

// PS: ancora da migliorare nel tener conto del caso di derivate successive nulle.

// In questo caso non può infatti essere utilizzato l'errore relativo.

// Per la funzione sotto indicata, il "bug" è per n=6,7,8,...

//---------------------------------------------------------------------

clc;clf;clear;

function f=funzione(x)

    f=4*x^5

endfunction

x0=2; // punto in cui calcolare la derivata

dx=0.1; // intervallo dx iniziale

errmin=0.001 // errore relativo richiesto

bisezioni=10 // numero massimo di bisezioni

n=6 // ordine della derivata

//

// La funzione "derivata" realizza una derivata n-esima e verifica la convergenza

//

// Gli INPUT sono:

// "x0" = il valore per cui calcolare la derivata

// "funct" = la funzione da derivare

// "dx" = il passo iniziale su cui calcolare il differenziale

// "minerrore" = l'errore relativo sotto il quale scendere

// "bisections" = il numero massimo di bisezioni consentito

// "n" = l'ordine della derivata

//

// Gli OUTPUT sono successioni di valori per passi sempre più piccoli:

// "f" = derivate

// "ddx" = passi dx

// "Ddx" = derivate a destra

// "Dsx" = derivate a sinistra

// "errore" = errori relativi

// "n" = l'ordine della derivata

//

function [f, ddx, Ddx, Dsx, errore, n]=derivata(x0, funct, dx, minerrore, bisections, n)

    function y=binomio(n, k)

        y=factorial(n)/(factorial(k)*factorial(n-k))

    endfunction

    i=0

    progressione=[]

    errore=minerrore+1

    while errore>minerrore && i<bisections

        i=i+1

        ddx(i)=2*(dx/n)/2^i

        Ddx(i)=0

        Dsx(i)=0

        for j=0:n

            Ddx(i)=Ddx(i)+binomio(n,j)*(-1)^j*funct(x0+(n-j)*ddx(i))/ddx(i)^n

            Dsx(i)=Dsx(i)+binomio(n,j)*(-1)^j*funct(x0-j*ddx(i))/ddx(i)^n

            errore(i)=abs((Ddx(i)-Dsx(i))/(Ddx(i)+Dsx(i)))

            f(i)=(Ddx(i)+Dsx(i))/2

        end 

        if i==bisections then

            disp("Stop at the "+string(bisections)+"-th bisection.")

        else end 

    end

endfunction

//

// Eseguo la funzione appena definita per i valori richiesti

//

[deriv,ddx,Ddx,Dsx,errore,n]=derivata(x0,funzione,dx,errmin,bisezioni,n)

disp("The derivative of order "+string(n)+" at point "+string(x0)+" is")

for k=1:size(deriv)(1)

    disp(string(deriv(k))+" +/- "+string(deriv(k).*errore(k)/2)+" for dx = "+string(ddx(k)))

end