L'algorithme de la Méthode de Simpson Adaptative pour le calcul d'intégral en Matlab

L'algorithme de la Méthode de Simpson Adaptative pour le calcul d'intégral en Matlab :

function [JSf ,nodes]=simpadpt (fun ,a,b,tol ,hmin,varargin ) %SIMPADPT calcul numérique de l’ integrale avec la % méthode de Simpson adaptative. % JSF = SIMPADPT (FUN ,A,B,TOL ,HMIN) tente d’approcher % l’intégrale de la fonction FUN de A à B avec une % erreur inférieure à TOL en utilisant par récurrence % la méthode adaptative de Simpson avec H>=HMIN. % La fonction Y = FUN(X) doit accepter en % entrée un vecteur X et retourner dans un vecteur Y, % les valeurs de l’intégrande en chaque composante de V. % FUN peut être une fonction inline , une fonction % anonyme ou définie par un m-file. % JSF = SIMPADPT (FUN ,A,B,TOL ,HMIN,P1,P2 ,...) appelle la % fonction FUN en passant les paramètres optionnels % P1,P2 ,... de la maniere suivante : FUN(X,P1,P2 ,...). % [JSF ,NODES] = SIMPADPT (...) renvoie la % distribution des noeuds. A=[a,b]; N=[]; S=[]; JSf = 0; ba = 2*(b - a); nodes=[]; while ~isempty(A), [deltaI,ISc]=caldeltai(A,fun ,varargin {:}); if abs(deltaI) < 15*tol*(A(2)-A(1))/ba; JSf = JSf + ISc; S = union(S,A); nodes = [nodes, A(1) (A(1)+A(2))*0.5 A(2)]; S = [S(1), S(end)]; A = N; N = []; elseif A(2)-A(1) < hmin JSf=JSf+ISc; S = union(S,A); S = [S(1), S(end)]; A=N; N=[]; warning(’Pas d’’integration trop petit’); else Am = (A(1)+A(2))*0.5; A = [A(1) Am]; N = [Am, b]; end end nodes=unique(nodes); return function [deltaI ,ISc]=caldeltai(A,fun ,varargin) L=A(2)-A(1); t=[0; 0.25; 0.5; 0.75; 1]; x=L*t+A(1); L=L/6; w=[1; 4; 1]; wp=[1;4;2;4;1]; fx=feval(fun ,x, varargin {:}).* ones(5 ,1); IS=L*sum(fx([1 3 5]).*w); ISc=0.5*L*sum(fx.*wp); deltaI=IS-ISc; return