tel/TPmatlab/equation_differentielle/TP5a/resout_theta_methode.m

   1 function Y=resout_theta_methode(theta,N,T,y0,fcn,nmax,epsilon)\r
   2 % \r
   3 %\r
   4 %\r
   5 %       resout_theta_methode : résolution d'une équation différentielle par la theta-méthode.\r
   6 %\r
   7 % *********************************************************\r
   8 %\r
   9 %    Y=resout_theta_methode(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 calculée en symbolique.\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 %      epsilon : précision souhaitée (pour la méthode de Newton)\r
  25 %      nmax : nombre maximal d'itérations (pour la méthode de Newton)\r
  26 %\r
  27 %       * variables de sortie :\r
  28 %     Y : valeurs calculée aux instants 0,h,2h,...,Nh=T.\r
  29 %\r
  30 % \r
  31 % \r
  32 % ************ Fonctions auxiliaires utilisées ************\r
  33 %\r
  34 %       aucune\r
  35 %\r
  36 % *********************************************************\r
  37 %\r
  38 \r
  39 % Contrôles d'entrée\r
  40 \r
  41 % nombre d'arguments\r
  42 if nargin~=7\r
  43    error('nombre d''arguments de la fonction incorrect');\r
  44 end\r
  45    \r
  46 % Corps d'algorithme\r
  47 h=T/N;\r
  48 Y=zeros(1,N+1);\r
  49 Y(1)=y0;\r
  50 auxiY=y0;\r
  51 syms t y;\r
  52 g=eval(fcn);\r
  53 fcnpsymb=diff(g,y);\r
  54 fcnp=char(fcnpsymb);\r
  55 yapp=y0;\r
  56 for i=1:N\r
  57    y=yapp;\r
  58    t=h*(i-1);\r
  59    auxidud=h*(1-theta)*eval(fcn);\r
  60    t=h*i;\r
  61    fx=-h*theta*eval(fcn)-auxidud;\r
  62    fxp=1-h*theta*eval(fcnp);\r
  63    test=((fx)~=0);\r
  64    n=0;\r
  65    v=yapp;\r
  66    while(test)\r
  67       oldv=v;\r
  68       v=v-fx/fxp;\r
  69       y=v;\r
  70       t=h*i;\r
  71       fx=v-yapp-h*theta*eval(fcn)-auxidud;\r
  72       fxp=1-h*theta*eval(fcnp);\r
  73       n=n+1;\r
  74       test=~((fx==0)|(n>=nmax)|(abs(v-oldv)<=epsilon));\r
  75    end   \r
  76    yapp=v;\r
  77    if (n>=nmax)\r
  78       error('attention, n>nmax pour la méthode de Newton : arrêt du programme');\r
  79    end\r
  80    Y(i+1)=yapp;\r
  81 end\r