tel/TPmatlab/equation_differentielle/TP5d/demo_trace_systeme_mecanique.m

   1 % ce script demo_trace_systeme_mecanique calcule les \r
   2 % solutions du système à deux degrés de liberté régi par le système : \r
   3 % x1''(t)+(k1+k2)x1-k2x2+c1x'1+d1(x'1)^3=f1\r
   4 % x2''(t)-k2x1+(k2+k3)x2+c2x'2+d2(x'2)^3=f2\r
   5 % déterminée par les schémas d'Euler implicite et de Runge Kutta 4.\r
   6 %\r
   7 %\r
   8 %  * les paramètres d'entrée sont : \r
   9 %    k : le vecteur contenant k1, k2 et k3,\r
  10 %    c : le vecteur contenant c1 et c2,\r
  11 %    d : le vecteur contenant d1, d2,\r
  12 %    y0 : le vecteur contenant les conditions initiales (x1(0), x1'(0), x2(0), x2'(0))\r
  13 %    fcn1 : chaine de caractères (de type 'f1(t)') représentant la fonction f1,  \r
  14 %    fcn2 : chaine de caractères (de type 'f2(t)') représentant la fonction f2,  \r
  15 %    N : un entier non nul,\r
  16 %    T : la borne supérieur de l'intervalle d'étude.\r
  17 %\r
  18 %  * le script affiche les approximation de x1, x1', x2 et x2' sur [0,T].\r
  19 \r
  20 clear all;\r
  21 %nombre de points pour l'affichage défini par ngra \r
  22 ngra=1000;\r
  23 \r
  24 % entree des paramètres\r
  25 k=input('entrez le vecteur [k1,k2,k3] : ');\r
  26 c=input('entrez  le vecteur [c1,c2]: ');\r
  27 d=input('entrez  le vecteur [d1,d2]: ');\r
  28 y0=input('entrez  le vecteur des conditions initiales [x1(0),x1''(0),x2(0),x2''(0)]: ');\r
  29 fcn1=input('entrez  la fonction f1 sous la forme ''f1(t)'': ');\r
  30 fcn2=input('entrez  la fonction f2 sous la forme ''f2(t)'': ');\r
  31 N=input('entrez  le nombre de points de discretisation : ');\r
  32 T=input('entrez  T : ');\r
  33 \r
  34 \r
  35 % définition de la fonction F telle que le système soit équivalent à Y'(t)=F(t,Y(t)) dans R^4.\r
  36 % F est defini en fonction de t et de y, vecteur à 4 composantes.\r
  37 A=-[0 -1 0 0;\r
  38    k(1)+k(2) c(1) -k(2) 0;\r
  39    0  0 0 -1;\r
  40    -k(2) 0 k(2)+k(3) c(2)];\r
  41 H='-[0;d(1)*(y(2))^3;0;d(2)*(y(4))^3]';\r
  42 G=['[0; ',fcn1,';','0;',fcn2,']'];\r
  43 fcn=[H,'+',G,'+A*y'];\r
  44 h=T/N;\r
  45 \r
  46 % calcul de la solution par Euler explicite.\r
  47 Yeuler=zeros(4,N+1);\r
  48 Yeuler(:,1)=y0';\r
  49 auxiY=y0';\r
  50 for i=1:N\r
  51    t=h*(i-1);\r
  52    y=auxiY;\r
  53    auxiY=auxiY+h*eval(fcn);\r
  54    Yeuler(:,i+1)=auxiY;\r
  55 end\r
  56 \r
  57 % calcul de la solution par Runge Kutta.\r
  58 Yrungekutta=zeros(4,N+1);\r
  59 Yrungekutta(:,1)=y0';\r
  60 auxiY=y0';\r
  61 for i=1:N\r
  62    t=h*(i-1);\r
  63    y=auxiY;\r
  64    k1=h*eval(fcn);\r
  65    t=t+h/2;\r
  66    y=auxiY+k1/2;\r
  67    k2=h*eval(fcn);\r
  68    y=auxiY+k2/2;\r
  69    k3=h*eval(fcn);\r
  70    t=h*i;\r
  71    y=auxiY+k3;\r
  72    k4=h*eval(fcn);\r
  73    auxiY=auxiY+(k1+2*k2+2*k3+k4)/6;\r
  74    Yrungekutta(:,i+1)=auxiY;\r
  75 end\r
  76 \r
  77 \r
  78 % préparation des vecteurs de tracé.\r
  79 if (N<ngra)\r
  80     x=0:h:T;\r
  81     yeuler=Yeuler;\r
  82     yrungekutta=Yrungekutta;\r
  83 else\r
  84     p=fix(N/ngra);\r
  85     yeuler=[Yeuler(:,1:p:N+1) Yeuler(:,N+1)];\r
  86     yrungekutta=[Yrungekutta(:,1:p:N+1) Yrungekutta(:,N+1)];\r
  87     x=0:h:T;\r
  88     x=[x(1:p:N+1) T];\r
  89 end\r
  90 \r
  91 % tracé\r
  92 clf;\r
  93 subplot(2,1,1); \r
  94 plot(x,yeuler(1,:),x,yeuler(3,:),x,yrungekutta(1,:),x,yrungekutta(3,:),[0 T],[0 0],'k');\r
  95 legend('x1 par Euler','x2 par Euler','x1 par RK4','x2 par RK4',0);\r
  96 subplot(2,1,2); \r
  97 plot(x,yeuler(2,:),x,yeuler(4,:),x,yrungekutta(2,:),x,yrungekutta(4,:),[0 T],[0 0],'k');\r
  98 legend('x1'' par Euler','x2'' par Euler','x1'' par RK4','x2'' par RK4',0);\r
  99 \r
 100 \r
 101 \r
 102 \r
 103 \r
 104 \r
 105 \r
 106 \r