📄 ch10_2.htm
字号:
<html><head><title> 阮奇-库达方法 </title><meta NAME="GENERATOR" CONTENT="Microsoft FrontPage 3.0"></head><body BACKGROUND="../img1/bg0000.gif" tppabs="http://webclass.ncu.edu.tw/~junwu/img/bg0000.gif"><script language="JAVASCRIPT"><!--if (navigator.onLine){document.write("<!-- Spidersoft WebZIP Ad Banner Insert -->");document.write("<TABLE width=100% border=0 cellpadding=0 cellspacing=0>");document.write("<TR>");document.write("<TD>");document.write("<ILAYER id=ad1 visibility=hidden height=60></ILAYER>");document.write("<NOLAYER>");document.write("<IFRAME SRC='http://www.spidersoft.com/ads/bwz468_60.htm' width=100% height=60 marginwidth=0 marginheight=0 hspace=0 vspace=0 frameborder=0 scrolling=no></IFRAME>");document.write("</NOLAYER>");document.write("</TD>");document.write("</TR>");document.write("</TABLE>");document.write("<!-- End of Spidersoft WebZIP Ad Banner Insert-->");} //--></script><!-- Spidersoft WebZIP Ad Banner Insert --><!-- End of Spidersoft WebZIP Ad Banner Insert--><h1><font SIZE="6" COLOR="#0000FF">10.2 阮奇-库达方法 </font></h1><hr><p>阮奇-库达 (Runge-Kutta) 方法是最通用的解 ODE 的方法，它可以依计算精确度的要求有低阶到高阶的各个 计算式，举例来说，一阶阮奇-库达法为 </p><p><img SRC="../img10/img00007.gif" tppabs="http://webclass.ncu.edu.tw/~junwu/img10/img00007.gif" WIDTH="89" HEIGHT="22"> </p><p>其实上式即是一阶的泰勒序数近似式，只不过令 <img SRC="../img10/img00008.gif" tppabs="http://webclass.ncu.edu.tw/~junwu/img10/img00008.gif" WIDTH="218" HEIGHT="22"> 。因为一阶阮奇-库达法 的精确度太低，所以在真正解ODE 时，最少须用二阶以上的方法。 <br></p><p>MATLAB应用阮奇-库达法的函数有<font COLOR="#FF0000">ode23</font>, <font COLOR="#FF0000">ode45</font>，其中<font COLOR="#FF0000">ode23</font>是同时以二阶及三阶阮奇-库达法求解，而<font COLOR="#FF0000">ode45</font> 则是以四阶及五阶阮奇-库达法求解。其语法为<font COLOR="#FF0000">ode23('dy',x0,xn,y0)</font>，其中 <font COLOR="#FF0000">dy</font> 为ODE中的等式右边的函数（如之 前介绍的<img SRC="../img10/img00009.gif" tppabs="http://webclass.ncu.edu.tw/~junwu/img10/img00009.gif" WIDTH="62" HEIGHT="22">），<font COLOR="#FF0000">x0</font>, <font COLOR="#FF0000">xn</font> 是要解ODE的区间 [<i>x0, xn</i>] 的二个端点，<font COLOR="#FF0000">y0</font>是起始值 (<i>y0=y</i>(<i>x0</i>))。而<font COLOR="#FF0000">ode45</font>的语法 与<font COLOR="#FF0000">ode23</font>相同。 <br></p><p>先前的四个 ODE 的解法如下： <br></p><p>例一、要在区间 [2, 4] 解以下的 ODE： </p><p><img SRC="../img10/img00010.gif" tppabs="http://webclass.ncu.edu.tw/~junwu/img10/img00010.gif" WIDTH="256" HEIGHT="25"> </p><p><font COLOR="#FF0000">% m-function, g1.m</font> </p><p><font COLOR="#FF0000">function dy=g1(x,y)</font> </p><p><font COLOR="#FF0000">dy=3*x.^2;<br></font></p><p><font COLOR="#FF0000">>> [x,num_y]=ode23('g1',2,4,0.5);</font> </p><p><font COLOR="#FF0000">>> anl_y=x.^3-7.5;</font> </p><p><font COLOR="#FF0000">>> plot(x,num_y,x,anl_y,'o')</font> </p><p><font COLOR="#FF0000">>> title('Solution of g1')</font> </p><p><font COLOR="#FF0000">>> xlabel('x'), ylabel('y=f(x)'), grid <br></font></p><p>例二、要在区间 [0, 5] 解以下的 ODE： </p><p><img SRC="../img10/img00011.gif" tppabs="http://webclass.ncu.edu.tw/~junwu/img10/img00011.gif" WIDTH="276" HEIGHT="23"> </p><p><font COLOR="#FF0000">% m-function, g2.m</font> </p><p><font COLOR="#FF0000">function dy=g2(x,y)</font> </p><p><font COLOR="#FF0000">dy=-0.131*y;<br></font></p><p><font COLOR="#FF0000">>> [x,num_y]=ode23('g2',0,5,4);</font> </p><p><font COLOR="#FF0000">>> anl_y=4*exp(-0.131*x);</font> </p><p><font COLOR="#FF0000">>> plot(x,num_y,x,anl_y,'o')</font> </p><p><font COLOR="#FF0000">>> title('Solution of g2')</font> </p><p><font COLOR="#FF0000">>> xlabel('x'), ylabel('y=f(x)'), grid <br></font></p><p>例三、要在区间 [0, 2] 解以下的 ODE： </p><p><img SRC="../img10/img00012.gif" tppabs="http://webclass.ncu.edu.tw/~junwu/img10/img00012.gif" WIDTH="312" HEIGHT="26"> </p><p><font COLOR="#FF0000">% m-function, g3.m</font> </p><p><font COLOR="#FF0000">function dy=g3(x,y)</font> </p><p><font COLOR="#FF0000">dy=2*x*cos(y)^2;<br></font></p><p><font COLOR="#FF0000">>> [x,num_y]=ode23('g3',0,2,pi/4);</font> </p><p><font COLOR="#FF0000">>> anl_y=atan(x.*x+1);</font> </p><p><font COLOR="#FF0000">>> plot(x,num_y,x,anl_y,'o')</font> </p><p><font COLOR="#FF0000">>> title('Solution of g3')</font> </p><p><font COLOR="#FF0000">>> xlabel('x'), ylabel('y=f(x)'), grid <br></font></p><p>例四、要在区间 [0, 3] 解以下的 ODE： </p><p><a NAME="work"><img SRC="../img10/img00013.gif" tppabs="http://webclass.ncu.edu.tw/~junwu/img10/img00013.gif" WIDTH="278" HEIGHT="26"></a> </p><p><font COLOR="#FF0000">% m-function, g4.m</font> </p><p><font COLOR="#FF0000">function dy=g4(x,y)</font> </p><p><font COLOR="#FF0000">dy=3*y+exp(2*x);<br></font></p><p><font COLOR="#FF0000">>> [x,num_y]=ode23('g4',0,3,3);</font> </p><p><font COLOR="#FF0000">>> anl_y=4*exp(3*x)-exp(2*x);</font> </p><p><font COLOR="#FF0000">>> plot(x,num_y,x,anl_y,'o')</font> </p><p><font COLOR="#FF0000">>> title('Solution of g4')</font> </p><p><font COLOR="#FF0000">>> xlabel('x'), ylabel('y=f(x)'), grid <br></font></p><p>如果将上述方法改以 <font COLOR="#FF0000">ode45</font> 计算，你可能无法察觉出其与<font COLOR="#FF0000">ode23</font>的解之间的差异，那是因为我们选的 ODE 代表的函数分布变化平缓，所以高阶方法就显现不出其优点。不过以二者在计算的误差上做比较，<font COLOR="#FF0000">ode45</font> 的误差量级会比 <font COLOR="#FF0000">ode23</font>要小。以下是几个例子： </p><p><font COLOR="#FF0000">% m-function, g1.m</font> </p><p><font COLOR="#FF0000">function dy=g1(x,y)</font> </p><p><font COLOR="#FF0000">dy=3*x.^2;<br></font></p><p><font COLOR="#FF0000">% m-file, odes1.m</font> </p><p><font COLOR="#FF0000">% Solve an ode using ode23 and ode45</font> </p><p><font COLOR="#FF0000">clg</font> </p><p><font COLOR="#FF0000">[x1,num_y1]=ode23('g1',2,4,0.5);</font> </p><p><font COLOR="#FF0000">anl_y1=x1.^3-7.5;</font> </p><p><font COLOR="#FF0000">error_1=abs(anl_y1-num_y1)./abs(anl_y1); % ode23 的计算误差</font> </p><p><font COLOR="#FF0000">[x2,num_y2]=ode45('g1',2,4,0.5);</font> </p><p><font COLOR="#FF0000">anl_y2=x2.^3-7.5; % 注意 x2 个数与 x1 不一定相同</font> </p><p><font COLOR="#FF0000">error_2=abs(anl_y2-num_y2)./abs(anl_y2); % ode45 的计算误差<br></font></p><p><font COLOR="#FF0000">hold on</font> </p><p><font COLOR="#FF0000">subplot(2,2,1)</font> </p><p><font COLOR="#FF0000">plot(x1,num_y1,x1,anl_y1,'o')</font> </p><p><font COLOR="#FF0000">title('ODE23 solution'), ylabel('y')</font> </p><p><font COLOR="#FF0000">subplot(2,2,2)</font> </p><p><font COLOR="#FF0000">plot(x1,error_y1) % 注意二种方法的误差</font> </p><p><font COLOR="#FF0000">title('ODE23 error'), ylabel('y') % ode23 的误差的量级为 1.e-16</font> </p><p><font COLOR="#FF0000">subplot(2,2,3)</font> </p><p><font COLOR="#FF0000">plot(x2,num_y2,x2,anl_y2,'o')</font> </p><p><font COLOR="#FF0000">title('ODE45 solution'), ylabel('y')</font> </p><p><font COLOR="#FF0000">subplot(2,2,4)</font> </p><p><font COLOR="#FF0000">plot(x1,error_y2)</font> </p><p><font COLOR="#FF0000">title('ODE45 error'), ylabel('y') % ode45 的解没有误差</font> </p><p><font COLOR="#FF0000">hold off<br></font></p><p>另一个例子： </p><p><font COLOR="#FF0000">% m-function, g5.m</font> </p><p><font COLOR="#FF0000">function dy=g5(x,y)</font> </p><p><font COLOR="#FF0000">dy=-y+2*cos(x);<br></font></p><p><font COLOR="#FF0000">% m-file, odes1.m</font> </p><p><font COLOR="#FF0000">% Solve an ode using ode23 and ode45</font> </p><p><font COLOR="#FF0000">clg</font> </p><p><font COLOR="#FF0000">[x1,num_y1]=ode23('g5',0,5,1);</font> </p><p><font COLOR="#FF0000">anl_y1=sin(x1)+cos(x1);</font> </p><p><font COLOR="#FF0000">error_1=abs(anl_y1-num_y1)./abs(anl_y1);</font> </p><p><font COLOR="#FF0000">[x2,num_y2]=ode45('g5',0,5,1);</font> </p><p><font COLOR="#FF0000">anl_y2=sin(x2)+cos(x2);</font> </p><p><font COLOR="#FF0000">error_2=abs(anl_y2-num_y2)./abs(anl_y2); <br></font></p><p><font COLOR="#FF0000">hold on</font> </p><p><font COLOR="#FF0000">subplot(2,2,1)</font> </p><p><font COLOR="#FF0000">plot(x1,num_y1,x1,anl_y1,'o')</font> </p><p><font COLOR="#FF0000">title('ODE23 solution'), ylabel('y')</font> </p><p><font COLOR="#FF0000">subplot(2,2,2)</font> </p><p><font COLOR="#FF0000">plot(x1,error_y1) % 注意二种方法的误差</font> </p><p><font COLOR="#FF0000">title('ODE23 error'), ylabel('y') % ode23 的误差的量级为 1.e-4</font> </p><p><font COLOR="#FF0000">subplot(2,2,3)</font> </p><p><font COLOR="#FF0000">plot(x2,num_y2,x2,anl_y2,'o')</font> </p><p><font COLOR="#FF0000">title('ODE45 solution'), ylabel('y')</font> </p><p><font COLOR="#FF0000">subplot(2,2,4)</font> </p><p><font COLOR="#FF0000">plot(x1,error_y2)</font> </p><p><font COLOR="#FF0000">title('ODE45 error'), ylabel('y') % ode45 的误差的量级为 1.e-6</font> </p><p><font COLOR="#FF0000">hold off<br></font></p><hr><p><a HREF="ch10_1.htm" tppabs="http://webclass.ncu.edu.tw/~junwu/ch10_1.htm"><img SRC="../img1/lastpage.gif" tppabs="http://webclass.ncu.edu.tw/~junwu/img/lastpage.gif" BORDER="0" WIDTH="42" HEIGHT="42"></a> <a HREF="ch10_3.htm" tppabs="http://webclass.ncu.edu.tw/~junwu/ch10_3.htm"><img SRC="../img1/nextpage.gif" tppabs="http://webclass.ncu.edu.tw/~junwu/img/nextpage.gif" BORDER="0" HSPACE="10" WIDTH="42" HEIGHT="42"></a> <a HREF="../index.html" tppabs="http://webclass.ncu.edu.tw/~junwu/index.html"><img SRC="../img1/outline.gif" tppabs="http://webclass.ncu.edu.tw/~junwu/img/outline.gif" BORDER="0" HSPACE="6" WIDTH="42" HEIGHT="42"></a> <br><font SIZE="2" COLOR="#FF00FF">上一页 下一页 讲义大纲 </font><layer src="http://www.spidersoft.com/ads/bwz468_60.htm" visibility="hidden" id="a1" width="600" onload="moveToAbsolute(ad1.pageX,ad1.pageY); a1.clip.height=60;visibility='show';"></layer> </p></body></html>⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -