prima implementazione di integratori

functions/integratore_adattivo_globale.m

@@ -0,0 +1,83 @@
+function [ area, err, ier ] = integratore_adattivo_globale( f, m, toll , maxite )
+%algoritmo dei trapezi composto, addattivo con controllo di errore globale
+%   input
+%   f = function handle
+%   m = matrice n x 2 ordinata in modo crescente per righe rispetto alla
+%   prima colonna. la prima colonna rappresenta un punto x, la seconda
+%   colonna la f(x). ogni valore della seconda colonna deve essere uguale a
+%   0 o essere la f(x) del punto x della prima colonna sulla medesima riga.
+%   toll = tolleranza assoluta voluta
+%   maxite = numero massimo di iterazioni
+%   output
+%   area = area finale
+%   err = errore assoluto
+%   ier = indice di errore
+    
+[x,~] = size(m);   
+i=1;
+k=i+1;
+j=x;
+ite = 0;
+ier = 0;
+report=[];
+    while (i <= j - 1) && (ite < maxite) 
+        [~,larea,lerr,m1] = static_quadratic(f,m(i:k,:));
+        m=[m(1:i-1,:);m1;m(k+1:j,:)]
+        ite = ite + 1;
+        j= j + 1;
+        report=[report;i,k+1,larea,lerr]
+        i = i + 2;
+        k = k + 2;
+    end
+    report = sortrows(report,4)
+    report = flipud (report)
+    [lenght,~] = size(report);
+    area=0;
+    err = 0;
+    for i=1:lenght
+      err= err + report(i,4);
+      area = area + report(i,3);
+    end
+    
+    while ( err > toll)&& (ite < maxite) 
+        i=report(1,1);
+        k=report(1,2) - 1;
+        [lenght,~] = size(report);
+        report=report(2:lenght,:)
+        [~,larea,lerr,m1] = static_quadratic(f,m(i:k,:));
+        m=[m(1:i-1,:);m1;m(k+1:j,:)]
+        ite = ite + 1;
+        j= j + 1;
+        report=[report;i,k+1,larea,lerr]
+        i=k+1;
+        k=k+2;
+        [~,larea,lerr,m1] = static_quadratic(f,m(i:k,:));
+        m=[m(1:i-1,:);m1;m(k+1:j,:)]
+        ite = ite + 1;
+        j= j + 1;
+        report=[report;i,k+1,larea,lerr];
+        report = sortrows(report,4)
+        report = flipud (report)
+        [lenght,~] = size(report);
+        area=0;
+        err = 0;
+        for i=1:lenght
+          err= err + report(i,4)
+          area = area + report(i,3)
+        end
+    end
+    
+    if ite >= maxite
+        disp('Warning: Massimo numero di step raggiunti')
+        ier = 1;
+
+    end
+    
+      
+    hold on
+    fplot (f,[m(1,1) m(j,1)],'k')
+    plot (m(:,1),m(:,2),'*r')
+    hold off
+    
+end
+