1.1.1. Realizzazione di programma in Scilab che integra numericamente e include criteri di convergenza

// Definisco una funzione che calcola l'integrale da x=a ad x=b

// seguendo un metodo che ricalca l'approccio di Riemann

//

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

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

clear;clc;clf;

function f=yy(x)  // funzione da integrare

    f=cos(x)

endfunction

//

// La funzione "integrale" realizza l'integrale

//

// Gli INPUT sono

// "funzione" = funzione da integrare

// "a" = minimo per la x

// "b" = massimo per la x

// "puntiiniziali" = # di punti utilizzati per il primo campionamento

// "minerrore" = errore relativo da raggiungere per ogni passo iniziale

// "bisezioni" = # massimo di bisezioni concesse per scendere sotto "minerrore"

//

// Gli OUTPUT sono:

//

// "f" = l'integrale calcolato

// "df" = l'incertezza associata

//

function [f, df]=integrale(funzione, a, b, puntiiniziali, minerrore, bisezioni)

    datax=a:(b-a)/puntiiniziali:b;

    datay=funzione(datax);

    dimx=size(datax)(2);

    f=0;

    df=0;

    for i=1:dimx-1

        sup=(datax(i+1)-datax(i))*max(datay(i),datay(i+1));

        inf=(datax(i+1)-datax(i))*min(datay(i),datay(i+1));

        j=0

        errore=0.5*(sup-inf)/(sup+inf)

        while (errore>minerrore && j<bisezioni) // operazioni di bisezione

            j=j+1

            xj=datax(i):(datax(i+1)-datax(i))/2^j:datax(i+1)

            yj=funzione(xj)

            supj=0;

            infj=0;

            for k=1:2^j

                supj=supj+(xj(k+1)-xj(k))*max(yj(k),yj(k+1));

                infj=infj+(xj(k+1)-xj(k))*min(yj(k),yj(k+1));

            end

            errore=0.5*(supj-infj)/(supj+infj);

            sup=supj;

            inf=infj;

            if (j==bisezioni) then printf("Stop alla "+string(bisezioni)+"-esima bisezione fra x="+string(datax(i))+" ed x="+string(datax(i+1))+"\n") end

        end

        f=f+0.5*(sup+inf)

        df=df+0.5*(sup-inf)

    end

    disp("L''integrale da "+string(a)+" a "+string(b)+" è "+string(f)+" piu'' o meno "+string(df))

    disp("")

endfunction

//

// Qui illustro il funzionamento testando per condizioni su errori relativi sempre più piccoli

[f,df]=integrale(yy,0,%pi/2,10,10^-1,12)

[f,df]=integrale(yy,0,%pi/2,10,10^-2,12)

[f,df]=integrale(yy,0,%pi/2,10,10^-3,12)

[f,df]=integrale(yy,0,%pi/2,10,10^-4,12)

[f,df]=integrale(yy,0,%pi/2,10,10^-7,12)