tel/TPmatlab/equation_differentielle/TP5a/resout_runge_kutta4.m

   1 function Y=resout_runge_kutta4(N,T,y0,fcn)\r
   2 % \r
   3 %\r
   4 %\r
   5 %       resout_runge_kutta4 : résolution d'une équation différentielle par la méthode de Runge Kutta 4 (explicite).\r
   6 %\r
   7 % *********************************************************\r
   8 %\r
   9 %    Y=resout_runge_kutta4(N,T,y0,fcn) calcule les valeurs aux instants 0,h,2h,...,Nh=T par la méthode\r
  10 %    de Runge Kutta d'ordre 4 (explicite) pour l'équation différentielle y't)=f(t,y(t)) sur [0,T] et y(0)=y0.\r
  11 %       \r
  12 %       * variables d'entrées :\r
  13 %      N : nombre de valeurs calculées ; \r
  14 %      T : borne supérieure de l'intervalle ;\r
  15 %      y0 : valeur de y à t=0 ;\r
  16 %      fcn : chaine de caractères (de type 'f(t,y)') représentant la fonction f.\r
  17 %\r
  18 %       * variables de sortie :\r
  19 %     Y : valeurs calculée aux instants 0,h,2h,...,Nh=T.\r
  20 %\r
  21 % \r
  22 % \r
  23 % ************ Fonctions auxiliaires utilisées ************\r
  24 %\r
  25 %       aucune\r
  26 %\r
  27 % *********************************************************\r
  28 %\r
  29 \r
  30 \r
  31 \r
  32 % Contrôles d'entrée\r
  33 \r
  34 % nombre d'arguments\r
  35 if nargin~=4\r
  36    error('nombre d''arguments de la fonction incorrect');\r
  37 end\r
  38    \r
  39 % Corps d'algorithme\r
  40 h=T/N;\r
  41 Y=zeros(1,N+1);\r
  42 Y(1)=y0;\r
  43 auxiY=y0;\r
  44 for i=1:N\r
  45    t=h*(i-1);\r
  46    y=auxiY;\r
  47    k1=h*eval(fcn);\r
  48    t=t+h/2;\r
  49    y=auxiY+k1/2;\r
  50    k2=h*eval(fcn);\r
  51    y=auxiY+k2/2;\r
  52    k3=h*eval(fcn);\r
  53    t=h*i;\r
  54    y=auxiY+k3;\r
  55    k4=h*eval(fcn);\r
  56    auxiY=auxiY+(k1+2*k2+2*k3+k4)/6;\r
  57    Y(i+1)=auxiY;\r
  58 end\r
  59 \r
  60 \r
  61 \r
  62 \r