tel/TPmatlab/equation_differentielle/TP5b/demo_log_erreur_schema.m

   1 % demo_log_erreur_schema : ce script permet de déterminer \r
   2 % l'ordre du schéma de façon graphique\r
   3 % et numérique des cinq méthodes suivantes : \r
   4 % * Euler explicite ;\r
   5 % * Euler implicite ;\r
   6 % * Runge Kutta explicite d'ordre 2\r
   7 % * theta_méthode\r
   8 % * Runge Kutta explicite d'ordre 4\r
   9 %\r
  10 % paramètres d'entrée : \r
  11 %      nmin et nmax : deux entiers entre lesquels varie n ;\r
  12 %      T : borne supérieure de l'intervalle ;\r
  13 %      y0 : valeur de y à t=0 ;\r
  14 %      theta est un réel dans [0,1] ;\r
  15 %      epsilon : précision souhaitée (pour la méthode de Newton, pour les méthodes implicites)\r
  16 %      nmax : nombre maximal d'itérations (pour la méthode de Newton, pour les méthodes implicites)\r
  17 %      fcn : chaine de caractères (de type 'f(t,y)') représentant la fonction f\r
  18 \r
  19 pmin=input('entrez l''exposant pmin tel que le plus petit entier de calcul soit 10^pmin : ');\r
  20 pmax=input('entrez l''exposant pmax tel que le plus grand entier de calcul soit 10^pmax : ');\r
  21 nc=input('entrez le nombre de calcul : ');\r
  22 T=input('entrez la borne supérieur de l''intervalle : ');\r
  23 y0=input('entrez la condition initiale : ');\r
  24 theta=input('entrez le paramètre theta : ');\r
  25 epsilon=input('entrez epsilon (pour la méthode de Newton, pour les méthodes implicites) : ');\r
  26 nmax=input('entrez le nombre maximum d''itérations (pour la méthode de Newton, pour les méthodes implicites) : ');\r
  27 fcn=input('entrez fcn (sous la forme ''f(t,y)'') : ');\r
  28 \r
  29 syms t y;\r
  30 g=eval(fcn);\r
  31 fcnpsymb=diff(g,y);\r
  32 fcnp=char(fcnpsymb);\r
  33 x=[];\r
  34 y1=[];\r
  35 y2=[];\r
  36 y3=[];\r
  37 y4=[];\r
  38 y5=[];\r
  39 for N=logspace(pmin,pmax,nc)\r
  40    NN=fix(N);\r
  41    x=[x log10(NN)];\r
  42    Ya1=resout_euler_explicite(NN,T,y0,fcn);\r
  43    Yb1=resout_euler_explicite(2*NN,T,y0,fcn);\r
  44    Ya2=resout_euler_implicite_bis(NN,T,y0,fcn,fcnp,nmax,epsilon);\r
  45    Yb2=resout_euler_implicite_bis(2*NN,T,y0,fcn,fcnp,nmax,epsilon);\r
  46    Ya3=resout_runge_kutta2(NN,T,y0,fcn); \r
  47    Yb3=resout_runge_kutta2(2*NN,T,y0,fcn); \r
  48    Ya4=resout_theta_methode_bis(theta,NN,T,y0,fcn,fcnp,nmax,epsilon);\r
  49    Yb4=resout_theta_methode_bis(theta,2*NN,T,y0,fcn,fcnp,nmax,epsilon);\r
  50    Ya5=resout_runge_kutta4(NN,T,y0,fcn); \r
  51    Yb5=resout_runge_kutta4(2*NN,T,y0,fcn); \r
  52    auxi1=-log10(max(abs(Ya1-Yb1(1:2:end)))); \r
  53    auxi2=-log10(max(abs(Ya2-Yb2(1:2:end)))); \r
  54    auxi3=-log10(max(abs(Ya3-Yb3(1:2:end)))); \r
  55    auxi4=-log10(max(abs(Ya4-Yb4(1:2:end)))); \r
  56    auxi5=-log10(max(abs(Ya5-Yb5(1:2:end)))); \r
  57    y1=[y1 auxi1]; \r
  58    y2=[y2 auxi2]; \r
  59    y3=[y3 auxi3]; \r
  60    y4=[y4 auxi4]; \r
  61    y5=[y5 auxi5]; \r
  62    disp(' ');\r
  63    disp(['calcul pour N=',int2str(NN),' terminé']);\r
  64 end\r
  65 [a1,b,r1]=regression_lineaire(x,y1);\r
  66 [a2,b,r2]=regression_lineaire(x,y2);\r
  67 [a3,b,r3]=regression_lineaire(x,y3);\r
  68 [a4,b,r4]=regression_lineaire(x,y4);\r
  69 [a5,b,r5]=regression_lineaire(x,y5);\r
  70 clf;\r
  71 hold on;\r
  72 plot(x,y1,x,y2,x,y3,x,y4,x,y5);\r
  73 legend('Euler explicite','Euler implicite','runge Kutta 2','theta-méthode','runge Kutta 4',0);\r
  74 title('détermination des ordres des schémas');\r
  75 hold off;\r
  76 disp(' ');\r
  77 disp('pour avoir une estimation des ordres, tapez sur une touche');\r
  78 pause;\r
  79 disp(' ');\r
  80 disp('Euler explicite');\r
  81 disp(['ordre=',num2str(a1)]);\r
  82 disp(['correlation=',num2str(r1)]);\r
  83 disp(' ');\r
  84 disp('Euler implicite');\r
  85 disp(['ordre=',num2str(a2)]);\r
  86 disp(['correlation=',num2str(r2)]);\r
  87 disp(' ');\r
  88 disp('runge Kutta 2');\r
  89 disp(['ordre=',num2str(a3)]);\r
  90 disp(['correlation=',num2str(r3)]);\r
  91 disp(' ');\r
  92 disp('theta-méthode');\r
  93 disp(['ordre=',num2str(a4)]);\r
  94 disp(['correlation=',num2str(r4)]);\r
  95 disp(' ');\r
  96 disp('runge Kutta 4');\r
  97 disp(['ordre=',num2str(a5)]);\r
  98 disp(['correlation=',num2str(r5)]);\r
  99 \r
 100 \r