Remerciez-le!

Remerciez @Admin pour avoir partagé cet document gratuitement, de la manière la plus simple, en partageant sur les réseaux sociaux.

Quelques m´ ethodes num´ eriques en math´ ematique programm´ ees sous SCILAB. K

Quelques m´ ethodes num´ eriques en math´ ematique programm´ ees sous SCILAB. K.Barty 29 avril 2004 Table des mati` eres 1 R´ esolution de syst` emes lin´ eaires 1 1.1 ´ Elimination de Gauss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Factorisation LU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Factorisation de Cholesky . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.4 M´ ethode de Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2 R´ esolution num´ erique d’´ equation 7 2.1 M´ ethode de la Dichotomie . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2 M´ ethode de la Section Dor´ ee . . . . . . . . . . . . . . . . . . . . . . . . . . 9 3 Interpolation de Lagrange 11 4 Programmation Dynamique 12 4.1 Un probl` eme de gestion de stock . . . . . . . . . . . . . . . . . . . . . . . . 12 1 R´ esolution de syst` emes lin´ eaires Dans cette partie nous sommes concern´ es par le probl` eme suivant : R´ esoudre l’´ equation Ax = b ; (1) sans pour autant inverser la matrice A ce qui comme nous le savons est num´ eriquement tr` es coˆ uteux. 1 1.1 ´ Elimination de Gauss function [A,b]=algo(A,b,k) n=size(A,’c’) for i=k+1:n b(i)=b(i)-(A(i,k)/A(k,k))*b(k) for j=k+1:n A(i,j)=A(i,j)-(A(i,k)/A(k,k))*A(k,j) end; for j=1:k A(i,j)=0 end end; endfunction function [A,b]=gauss(A,b) // Methode de Gauss // On suppose que tous les pivots A(k,k) sont differents de 0. // En entree il y a les coefficients de la matrice A et du vecteur b qui // d´ eterminent le syst` eme lineaire Ax=b // En sortie le couple (A,b) determine un syt` eme lineaire equivalent, // mais numeriquement plus simple a resoudre, la matrice A etant // triangulaire superieure. n=size(A,’c’) for k=1:n [A,b]=algo(A,b,k) end; endfunction; A=[3 -4 0;12 -9 -1;24 -4 8] ←D´ eﬁnition de la A = ! 3. - 4. 0. ! ! 12. - 9. - 1. ! ! 24. - 4. 8. ! //matrice A b=[12;-6;23] ←Second membre du syst` eme b = ! 12. ! 2 ! - 6. ! ! 23. ! getf gauss.sci ←On charge la fonction gauss //dans scilab [A1,b1]=gauss(A,b) ←Appel de la fonction gauss b1 = ! 12. ! ! - 54. ! ! 143. ! A1 = ! 3. - 4. 0. ! ! 0. 7. - 1. ! ! 0. 0. 12. ! Le syst` eme : A1x = b1 ; est plus facile ` a r´ esoudre que (1) car la matrice A1 est triangulaire. 1.2 Factorisation LU On factorise A en un produit d’une matrice triangulaire inf´ erieure L et d’une matrice triangulaire sup´ erieure U, autrement dit A = LU. Cela permet une r´ esolution num´ erique plus facile en deux ´ etapes. La premi` ere ´ etape consiste ` a r´ esoudre le probl` eme : Lx′ = b ; la seconde ´ etape ` a r´ esoudre le probl` eme : Ux = x′. Le vecteur x est alors solution de l’´ equation (1). function [L,U]=facLU(A) // Factorisation LU d’une matrice inversible A // 3 U=zeros(A); L=eye(A); // U(1,:)=A(1,:); // dimmat=size(A,’r’); // for i=2:dimmat for j=1:i-1 L(i,j)=(A(i,j)-L(i,1:j-1)*U(1:j-1,j))/U(j,j); end for j=i:dimmat U(i,j)=A(i,j)-L(i,1:i-1)*U(1:i-1,j); end end endfunction A=[3 -4 0;12 -9 -1;24 -4 8] ←D´ eﬁnition de la A = ! 3. - 4. 0. ! ! 12. - 9. - 1. ! ! 24. - 4. 8. ! //matrice A getf facLU.sci ←On charge la fonction facLU //dans scilab [L,U]=facLU(A) ←Appel de la fonction facLU U = ! 3. - 4. 0. ! ! 0. 7. - 1. ! ! 0. 0. 12. ! L = ! 1. 0. 0. ! ! 4. 1. 0. ! ! 8. 4. 1. ! 4 1.3 Factorisation de Cholesky On d´ etermine la factorisation de Cholesky d’une matrice sym´ etrique A. Autrement dit on calcule la matrice C telle que A = CC′. function [C]=facCho(A) // Factorisation de Cholesky de la matrice A // C=zeros(A); C(1,1)=sqrt(A(1,1)); dimmat=size(A,’r’); for i=2:dimmat for j=1:i-1 C(i,j)=(A(i,j)-(C(i,1:j-1)*C(j,1:j-1)’))/C(j,j); end C(i,i)=sqrt(A(i,i)-(C(i,1:i-1)*C(i,1:i-1)’)); end endfunction A=[3 -4 0;12 -9 -1;24 -4 8] ←D´ eﬁnition de la A = ! 3. - 4. 0. ! ! 12. - 9. - 1. ! ! 24. - 4. 8. ! //matrice A b=[12;-6;23] ←Second membre du syst` eme b = ! 12. ! ! - 6. ! ! 23. ! getf gauss.sci ←On charge la fonction gauss //dans scilab [A1,b1]=gauss(A,b) ←Appel de la fonction gauss b1 = 5 ! 12. ! ! - 54. ! ! 143. ! A1 = ! 3. - 4. 0. ! ! 0. 7. - 1. ! ! 0. 0. 12. ! 1.4 M´ ethode de Jacobi La m´ ethode de Jacobi est une technique it´ erative pour r´ esoudre un syst` eme lin´ eaire. function [x]=Jacobi(A,b,x0,k) // // On veut resoudre le systeme Ax=b // de maniere recursive // x0 condition initiale de l’algorithme de Jacobi // k nombre d’it´ eration // n=size(A,’c’) x=x0; for p=1:k for i=1:n x(i)=(b(i)-A(i,1:n)*x+A(i,i)*x(i))/A(i,i) end end endfunction A=[264 110 98;110 105 21;98 21 54] A = ! 264. 110. 98. ! ! 110. 105. 21. ! ! 98. 21. 54. ! 6 b=[12;8;-5] b = ! 12. ! ! 8. ! ! - 5. ! x0=[0;0;0] x0 = ! 0. ! ! 0. ! ! 0. ! getf Jacobi.sci [x1]=Jacobi(A,b,x0,1000) ←Appel de x1 = ! 0.4176994 ! ! - 0.2074027 ! ! - 0.7699830 ! // la fonction //Jacobi norm(A*x1-b) ans = 3.202D-15 2 R´ esolution num´ erique d’´ equation 2.1 M´ ethode de la Dichotomie Nous allons pr´ esenter la m´ ethode de Dichotomie pour minimiser une fonction unimo- dale. function [intervalle,iter]=dicho(fonction,intervalle,epsilon,iter) // Methode de dichotomie utilisant la d´ ediv´ ee 7 // En entree // fonction : fonction supposee unimodale de la forme // y=fonction(x). // intervalle : vecteur 1x2, intervalle de depart pour l’optimisation. // epsilon : reel positif, precision de l’optimisation. // iter : entier, nombre d’iteration pour l’algotithme. // En sortie // intervalle : vecteur 1x2, resultat de la dichotomie. // iter : nombre d’iterations restantes a effectuer. // // Test de l’intervalle de d´ epart // init=iter; binf=min(intervalle); bsup=max(intervalle); if norm(binf-bsup)==%inf then error(’une des bornes de l’’intervalle est infinie’), end; bdemi = (binf+bsup)/2; f_binf = fonction(binf); f_bsup = fonction(bsup); f_bdemi = fonction(bdemi); u=file(’open’,’results_dicho’,’unknown’); // Methode de dichotomie while ((abs(bsup-binf)>= epsilon) & (iter >= 1)) b1quart = (bdemi+binf)/2; b3quart = (bsup+bdemi)/2; f_b1quart = fonction(b1quart); f_b3quart = fonction(b3quart); tab=[f_binf,f_b1quart,f_bdemi,f_b3quart,f_bsup]; [a,b]=gsort(tab,’c’,’i’); if (b(1) == 1) | (b(1) == 2) then bsup = bdemi; bdemi = b1quart; f_bsup = f_bdemi; f_bdemi = f_b1quart; elseif b(1) == 3 then binf = b1quart; bsup = b3quart; f_binf = f_b1quart; f_bsup = f_b3quart; elseif (b(1) == 4) | (b(1) == 5) then binf = bdemi; 8 bdemi = b3quart; f_binf = f_bdemi; f_bdemi = f_b3quart; end fprintf(u,’iter = %2.0f binf = %10.8f bsup = %10.8f precision %10.8f’,... init-iter+1,binf,bsup,abs(bsup-binf)); iter=iter-1; end intervalle=[binf,bsup]; file(’close’,u) endfunction deff("[y]=func(x)","y=x^2-4*x+6") ←D´ eﬁnition de la //fonction func getf dicho.sci ←On charge la fonction dicho //dans scilab [intervalle,iter]=dicho(func,[-20,10],0.001,10) iter = 0. intervalle = ! 1.9873047 2.0166016 ! 2.2 M´ ethode de la Section Dor´ ee C’est une variante de la m´ ethode de Dichotomie, elle permet entre autre de n’avoir ` a ´ evaluer qu’une seule fois la fonction ` a minimiser ` a chaque it´ eration. function [intervalle,iter]=section_doree(fonction,intervalle,epsilon,iter) // init=iter; binf=min(intervalle); bsup=max(intervalle); tau=(1+sqrt(5))/2; tau=1/tau; if norm(binf-bsup)==%inf then error(’une des bornes de l’’intervalle est infinie’), end; b1 = bsup-tau*(bsup-binf); 9 b2 = binf+tau*(bsup-binf); f_binf = fonction(binf); f_bsup = fonction(bsup); f_b1 = fonction(b1); f_b2 = fonction(b2); u=file(’open’,’results_doree’,’unknown’); // Methode de la section doree while ((abs(binf-bsup)>= epsilon) & (iter >= 1)) tab=[f_b1,f_b2]; [a,b]=gsort(tab,’c’,’i’); if b(1)==1 then bsup=b2; b2=b1; f_b2=f_b1; b1=binf+bsup-b2; f_b1=fonction(b1); else binf=b1; b1=b2; f_b1=f_b2; b2=binf+bsup-b1; f_b2=fonction(b2); end fprintf(u,’iter = %2.0f binf = %10.8f bsup = %10.8f precision %10.8f’,... init-iter+1,binf,bsup,abs(bsup-binf)); iter=iter-1; end file(’close’,u) intervalle=[binf,bsup]; endfunction deff("[y]=func(x)","y=x^2-4*x+6") ←D´ eﬁnition de la //fonction func getf section_doree.sci ←On charge la fonction dicho //dans scilab [intervalle,iter]=section_doree(func,[-20,10],0.001,10) iter = 10 0. intervalle = ! 1.8847051 2.1286236 ! 3 Interpolation de Lagrange ´ Etant donn´ e une fonction f : R →R et x0, . . . , xn n + 1 points dans R. Lagrange ` a r´ esolu le probl` eme de d´ eterminer uploads/Litterature/ tp-scilab.pdf