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 |
Aucun commentaire:
Enregistrer un commentaire