tel/TPmatlab/equation_differentielle/TP5b/resout_euler_implicite_bis.m

   1 function Y=resout_euler_implicite_bis(N,T,y0,fcn,fcnp,nmax,epsilon)\r
   2 % \r
   3 %\r
   4 %\r
   5 %       resoud_euler_implicite_bis : 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_bis(N,T,y0,fcn,fcnp,nmax,epsilon) calcule les valeurs aux instants 0,h,2h,...,Nh=T \r
  10 %    par la méthode 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 entrée comme paramètre.\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 \r
  21 %      fcnp : chaine de caractères (de type 'f(t,y)') représentant la fonction df/dy\r
  22 %      epsilon : précision souhaitée (pour la méthode de Newton)\r
  23 %      nmax : nombre maximal d'itérations (pour la méthode de Newton)\r
  24 %\r
  25 %       * variables de sortie :\r
  26 %     Y : valeurs calculée aux instants 0,h,2h,...,Nh=T.\r
  27 %\r
  28 % \r
  29 % \r
  30 % ************ Fonctions auxiliaires utilisées ************\r
  31 %\r
  32 %       aucune\r
  33 %\r
  34 % *********************************************************\r
  35 %\r
  36 \r
  37 % Contrôles d'entrée\r
  38 \r
  39 % nombre d'arguments\r
  40 if nargin~=7\r
  41    error('nombre d''arguments de la fonction incorrect');\r
  42 end\r
  43    \r
  44 % Corps d'algorithme\r
  45 h=T/N;\r
  46 Y=zeros(1,N+1);\r
  47 Y(1)=y0;\r
  48 auxiY=y0;\r
  49 yapp=y0;\r
  50 for i=1:N\r
  51    t=h*i;\r
  52    y=yapp;\r
  53    fx=-h*eval(fcn);\r
  54    fxp=1-h*eval(fcnp);\r
  55    test=((fx)~=0);\r
  56    n=0;\r
  57    v=yapp;\r
  58    while(test)\r
  59       oldv=v;\r
  60       v=v-fx/fxp;\r
  61       y=v;\r
  62       fx=v-yapp-h*eval(fcn);\r
  63       fxp=1-h*eval(fcnp);\r
  64       n=n+1;\r
  65       test=~((fx==0)|(n>=nmax)|(abs(v-oldv)<=epsilon));\r
  66    end   \r
  67    yapp=v;\r
  68    if (n>=nmax)\r
  69       error('attention, n>nmax pour la méthode de Newton : arrêt du programme');\r
  70    end\r
  71    Y(i+1)=yapp;\r
  72 end\r
  73 \r
  74 \r
  75 \r
  76 \r