11-2 ODE ���O�򥻥Ϊk

¦b¨Ï¥Î ODE «ü¥O®É¡A§Ú­Ì¥²¶·¥ý±N­n¨D¸Ñªº ODE ªí¥Ü¦¨¤@­Ó¨ç¦¡¡A¦¹¨ç¦¡ªº¿é¤J¬° t¡]®É¶¡¡^¤Î y¡]ª¬ºAÅܼơAState Variables¡^¡A¿é¥X«h¬° dy¡]ª¬ºAÅܼƪº·L¤À­È¡^¡C¦pªG¦¹ ODE ¨ç¦¡¡]©Î¤SºÙ ODE Àɮס^ªºÀɦW¬° odeFile.m¡A«h§Ú­Ì©I¥s ODE «ü¥Oªº®æ¦¡¦p¤U¡G

[t, y] = solver ('odeFile', [t0, t1], y0);

¨ä¤¤ [t0, t1] ¬O¿n¤Àªº®É¶¡°Ï¶¡¡Ay0 ¥Nªí°_©l±ø¥ó¡]Initial Conditions¡^¡Asolver ¬O«eªí©Ò¦Cªº¦UºØ ODE «ü¥O¡At ¬O¿é¥Xªº®É¶¡¦V¶q¡Ay «h¬O¹ïÀ³ªºª¬ºAÅܼƦV¶q¡C

¥H van der Pol ·L¤À¤èµ{¬°¨Ò¡A¨ä¤èµ{¦¡¬°¡G

$$y''-\mu (1-y^2)y'+y=0$$

­º¥ý§Ú­Ì¥²¶·±N¥¦¤Æ¦¨¼Ð·Ç®æ¦¡¡A¥O $y_1=y$¡A$y_2=y'$¡A«h¤W­z·L¤À¤èµ{¦¡¥i¥H§ï¼g¬°¨â­Ó¤@¶¥ªº·L¤À¤èµ{¦¡¡G

$$ \left\{ \begin{array}{rcl} y_1' & = & y_2 \\ y_2' & = & \mu (1-y_1^2)y_2-y_1\\ \end{array} \right. $$

¦¹ºØ¼Ð·Ç®æ¦¡¥i¥Î¦V¶q¨Óªí¥Ü¦¨¤@¯ë¤Æªº¼Æ¾Ç¦¡¡G

$$\vec{y}' = \vec{F}(\vec{y}, t)$$

¨ä¤¤ $\vec{y}'$ ¬°¤@¦V¶q¡A¥Nªíª¬ºAÅܼơC¨Ï¥Î¦¹¼Ð·Ç®æ¦¡«á¡A¨Ã°²³] $\mu=1$¡A§Ú­Ìªº ODE Àɮס]vdp1.m¡^¥iÅã¥Ü¦p¤U¡G

11-±`·L¤À¤èµ{¦¡/vdp1.mfunction dy = vdp1(t, y) mu = 1; dy = [y(2); mu*(1-y(1)^2)*y(2)-y(1)];

¤@¥¹¦³¤F ODE ÀɮסA§Ú­Ì§Y¥i¿ï¥Î«e­z¤§ ODE «ü¥O¨Ó¥º¸Ñ¡C¦b $\mu=1$ ®É¡Avan der Pol ¤èµ{¦¡¨Ã«D Stiff ¨t²Î¡A¦]¦¹¡A§Ú­Ì¥i¿ï¥Î ode45 ¨Ó¨D¸Ñ¡A³Ì²³æªº§@ªk¡A´N¬O¨Ï¥Î¦¹«ü¥O¨Óµe¥X¿n¤Àªºµ²ªG¡A½d¨Ò¦p¤U¡G

Example 1: 11-±`·L¤À¤èµ{¦¡/odeBasic01.mode45('vdp1', [0 25], [3 3]');

¨ä¤¤ [0, 25] ¥Nªí¿n¤Àªº®É¶¡°Ï¶¡¡A[3 3]¡¦ «h¥Nªí°_©l±ø¥ó¡]¥²¶·¥H¦æ¦V¶q¨Óªí¥Ü¡^¡C¥Ñ©ó¨S¦³¿é¥XÅܼơA©Ò¥H¤W­zµ{¦¡°õ¦æµ²§ô«á¡AMATLAB ¥u·|µe¥Xª¬ºAÅܼƹï®É¶¡ªº¹Ï§Î¡A¦p¤W¹Ï¡C­Y­n¨ú±o¿n¤À©Ò±oªºª¬ºAÅܼƤιïÀ³ªº®É¶¡¡A¥i¥H¥[¤W¿é¥XÅܼƥH¨ú±o³o¨Ç¸ê®Æ¡A½d¨Ò¦p¤U¡G

Example 2: 11-±`·L¤À¤èµ{¦¡/odeGetData01.m[t, y] = ode45('vdp1', [0 25], [3 3]'); plot(t, y(:,1), t, y(:,2)); xlabel('Time t'); ylabel('Solution y(t) and y''(t)'); legend('y(t)', 'y''(t)');

§Ú­Ì¤]¥i¥Hµe¥X $y(t)$ ¤Î $y'(t)$ ¦b¬Û¦ì¥­­±¡]Phase Plane ¡^ªº¹B°Ê±¡ªp¡A¨Ò¦p¡G

Example 3: 11-±`·L¤À¤èµ{¦¡/odephasePlot01.m[t, y] = ode45('vdp1', [0 25], [3 3]'); plot(y(:,1), y(:,2), '-o'); xlabel('y(t)'); ylabel('y''(t)');

¦b¤W¨Ò¤¤¡A·í $\mu$ ­È¶V¨Ó¶V¤j®É¡Avan der Pol ¤èµ{¦¡´Nº¥º¥Åܦ¨¤@­Ó Stiff ªº¨t²Î¡A¦¹®É­Y­n¸Ñ¦¹¤èµ{¦¡¡A´N¥²¶·§ï¥Î±Mªù¹ï¥I Stiff ¨t²Îªº«ü¥O¡C

­º¥ý¡A§Ú­Ì±N $\mu$ ­È§ï¦¨ 1000¡A ODE ÀÉ®×­n§ï¼g¦p¤U¡]vdp2.m¡^¡G

11-±`·L¤À¤èµ{¦¡/vdp2.mfunction dy = vdp2(t, y) mu = 1000; dy = [y(2); mu*(1-y(1)^2)*y(2)-y(1)];

¦¹®É§Y¥i¿ï¥Î±Mªù¹ï¥I Stiff ¨t²Îªº ODE «ü¥O¡A¨Ò¦p ode15s¡A¨Ó¨D¸Ñ¦¹¨t²Î¨Ã§@¹ÏÅã¥Ü¡A¦p¤U¡G

Example 4: 11-±`·L¤À¤èµ{¦¡/ode15s01.m[t, y]= ode15s('vdp2', [0 3000], [2 1]'); subplot(2,1,1); plot(t, y(:,1), '-o'); xlabel('Time t'); ylabel('y(t)'); subplot(2,1,2); plot(t, y(:,2), '-o'); xlabel('Time t'); ylabel('y''(t)'); % ª`·N³æ¤Þ¸¹¡u'¡vªº­«ÂШϥÎ

¥Ñ¤W¹Ï¥i¬Ý¥X¡A$y'(t)$ ªºÅܤƬ۷í¼@¯P¡]¶W¹L $\pm 1000$¡^¡A³o´N¬O Stiff ¨t²Îªº¯S¦â¡C

­Y­nµe¥X¤Gºû¥­­±¬Û¦ì¹Ï¡A¥i¥H¨Ï¥Î¤U¦C½d¨Ò¡G

Example 5: 11-±`·L¤À¤èµ{¦¡/ode15s02.m[t, y]= ode15s('vdp2', [0 3000], [2 1]'); subplot(1,1,1); plot(y(:, 1), y(:, 2), '-o'); xlabel('y(t)'); ylabel('y''(t)')

­Y­n²£¥Í¦b¬Y¨Ç¯S©w®É¶¡ÂIªºª¬ºAÅܼƭȡA«h§Ú­Ì©I¥s ODE «ü¥Oªº®æ¦¡¥i§ï¦¨¡G

[t, y] = solver('odeFile', [t0, t1, ..., tn], y0);

¨ä¤¤ [t0, t1, ..., tn] §Y¬O¯S©w®É¶¡ÂI©Ò§Î¦¨ªº¦V¶q¡A½ÐŪªÌ¦Û¦æ¸Õ¸Õ¬Ý¡C


MATLABµ{¦¡³]­p¡G¶i¶¥½g