tel/TPmatlab/equation_differentielle/TP5b/resout_theta_methode_bis.m

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