Algorithme de la Méthode de Hamming pour résoudre les équations différentielles vectorielles en Matlab

Algorithme de la Méthode de Hamming pour résoudre les équations différentielles vectorielles en Matlab :

function [t,y] = ode_Ham(f,tspan,y0,N,KC,varargin) % Hamming method to solve vector d.e. y’(t) = f(t,y(t)) % for tspan = [t0,tf] and with the initial value y0 and N time steps % using the modifier based on the error estimate depending on KC = 1/0 if nargin < 5, KC = 1; end %with modifier by default if nargin < 4 | N <= 0, N = 100; end %default maximum number of iterations if nargin < 3, y0 = 0; end %default initial value y0 = y0(:)’; end %make it a row vector h = (tspan(2)-tspan(1))/N; %step size tspan0 = tspan(1)+[0 3]*h; [t,y] = ode_RK4(f,tspan0,y0,3,varargin{:}); %Initialize by Runge-Kutta t = [t(1:3)’ t(4):h:tspan(2)]’; for k = 2:4, F(k - 1,:) = feval(f,t(k),y(k,:),varargin{:}); end p = y(4,:); c = y(4,:); h34 = h/3*4; KC11 = KC*112/121; KC91 = KC*9/121; h312 = 3*h*[-1 2 1]; for k = 4:N p1 = y(k - 3,:) + h34*(2*(F(1,:) + F(3,:)) - F(2,:)); %Eq.(6.4.9a) m1 = p1 + KC11*(c - p); %Eq.(6.4.9b) c1 = (-y(k - 2,:) + 9*y(k,:) +... h312*[F(2:3,:); feval(f,t(k + 1),m1,varargin{:})])/8; %Eq.(6.4.9c) y(k+1,:) = c1 - KC91*(c1 - p1); %Eq.(6.4.9d) p = p1; c = c1; %update the predicted/corrected values F = [F(2:3,:); feval(f,t(k + 1),y(k + 1,:),varargin{:})]; end

Algorithme de la Méthode d'Adams-Bashforth-Moulton pour résoudre les équations différentielles vectorielles en Matlab

Algorithme de la Méthode d'Adams-Bashforth-Moulton pour résoudre les équations différentielles vectorielles en Matlab :

function [t,y] = ode_ABM(f,tspan,y0,N,KC,varargin) %Adams-Bashforth-Moulton method to solve vector d.e. y’(t) = f(t,y(t)) % for tspan = [t0,tf] and with the initial value y0 and N time steps % using the modifier based on the error estimate depending on KC = 1/0 if nargin < 5, KC = 1; end %with modifier by default if nargin < 4 | N < = 0, N = 100; end %default maximum number of iterations y0 = y0(:)’; %make it a row vector h = (tspan(2) - tspan(1))/N; %step size tspan0 = tspan(1)+[0 3]*h; [t,y] = rk4(f,tspan0,y0,3,varargin{:}); %initialize by Runge-Kutta t = [t(1:3)’ t(4):h:tspan(2)]’; for k = 1:4, F(k,:) = feval(f,t(k),y(k,:),varargin{:}); end p = y(4,:); c = y(4,:); KC22 = KC*251/270; KC12 = KC*19/270; h24 = h/24; h241 = h24*[1 -5 19 9]; h249 = h24*[-9 37 -59 55]; for k = 4:N p1 = y(k,:) +h249*F; %Eq.(6.4.8a) m1 = pk1 + KC22*(c-p); %Eq.(6.4.8b) c1 = y(k,:)+ ... h241*[F(2:4,:); feval(f,t(k + 1),m1,varargin{:})]; %Eq.(6.4.8c) y(k + 1,:) = c1 - KC12*(c1 - p1); %Eq.(6.4.8d) p = p1; c = c1; %update the predicted/corrected values F = [F(2:4,:); feval(f,t(k + 1),y(k + 1,:),varargin{:})]; End

Algorithme de la Méthode de Runge-Kutta pour résoudre les équations différentielles vectorielles en Matlab

Algorithme de la Méthode de Runge-Kutta pour résoudre les équations différentielles vectorielles en Matlab :

function [t,y] = ode_RK4(f,tspan,y0,N,varargin) %Runge-Kutta method to solve vector differential eqn y’(t) = f(t,y(t)) % for tspan = [t0,tf] and with the initial value y0 and N time steps if nargin < 4 | N <= 0, N = 100; end if nargin < 3, y0 = 0; end y(1,:) = y0(:)’; %make it a row vector h = (tspan(2) - tspan(1))/N; t = tspan(1)+[0:N]’*h; for k = 1:N f1 = h*feval(f,t(k),y(k,:),varargin{:}); f1 = f1(:)’; %(6.3.2a) f2 = h*feval(f,t(k) + h/2,y(k,:) + f1/2,varargin{:}); f2 = f2(:)’;%(6.3.2b) f3 = h*feval(f,t(k) + h/2,y(k,:) + f2/2,varargin{:}); f3 = f3(:)’;%(6.3.2c) f4 = h*feval(f,t(k) + h,y(k,:) + f3,varargin{:}); f4 = f4(:)’; %(6.3.2d) y(k + 1,:) = y(k,:) + (f1 + 2*(f2 + f3) + f4)/6; %Eq.(6.3.1) end

Algorithme de la Méthode de Heun pour résoudre les équations différentielles vectorielles en Matlab

Algorithme de la Méthode de Heun pour résoudre les équations différentielles vectorielles en Matlab :

function [t,y] = ode_Heun(f,tspan,y0,N) %Heun method to solve vector differential equation y’(t) = f(t,y(t)) % for tspan = [t0,tf] and with the initial value y0 and N time steps if nargin<4 | N <= 0, N = 100; end if nargin<3, y0 = 0; end h = (tspan(2) - tspan(1))/N; %stepsize t = tspan(1)+[0:N]’*h; %time vector y(1,:) = y0(:)’; %always make the initial value a row vector for k = 1:N fk = feval(f,t(k),y(k,:)); y(k+1,:) = y(k,:)+h*fk; %Eq.(6.2.3) y(k+1,:) = y(k,:) +h/2*(fk +feval(f,t(k+1),y(k+1,:))); %Eq.(6.2.4) end

Algorithme de la Méthode d'Euler pour résoudre les équations différentielles de premier ordre en Matlab

Algorithme de la Méthode d'Euler pour résoudre les équations différentielles de premier ordre en Matlab :

%nm610: Euler method to solve a 1st-order differential equation clear, clf a = 1; r = 1; y0 = 0; tf = 2; t = [0:0.01:tf]; yt=1- exp(-a*t); %Eq.(6.1.5): true analytical solution plot(t,yt,’k’), hold on klasts = [8 4 2]; hs = tf./klasts; y(1) = y0; for itr = 1:3 %with various step size h = 1/8,1/4,1/2 klast = klasts(itr); h = hs(itr); y(1)=y0; for k = 1:klast y(k + 1) = (1 - a*h)*y(k) +h*r; %Eq.(6.1.3): plot([k - 1 k]*h,[y(k) y(k+1)],’b’, k*h,y(k+1),’ro’) if k < 4, pause; end end end

Algorithme de la Méthode de la Sécante en Matlab

Algorithme de la Méthode de la Sécante en Matlab :

function [x,fx,xx] = secant(f,x0,TolX,MaxIter,varargin) % solve f(x) = 0 by using the secant method. %input : f = ftn to be given as a string ’f’ if defined in an M-file % x0 = the initial guess of the solution % TolX = the upper limit of |x(k) - x(k - 1)| % MaxIter = the maximum # of iteration %output: x = the point which the algorithm has reached % fx = f(x(last)), xx = the history of x h = 1e-4; h2 = 2*h; TolFun=eps; xx(1) = x0; fx = feval(f,x0,varargin{:}); for k = 1: MaxIter if k <= 1, dfdx = (feval(f,xx(k) + h,varargin{:})-... feval(f,xx(k) - h,varargin{:}))/h2; else dfdx = (fx - fx0)/dx; end dx = -fx/dfdx; xx(k + 1) = xx(k) + dx; %Eq.(4.5.2) fx0 = fx; fx = feval(f,xx(k+1)); if abs(fx) < TolFun | abs(dx) < TolX, break; end end x = xx(k + 1); if k == MaxIter, fprintf(’The best in %d iterations\n’,MaxIter), end

Algorithme de la Méthode de la Position Fausse en Matlab

Algorithme de la Méthode de la Position Fausse en Matlab :

function [x,err,xx] = falsp(f,a,b,TolX,MaxIter) %bisct.m to solve f(x)=0 by using the false position method. %input : f = ftn to be given as a string ’f’ if defined in an M-file % a/b = initial left/right point of the solution interval % TolX = upperbound of error(max(|x(k)-a|,|b-x(k)|)) % MaxIter = maximum # of iterations %output: x = point which the algorithm has reached % err = max(x(last)-a|,|b-x(last)|) % xx = history of x TolFun = eps; fa = feval(f,a); fb=feval(f,b); if fa*fb > 0, error(’We must have f(a)f(b)<0!’); end for k = 1: MaxIter xx(k) = (a*fb-b*fa)/(fb-fa); %Eq.(4.3.1) fx = feval(f,xx(k)); err = max(abs(xx(k) - a),abs(b - xx(k))); if abs(fx) < TolFun | err<TolX, break; elseif fx*fa > 0, a = xx(k); fa = fx; else b = xx(k); fb = fx; end end x = xx(k); if k == MaxIter, fprintf(’The best in %d iterations\n’,MaxIter), end

Algorithme d'Itération de Gauss–Seidel en Matlab

Algorithme d'Itération de Gauss–Seidel en Matlab :

function X = gauseid(A,B,X0,kmax) %This function finds x = A^-1 B by Gauss–Seidel iteration. if nargin < 4, tol = 1e-6; kmax = 100; elseif kmax < 1, tol = max(kmax,1e-16); kmax = 1000; else tol = 1e-6; end if nargin < 4, tol = 1e-6; kmax = 100; end if nargin < 3, X0 = zeros(size(B)); end NA = size(A,1); X = X0; for k = 1: kmax X(1,:) = (B(1,:)-A(1,2:NA)*X(2:NA,:))/A(1,1); for m = 2:NA-1 tmp = B(m,:)-A(m,1:m-1)*X(1:m - 1,:)-A(m,m + 1:NA)*X(m + 1:NA,:); X(m,:) = tmp/A(m,m); %Eq.(2.5.4) end X(NA,:) = (B(NA,:)-A(NA,1:NA - 1)*X(1:NA - 1,:))/A(NA,NA); if nargout == 0, X, end %To see the intermediate results if norm(X - X0)/(norm(X0) + eps)<tol, break; end X0 = X; end

Algorithme d'Itérations de Jacobi en Matlab

Algorithme d'Itérations de Jacobi en Matlab :

function X = jacobi(A,B,X0,kmax) %This function finds a soltuion to Ax = B by Jacobi iteration. if nargin < 4, tol = 1e-6; kmax = 100; %called by jacobi(A,B,X0) elseif kmax < 1, tol = max(kmax,1e-16); kmax = 100; %jacobi(A,B,X0,tol) else tol = 1e-6; %jacobi(A,B,X0,kmax) end if nargin < 3, X0 = zeros(size(B)); end NA = size(A,1); X = X0; At = zeros(NA,NA); for m = 1:NA for n = 1:NA if n ~= m, At(m,n) = -A(m,n)/A(m,m); end end Bt(m,:) = B(m,:)/A(m,m); end for k = 1: kmax X = At*X + Bt; %Eq. (2.5.3) if nargout == 0, X, end %To see the intermediate results if norm(X - X0)/(norm(X0) + eps) < tol, break; end X0 = X; end

Algorithme de la Factorisation LU en Matlab

Algorithme de la Factorisation LU en Matlab :

function [L,U,P] = lu_dcmp(A) %This gives LU decomposition of A with the permutation matrix P % denoting the row switch(exchange) during factorization NA = size(A,1); AP = [A eye(NA)]; %augment with the permutation matrix. for k = 1:NA - 1 %Partial Pivoting at AP(k,k) [akx, kx] = max(abs(AP(k:NA,k))); if akx < eps error(’Singular matrix and No LU decomposition’) end mx = k+kx-1; if kx > 1 % Row change if necessary tmp_row = AP(k,:); AP(k,:) = AP(mx,:); AP(mx,:) = tmp_row; end % LU decomposition for m = k + 1: NA AP(m,k) = AP(m,k)/AP(k,k); %Eq.(2.4.8.2) AP(m,k+1:NA) = AP(m,k + 1:NA)-AP(m,k)*AP(k,k + 1:NA); %Eq.(2.4.9) end end P = AP(1:NA, NA + 1:NA + NA); %Permutation matrix for m = 1:NA for n = 1:NA if m == n, L(m,m) = 1.; U(m,m) = AP(m,m); elseif m > n, L(m,n) = AP(m,n); U(m,n) = 0.; else L(m,n) = 0.; U(m,n) = AP(m,n); end end end if nargout == 0, disp(’L*U = P*A with’); L,U,P, end %You can check if P’*L*U = A?

Algorithme d’Élimination de Gauss en Matlab

Algorithme d’Élimination de Gauss en Matlab :

function x = gauss(A,B) %The sizes of matrices A,B are supposed to be NA x NA and NA x NB. %This function solves Ax = B by Gauss elimination algorithm. NA = size(A,2); [NB1,NB] = size(B); if NB1 ~= NA, error(’A and B must have compatible dimensions’); end N = NA + NB; AB = [A(1:NA,1:NA) B(1:NA,1:NB)]; % Augmented matrix epss = eps*ones(NA,1); for k = 1:NA %Scaled Partial Pivoting at AB(k,k) by Eq.(2.2.20) [akx,kx] = max(abs(AB(k:NA,k))./ ... max(abs([AB(k:NA,k + 1:NA) epss(1:NA - k + 1)]’))’); if akx < eps, error(’Singular matrix and No unique solution’); end mx = k + kx - 1; if kx > 1 % Row change if necessary tmp_row = AB(k,k:N); AB(k,k:N) = AB(mx,k:N); AB(mx,k:N) = tmp_row; end % Gauss forward elimination AB(k,k + 1:N) = AB(k,k+1:N)/AB(k,k); AB(k,k) = 1; %make each diagonal element one for m = k + 1: NA AB(m,k+1:N) = AB(m,k+1:N) - AB(m,k)*AB(k,k+1:N); %Eq.(2.2.5) AB(m,k) = 0; end end %backward substitution for a upper-triangular matrix eqation % having all the diagonal elements equal to one x(NA,:) = AB(NA,NA+1:N); for m = NA-1: -1:1 x(m,:) = AB(m,NA + 1:N)-AB(m,m + 1:NA)*x(m + 1:NA,:); %Eq.(2.2.7) end

Algorithme de la Méthode de Newmark pour la résolution d'équations différentielles en Matlab

Algorithme de la Méthode de Newmark pour la résolution d'équations différentielles en Matlab :

function [t,u]=newmark (odefun,tspan,y0,Nh,param ,... varargin) %NEWMARK résout une équation différentielle du second % ordre avec la méthode de Newmark % [T,Y]=NEWMARK (ODEFUN,TSPAN,Y0,NH,PARAM) avec TSPAN = % [T0 TF] intègre le système d’ équations différen - % tielles y’’=f(t,y,y’) du temps T0 au temps TF avec % la condition initiale Y0=(y(t0),y’(t0) en utilisant % la méthode de Newmark sur une grille de NH % intervalles équidistribués. % PARAM contient les paramètres zeta et theta. % La fonction ODEFUN(T,Y) doit retourner un vecteur % contenant les évaluations de f(t,y) et de même % dimension que Y. Chaque ligne de la solution Y % correspond à un temps contenu dans le vecteur % colonne T. tt=linspace (tspan(1),tspan(2),Nh+1); y=y0(:); u=y.’; global glob_h glob_t glob_y glob_odefun; global glob_zeta glob_theta glob_varargin glob_fn; glob_h =(tspan(2)-tspan (1))/Nh; glob_y=y; glob_odefun=odefun; glob_zeta = param(1); glob_theta = param(2); glob_varargin=varargin ; if ( exist(’OCTAVE_VERSION’) ) o_ver=OCTAVE_VERSION; version =str2num ([o_ver(1),o_ver(3),o_ver (5)]); end if ( ~exist( ’OCTAVE_VERSION’ ) | version >= 320 ) options=optimset ; options.Display =’off’; options.TolFun =1.e-12; options.MaxFunEvals=10000; end glob_fn =feval(odefun,tt(1),glob_y,varargin {:}); for glob_t=tt(2:end) if ( exist( ’OCTAVE_VERSION’ ) & version < 320 ) w = fsolve(’newmarkfun’, glob_y ); else w = fsolve(@(w) newmarkfun(w),glob_y,options ); end glob_fn =feval(odefun,glob_t,w,varargin {:}); u = [u; w.’]; glob_y = w; end t=tt; clear glob_h glob_t glob_y glob_odefun; clear glob_zeta glob_theta glob_varargin glob_fn ; end function z=newmarkfun(w) global glob_h glob_t glob_y glob_odefun; global glob_zeta glob_theta glob_varargin glob_fn ; fn1=feval(glob_odefun ,glob_t,w, glob_varargin{:}); z(1)=w(1) - glob_y (1) -glob_h*glob_y (2) -... glob_h ^2*( glob_zeta*fn1+(0.5 -glob_zeta)*glob_fn ); z(2)=w(2) - glob_y (2) -... glob_h *((1 -glob_theta)*glob_fn +glob_theta*fn1); end

Algorithme de la Méthode de Prédicteur-Correcteur pour la résolution d'équations différentielles en Matlab

Algorithme de la Méthode de Prédicteur-Correcteur pour la résolution d'équations différentielles en Matlab :

function [t,u]=predcor (odefun,tspan,y0,Nh ,... predictor ,corrector ,varargin ) %PREDCOR Résout une équation différentielle avec une % méthode predicteur -correcteur % [T,Y]=PREDCOR (ODEFUN,TSPAN,Y0,NH,PRED,CORR) avec % TSPAN=[T0 TF] % intègre le système d’équations différentielles % y’=f(t,y) du temps T0 au temps TF avec la condition % initiale Y0 en utilisant une méthode générale % prédicteur -correcteur sur une grille de NH % intervalles équidistribués. La fonction ODEFUN(T,Y) % doit retourner un vecteur correspondant à f(t,y) % de même dimension que Y. % Chaque ligne de la solution Y correspond % à un temps du vecteur colonne T. % [T,Y]=PREDCOR (ODEFUN,TSPAN,Y0,NH,PRED,CORR,P1 ,..) % passe les paramètres supplémentaires P1,P2 ,.. aux % fonctions ODEFUN, PRED et COOR de la manière % CORR(T,Y,P1,P2...). h=(tspan(2)-tspan(1))/Nh; y=y0(:); w=y; u=y.’; tt=linspace (tspan(1),tspan(2),Nh+1); for t=tt(1:end -1) fn = feval(odefun ,t,w,varargin {:}); upre = feval(predictor ,t,w,h,fn); w = feval(corrector ,t+h,w,upre,h,odefun ,... fn,varargin {:}); u = [u; w.’]; end t = tt; end

Algorithme de la Méthode de Crank-Nicolson pour la résolution d'équations différentielles en Matlab

Algorithme de la Méthode de Crank-Nicolson pour la résolution d'équations différentielles en Matlab :

function [t,u]=cranknic(odefun,tspan,y0,Nh,varargin ) %CRANKNIC Résout une équation différentielle avec la % méthode de Crank-Nicolson . % [T,Y]=CRANKNIC (ODEFUN,TSPAN,Y0 ,NH) avec % TSPAN=[T0,TF] % intègre le système d’équations différentielles % y’=f(t,y) du temps T0 au temps TF avec la condition % initiale Y0 en utilisant la méthode de % Crank-Nicolson sur une grille de NH intervalles % équidistribués. La fonction ODEFUN(T,Y) doit % retourner un vecteur correspondant à f(t,y) % de même dimension que Y. % Chaque ligne de la solution Y correspond % à un temps du vecteur colonne T. % [T,Y] = CRANKNIC (ODEFUN,TSPAN,Y0,NH,P1,P2 ,...) % passe les paramètres supplémentaires P1,P2 ,.. à % la fonction ODEFUN de la manière suivante: % ODEFUN(T,Y,P1 ,P2...). tt=linspace (tspan(1),tspan(2),Nh+1); y=y0(:); % crée toujours un vecteur colonne u=y.’; global glob_h glob_t glob_y glob_odefun; glob_h =(tspan(2)-tspan (1))/Nh; glob_y=y; glob_odefun=odefun; if ( exist(’OCTAVE_VERSION’) ) o_ver=OCTAVE_VERSION; version =str2num ([o_ver(1),o_ver(3),o_ver (5)]); end if( ~exist(’OCTAVE_VERSION’) | version >= 320 ) options=optimset ; options.Display =’off’; options.TolFun =1.e-12; options.MaxFunEvals=10000; end for glob_t=tt(2:end) if ( exist(’OCTAVE_VERSION’) & version < 320 ) w = fsolve(’cranknicfun’,glob_y); else w = fsolve(@(w) cranknicfun(w),glob_y,options ); end u = [u; w.’]; glob_y = w; end t=tt; clear glob_h glob_t glob_y glob_odefun; end function z= cranknicfun(w) global glob_h glob_t glob_y glob_odefun; z=w - glob_y - ... 0.5* glob_h*(feval(glob_odefun ,glob_t,w) + ... feval(glob_odefun ,glob_t-glob_h ,glob_y)); end

Algorithme de la Méthode d’Euler Implicite pour la résolution d'équation différentielles en Matlab

Algorithme de la Méthode d’Euler Implicite pour la résolution d'équation différentielles en Matlab :

function [t,u]=beuler(odefun,tspan,y0,Nh,varargin ) %BEULER Résout une équation différentielle avec la % méthode d’Euler implicite. % [T,Y]=BEULER(ODEFUN,TSPAN,Y0,NH) avec TSPAN=[T0,TF] % intègre le système d’équations différentielles % y’=f(t,y) du temps T0 au temps TF avec la condition % initiale Y0 en utilisant la méthode d’Euler % implicite sur une grille de NH intervalles % équidistribués. La fonction ODEFUN(T,Y) doit % retourner un vecteur, correspondant à f(t,y), % de même dimension que Y. % Chaque ligne de la solution Y correspond % à un temps du vecteur colonne T. % [T,Y] = BEULER(ODEFUN,TSPAN,Y0 ,NH,P1,P2 ,...) passe % les paramètres supplémentaires P1,P2 ,.. à la % fonction ODEFUN de la manière suivante : % ODEFUN(T,Y,P1 ,P2...). tt=linspace (tspan(1),tspan(2),Nh+1); y=y0(:); % crée toujours un vecteur colonne u=y.’; global glob_h glob_t glob_y glob_odefun; glob_h =(tspan(2)-tspan (1))/Nh; glob_y=y; glob_odefun=odefun; glob_t=tt(2); if ( exist(’OCTAVE_VERSION’) ) o_ver=OCTAVE_VERSION; version =str2num ([o_ver(1),o_ver(3),o_ver (5)]); end if ( ~exist(’OCTAVE_VERSION’) | version >= 320 ) options =optimset ; options .Display=’off’; options .TolFun=1.e -12; options .MaxFunEvals=10000; end for glob_t=tt(2:end) if ( exist(’OCTAVE_VERSION’) & version < 320 ) w = fsolve(’beulerfun’,glob_y); else w = fsolve(@(w) beulerfun(w),glob_y ,options ); end u = [u; w.’]; glob_y = w; end t=tt; clear glob_h glob_t glob_y glob_odefun; end function [z]=beulerfun(w) global glob_h glob_t glob_y glob_odefun; z=w-glob_y-glob_h*feval(glob_odefun ,glob_t,w); end

Algorithme de la Méthode d’Euler Explicite pour la résolution d'équation différentielles en Matlab

Algorithme de la Méthode d’Euler Explicite pour la résolution d'équation différentielles en Matlab :

function [t,u]=feuler(odefun,tspan,y0,Nh,varargin ) %FEULER Résout une équation différentielle avec la % méthode d’Euler explicite. % [T,Y]=FEULER(ODEFUN,TSPAN,Y0,NH) avec TSPAN=[T0,TF] % intègre le système d’équations différentielles % y’=f(t,y) du temps T0 au temps TF avec la condition % initiale Y0 en utilisant la méthode d’Euler % explicite sur une grille de NH intervalles % équidistribués. La fonction ODEFUN(T,Y) doit % retourner un vecteur, correspondant à f(t,y), % de même dimension que Y. % Chaque ligne de la solution Y correspond % à un temps du vecteur colonne T. % [T,Y] = FEULER(ODEFUN,TSPAN,Y0 ,NH,P1,P2 ,...) passe % les paramètres supplémentaires P1,P2 ,.. à la % fonction ODEFUN de la maniere suivante: % ODEFUN(T,Y,P1 ,P2...). h=(tspan(2)-tspan(1))/Nh; y=y0(:); % crée toujours un vecteur colonne w=y; u=y.’; tt=linspace (tspan(1),tspan(2),Nh+1); for t = tt(1:end -1) w=w+h*feval(odefun,t,w,varargin {:}); u = [u; w.’]; end t=tt; return

Algorithme de la Méthode des Itérations QR pour le calcul des valeurs propres d'une matrice en Matlab

Algorithme de la Méthode des Itérations QR pour le calcul des valeurs propres d'une matrice en Matlab :

function D=qrbasic (A,tol ,nmax) %QRBASIC calcule les valeurs propres de la matrice A. % D=QRBASIC(A,TOL ,NMAX) calcule par itérations QR % toutes les valeurs propres de A avec une tolérance % TOL en NMAX itérations au maximum. La convergence % de cette méthode n’est pas toujours garantie . [n,m]=size(A); if n ~= m, error(’Matrices carrées seulement’); end T = A; niter = 0; test = norm(tril(A,-1),inf); while niter <= nmax & test >= tol [Q,R]=qr(T); T = R*Q; niter = niter + 1; test = norm(tril(T,-1),inf); end if niter > nmax warning([’La méthode ne converge pas dans le ’ ’nombre d’’itérations maximum voulu\n’]); else fprintf([’La methode converge en ’ ... ’%i itérations\n’],niter); end D = diag(T); return

Algorithme des Disques de Gershgorin en Matlab

Algorithme des Disques de Gershgorin en Matlab :

function gershcircles(A) %GERSHCIRCLES trace les disques de Gershgorin % GERSHCIRCLES(A) trace les disques de Gershgorin % pour la matrice carrée A et sa transposée. n = size(A); if n(1) ~= n(2) error(’Matrices carrées seulement’); else n=n(1); circler=zeros(n ,201); circlec=circler; end center = diag(A); radiic = sum(abs(A-diag(center ))); radiir = sum(abs(A’-diag(center ))); one = ones (1,201); cosisin = exp(i*[0:pi/100:2* pi]); figure (1); title(’Disques des lignes ’); xlabel(’Re’); ylabel(’Im’); figure (2); title(’Disques des colonnes ’); xlabel(’Re’); ylabel(’Im’); for k = 1:n circlec(k,:) = center(k)*one + radiic(k)*cosisin ; circler(k,:) = center(k)*one + radiir(k)*cosisin ; figure (1); patch(real(circler (k,:)) ,imag(circler (k,:)) ,’red’); hold on plot(real(circler(k,:)) ,imag(circler(k,:)) ,’k-’ ,... real(center(k)),imag(center(k)),’kx’); figure (2); patch(real(circlec (k,:)) ,imag(circlec (k,:)) ,’green’); hold on plot(real(circlec(k,:)) ,imag(circlec(k,:)) ,’k-’ ,... real(center(k)),imag(center(k)),’kx’); end for k = 1:n figure (1); plot(real(circler(k,:)) ,imag(circler(k,:)) ,’k-’ ,... real(center(k)),imag(center(k)),’kx’); figure (2); plot(real(circlec(k,:)) ,imag(circlec(k,:)) ,’k-’ ,... real(center(k)),imag(center(k)),’kx’); end figure (1); axis image; hold off; figure (2); axis image; hold off return

Algorithme de la Méthode de la Puissance Inverse avec Décalage pour la recherche de valeurs propres en Matalab

Algorithme de la Méthode de la Puissance Inverse avec Décalage pour la recherche de valeurs propres en Matalab :

function [lambda ,x,iter]=invshift (A,mu,tol ,nmax,x0) %INVSHIFT Evalue numériquement une valeur propre d’une % matrice % LAMBDA=INVSHIFT (A) calcule la valeur propre de A % de module minimum avec la méthode de la puissance % inverse % LAMBDA=INVSHIFT (A,MU) calcule la valeur propre de A % la plus proche d’un nombre réel ou complexe MU % LAMBDA=INVSHIFT (A,MU,TOL ,NMAX,X0) utilise une % tolérance sur l’erreur absolue TOL (1e-6 par défaut) % un nombre maximum d’itérations NMAX (100 par défaut) % en démarrant d’un vecteur initial X0. % [LAMBDA,V,ITER]=INVSHIFT(A,MU,TOL ,NMAX,X0) retourne % aussi un vecteur propre V tel que A*V=LAMBDA*V et % l’itération à laquelle V est calculée. [n,m]=size(A); if n ~= m, error(’Matrices carrées seulement’); end if nargin == 1 x0 = rand(n,1); nmax = 100; tol = 1.e-06; mu = 0; elseif nargin == 2 x0 = rand(n,1); nmax = 100; tol = 1.e-06; end [L,U]=lu(A-mu*eye(n)); if norm(x0) == 0 x0 = rand(n ,1); end x0=x0/norm(x0); z0=L\x0; pro=U\z0; lambda=x0’*pro; err=tol*abs(lambda )+1; iter=0; while err >tol*abs(lambda)&abs(lambda )~=0&iter <=nmax x = pro; x = x/norm(x); z=L\x; pro=U\z; lambdanew = x’*pro; err = abs( lambdanew - lambda); lambda = lambdanew; iter = iter + 1; end lambda = 1/lambda + mu; return

Algorithme de la Méthode de Puissance pour le calcul des Valeurs Propres d'une Matrice en Matlab

Algorithme de la Méthode de Puissance pour le calcul des Valeurs Propres d'une Matrice en Matlab :

function [lambda ,x,iter]=eigpower (A,tol ,nmax,x0) %EIGPOWER Evalue numériquement une valeur propre % d’une matrice % LAMBDA=EIGPOWER (A) calcule avec la méthode de la % puissance la valeur propre de A de module maximal % à partir d’une donnée initial qui par défaut est % le vecteur constitué de 1 % LAMBDA=EIGPOWER (A,TOL ,NMAX,X0) utilise la tolérance % TOL pour l’erreur absolue (1.e-6 par défaut), un % nombre maximal d’itérations NMAX (100 par défaut), % et démarre d’un vecteur initial X0. % [LAMBDA,V,ITER]=EIGPOWER (A,TOL ,NMAX,X0) retourne % aussi le vecteur propre V tel que A*V=LAMBDA*V et le % numéro de l’itération à laquelle V a été calculé . [n,m] = size(A); if n ~= m, error(’Matrices carrées seulement’); end if nargin == 1 tol = 1.e -06; x0 = ones(n,1); nmax = 100; end x0 = x0/norm(x0); pro = A*x0; lambda = x0 ’*pro; err = tol*abs(lambda) + 1; iter = 0; while err >tol*abs(lambda) & abs(lambda )~=0 & iter<=nmax x = pro; x = x/norm(x); pro = A*x; lambdanew = x’*pro; err = abs( lambdanew - lambda); lambda = lambdanew; iter = iter + 1; end return

Algorithme de la Méthode Itérative Générale pour la résolution des systèmes d’équations linéaires de type A*X=B en Matlab

Algorithme de la Méthode Itérative Générale pour la résolution des systèmes d’équations linéaires de type A*X=B en Matlab :

function [x, iter]= itermeth (A,b,x0,nmax,tol ,P) %ITERMETH Méthode itérative générale % X = ITERMETH(A,B,X0,NMAX,TOL ,P) tente de résoudre le % système d’équations linéaires A*X=B d’inconnue X. % La matrice A, de taille NxN , doit etre inversible et % le second membre B doit être de longueur N. % P=’J’ sélectionne la methode de Jacobi, P=’G’ celle % de Gauss-Seidel. Autrement , P est une matrice N x N % qui joue le rôle de préconditionneur dans la methode % de Richardson dynamique. % Les itérations s’arrêtent quand le rapport entre la % norme du k-ème residu et celle du résidu initial est % inférieure ou égale à TOL , le nombre d’itérations % effectuées est alors renvoyé dans ITER. % NMAX est le nombre maximum d’itérations. Si P % n’est pas défini , c’est la méthode du Gradient à % pas optimal qui est utilisée [n,n]=size(A); if nargin == 6 if ischar(P)==1 if P==’J’ L=diag(diag(A)); U=eye(n); beta=1; alpha =1; elseif P == ’G’ L=tril(A); U=eye(n); beta=1; alpha =1; end else [L,U]=lu(P); beta = 0; end else L = eye(n); U = L; beta = 0; end iter = 0; x = x0; r = b - A * x0; r0 = norm(r); err = norm (r); while err > tol & iter < nmax iter = iter + 1; z = L\r; z = U\z; if beta == 0 alpha = z’*r/(z’*A*z); end x = x + alpha*z; r = b - A * x; err = norm (r) / r0; end return

Algorithme de la Factorisation de Gauss en Matlab

Algorithme de la Factorisation de Gauss en Matlab :

function A=lugauss (A) %LUGAUSS Factorisation LU sans pivot. % A = LUGAUSS(A) stocke une matrice triangulaire % supérieure dans la partie triangulaire supérieure de % A et une matrice triangulaire inférieure dans la % partie strictement triangulaire inférieure A (les % termes diagonaux de L valant 1). [n,m]=size(A); if n ~= m error(’A n’’est pas une matrice carrée’); else for k = 1:n-1 for i = k+1:n A(i,k) = A(i,k)/A(k,k); if A(k,k) == 0, error(’Elément diagonal nul’); end j = [k+1:n]; A(i,j) = A(i,j) - A(i,k)*A(k,j); end end end return

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

Algorithme de calcul de Spline d’Interpolation Cubique en Matlab

Algorithme de calcul de Spline d’Interpolation Cubique en Matlab :

function s= cubicspline(x,y,zi,type,der) %CUBICSPLINE calcule une spline cubique % S=CUBICSPLINE(X,Y,ZI) calcule la valeur aux abscisses % ZI de la spline d’interpolation cubique naturelle qui % interpole les valeurs Y aux noeuds X. % S=CUBICSPLINE(X,Y,ZI,TYPE,DER) si TYPE=0 calcule la % valeur aux abscisses ZI de la spline cubique % interpolant les valeurs Y et dont la dérivée % première aux extrémités vaut DER(1) et DER(2). % Si TYPE=1 alors DER(1) et DER(2) sont les valeurs de % la dérivée seconde aux extrémités. [n,m]=size(x); if n == 1 x = x’; y = y’; n = m; end if nargin == 3 der0 = 0; dern = 0; type = 1; else der0 = der(1); dern = der(2); end h = x(2:end)-x(1:end -1); e = 2*[h(1); h(1:end -1)+h(2:end); h(end)]; A = spdiags ([[h; 0] e [0; h]],-1:1,n,n); d = (y(2:end)-y(1:end -1))./h; rhs = 3*(d(2:end)-d(1:end -1)); if type == 0 A(1 ,1) = 2*h(1); A(1 ,2) = h(1); A(n,n) = 2*h(end); A(end ,end -1) = h(end); rhs = [3*(d(1)-der0); rhs; 3*(dern-d(end))]; else A(1 ,:) = 0; A(1 ,1) = 1; A(n,:) = 0; A(n,n) = 1; rhs = [der0; rhs; dern]; end S = zeros(n,4); S(: ,3) = A\rhs; for m = 1:n-1 S(m,4) = (S(m+1,3)-S(m,3))/3/ h(m); S(m,2) = d(m) - h(m)/3*(S(m + 1 ,3)+2*S(m ,3)); S(m,1) = y(m); end S = S(1:n-1, 4: -1:1); pp = mkpp(x,S); s = ppval(pp,zi); return

Algorithme de la Méthode de Newton-Hörner pour le calcul des racines d'un polynôme en Matlab

Algorithme de la Méthode de Newton-Hörner pour le calcul des racines d'un polynôme en Matlab :

function [roots,iter]=newtonhorner(a,x0,tol ,nmax) %NEWTONHORNER méthode de Newton-Horner % [ROOTS,ITER]=NEWTONHORNER(A,X0) calcule les racines % du polynôme % P(X) = A(1)*X^N + A(2)*X^(N -1)+...+ A(N)*X + A(N+1) % en utilisant la méthode de Newton -Horner démarrant % d’une donnée initiale X0. L’algorithme s’arrête % après 100 iterations ou quand la valeur absolue de % la différence entre deux itérées consécutives est % plus petite que 1.e-04. % [ROOTS,ITER]=NEWTONHORNER(A,X0,TOL ,NMAX) permet de % définir la tolérance pour le critère d’arrêt et le % nombre maximal d’itérations. if nargin == 2 tol = 1.e-04; nmax = 100; elseif nargin == 3 nmax = 100; end n=length(a)-1; roots = zeros(n,1); iter = zeros(n,1); for k = 1:n % Itération de Newton niter = 0; x = x0; diff = tol + 1; while niter < nmax & diff >= tol [pz,b] = horner(a,x); [dpz ,b] = horner(b,x); xnew = x - pz/dpz; diff = abs(xnew-x); niter = niter + 1; x = xnew; end if (niter==nmax & diff> tol) fprintf ([’Ne converge pas après avoir atteint ’ ,... ’le nombre maximum d’’itérations\n’]); end % Déflation [pz,a] = horner(a,x); roots(k) = x; iter(k) = niter; end return