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