tel/TPmatlab/equation_differentielle/TP5a/resout_euler_implicite.m

   1 function Y=resout_euler_implicite(N,T,y0,fcn,nmax,epsilon)\r
   2 % \r
   3 %\r
   4 %\r
   5 %       resoud_euler_implicite : résolution d'une équation différentielle par la méthode d'Euler implicite.\r
   6 %\r
   7 % *********************************************************\r
   8 %\r
   9 %    X=resout_euler_implicite(N,T,y0,fcn,nmax,epsilon) calcule les valeurs aux instants 0,h,2h,...,Nh=T par la méthode\r
  10 %    d'Euler implicite pour l'équation différentielle y't)=f(t,y(t)) sur [0,T] et y(0)=y0.\r
  11 %    A chaque itération, l'équation y_{n+1}=y_n+hf(t_{n+1},y_{n+1}) est résolue par la méthode de Newton\r
  12 %    où la valeur initiale est choisie égale à y_n.\r
  13 %    \r
  14 %    La dérivée partielle df/dy est calculée en symbolique.\r
  15 %\r
  16 %       * variables d'entrées :\r
  17 %      N : nombre de valeurs calculées ; \r
  18 %      T : borne supérieure de l'intervalle ;\r
  19 %      y0 : valeur de y à t=0 ;\r
  20 %      fcn : chaine de caractères (de type 'f(t,y)') représentant la fonction f\r
  21 %      epsilon : précision souhaitée (pour la méthode de Newton)\r
  22 %      nmax : nombre maximal d'itérations (pour la méthode de Newton)\r
  23 %\r
  24 %       * variables de sortie :\r
  25 %     Y : valeurs calculée aux instants 0,h,2h,...,Nh=T.\r
  26 %\r
  27 % \r
  28 % \r
  29 % ************ Fonctions auxiliaires utilisées ************\r
  30 %\r
  31 %       aucune\r
  32 %\r
  33 % *********************************************************\r
  34 %\r
  35 \r
  36 % Contrôles d'entrée\r
  37 \r
  38 % nombre d'arguments\r
  39 if nargin~=6\r
  40    error('nombre d''arguments de la fonction incorrect');\r
  41 end\r
  42    \r
  43 % Corps d'algorithme\r
  44 h=T/N;\r
  45 Y=zeros(1,N+1);\r
  46 Y(1)=y0;\r
  47 auxiY=y0;\r
  48 syms t y;\r
  49 g=eval(fcn);\r
  50 fcnpsymb=diff(g,y);\r
  51 fcnp=char(fcnpsymb);\r
  52 yapp=y0;\r
  53 for i=1:N\r
  54    t=h*i;\r
  55    y=yapp;\r
  56    fx=-h*eval(fcn);\r
  57    fxp=1-h*eval(fcnp);\r
  58    test=((fx)~=0);\r
  59    n=0;\r
  60    v=yapp;\r
  61    while(test)\r
  62       oldv=v;\r
  63       v=v-fx/fxp;\r
  64       y=v;\r
  65       fx=v-yapp-h*eval(fcn);\r
  66       fxp=1-h*eval(fcnp);\r
  67       n=n+1;\r
  68       test=~((fx==0)|(n>=nmax)|(abs(v-oldv)<=epsilon));\r
  69    end   \r
  70    yapp=v;\r
  71    if (n>=nmax)\r
  72       error('attention, n>nmax pour la méthode de Newton : arrêt du programme');\r
  73    end\r
  74    Y(i+1)=yapp;\r
  75 end\r
  76 \r
  77 \r
  78 \r
  79 \r