📄 ch10_2.htm
字号:
<HTML>
<HEAD>
<TITLE> 阮奇-库达方法 </TITLE>
<META NAME="GENERATOR" CONTENT="Internet Assistant for Microsoft Word 2.0z">
</HEAD>
<BODY BACKGROUND="bg0000.gif" tppabs="http://166.111.167.223/computer/cai/matlabjc/img/bg0000.gif">
<H1><FONT SIZE=6 COLOR=#0000FF>10.2 阮奇-库达方法 </FONT><FONT FACE="新细明体"></FONT>
</H1>
<HR>
<P>
阮奇-库达 (Runge-Kutta) 方法是最通用的解 ODE 的方法，它可以依计算精确度的要求有低阶到高阶的各个
计算式，举例来说，一阶阮奇-库达法为
<P>
<IMG SRC="img00007-4.gif" tppabs="http://166.111.167.223/computer/cai/matlabjc/img10/img00007.gif">
<P>
其实上式即是一阶的泰勒序数近似式，只不过令 <IMG SRC="img00008-4.gif" tppabs="http://166.111.167.223/computer/cai/matlabjc/img10/img00008.gif">
。因为一阶阮奇-库达法 的精确度太低，所以在真正解ODE 时，最少须用二阶以上的方法。
<BR>
<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="img00009-4.gif" tppabs="http://166.111.167.223/computer/cai/matlabjc/img10/img00009.gif">），<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>
先前的四个 ODE 的解法如下： <BR>
<P>
例一、要在区间 [2, 4] 解以下的 ODE：
<P>
<IMG SRC="img00010-4.gif" tppabs="http://166.111.167.223/computer/cai/matlabjc/img10/img00010.gif">
<P>
<FONT COLOR=#FF0000>% m-function, g1.m</FONT>
<P>
<FONT COLOR=#FF0000>function dy=g1(x,y)</FONT>
<P>
<FONT COLOR=#FF0000>dy=3*x.^2;<BR>
</FONT>
<P>
<FONT COLOR=#FF0000>>> [x,num_y]=ode23('g1',2,4,0.5);</FONT>
<P>
<FONT COLOR=#FF0000>>> anl_y=x.^3-7.5;</FONT>
<P>
<FONT COLOR=#FF0000>>> plot(x,num_y,x,anl_y,'o')</FONT>
<P>
<FONT COLOR=#FF0000>>> title('Solution of g1')</FONT>
<P>
<FONT COLOR=#FF0000>>> xlabel('x'), ylabel('y=f(x)'), grid
<BR>
</FONT>
<P>
例二、要在区间 [0, 5] 解以下的 ODE：
<P>
<IMG SRC="img00011-3.gif" tppabs="http://166.111.167.223/computer/cai/matlabjc/img10/img00011.gif">
<P>
<FONT COLOR=#FF0000>% m-function, g2.m</FONT>
<P>
<FONT COLOR=#FF0000>function dy=g2(x,y)</FONT>
<P>
<FONT COLOR=#FF0000>dy=-0.131*y;<BR>
</FONT>
<P>
<FONT COLOR=#FF0000>>> [x,num_y]=ode23('g2',0,5,4);</FONT>
<P>
<FONT COLOR=#FF0000>>> anl_y=4*exp(-0.131*x);</FONT>
<P>
<FONT COLOR=#FF0000>>> plot(x,num_y,x,anl_y,'o')</FONT>
<P>
<FONT COLOR=#FF0000>>> title('Solution of g2')</FONT>
<P>
<FONT COLOR=#FF0000>>> xlabel('x'), ylabel('y=f(x)'), grid
<BR>
</FONT>
<P>
例三、要在区间 [0, 2] 解以下的 ODE：
<P>
<IMG SRC="img00012-4.gif" tppabs="http://166.111.167.223/computer/cai/matlabjc/img10/img00012.gif">
<P>
<FONT COLOR=#FF0000>% m-function, g3.m</FONT>
<P>
<FONT COLOR=#FF0000>function dy=g3(x,y)</FONT>
<P>
<FONT COLOR=#FF0000>dy=2*x*cos(y)^2;<BR>
</FONT>
<P>
<FONT COLOR=#FF0000>>> [x,num_y]=ode23('g3',0,2,pi/4);</FONT>
<P>
<FONT COLOR=#FF0000>>> anl_y=atan(x.*x+1);</FONT>
<P>
<FONT COLOR=#FF0000>>> plot(x,num_y,x,anl_y,'o')</FONT>
<P>
<FONT COLOR=#FF0000>>> title('Solution of g3')</FONT>
<P>
<FONT COLOR=#FF0000>>> xlabel('x'), ylabel('y=f(x)'), grid
<BR>
</FONT>
<P>
例四、要在区间 [0, 3] 解以下的 ODE：
<P>
<A NAME="work"><IMG SRC="img00013-4.gif" tppabs="http://166.111.167.223/computer/cai/matlabjc/img10/img00013.gif"></A>
<P>
<FONT COLOR=#FF0000>% m-function, g4.m</FONT>
<P>
<FONT COLOR=#FF0000>function dy=g4(x,y)</FONT>
<P>
<FONT COLOR=#FF0000>dy=3*y+exp(2*x);<BR>
</FONT>
<P>
<FONT COLOR=#FF0000>>> [x,num_y]=ode23('g4',0,3,3);</FONT>
<P>
<FONT COLOR=#FF0000>>> anl_y=4*exp(3*x)-exp(2*x);</FONT>
<P>
<FONT COLOR=#FF0000>>> plot(x,num_y,x,anl_y,'o')</FONT>
<P>
<FONT COLOR=#FF0000>>> title('Solution of g4')</FONT>
<P>
<FONT COLOR=#FF0000>>> xlabel('x'), ylabel('y=f(x)'), grid
<BR>
</FONT>
<P>
如果将上述方法改以 <FONT COLOR=#FF0000>ode45</FONT> 计算，你可能无法察觉出其与<FONT COLOR=#FF0000>ode23</FONT>的解之间的差异，那是因为我们选的
ODE 代表的函数分布变化平缓，所以高阶方法就显现不出其优点。不过以二者在计算的误差上做比较，<FONT COLOR=#FF0000>ode45</FONT>
的误差量级会比 <FONT COLOR=#FF0000>ode23</FONT>要小。以下是几个例子：
<P>
<FONT COLOR=#FF0000>% m-function, g1.m</FONT>
<P>
<FONT COLOR=#FF0000>function dy=g1(x,y)</FONT>
<P>
<FONT COLOR=#FF0000>dy=3*x.^2;<BR>
</FONT>
<P>
<FONT COLOR=#FF0000>% m-file, odes1.m</FONT>
<P>
<FONT COLOR=#FF0000>% Solve an ode using ode23 and ode45</FONT>
<P>
<FONT COLOR=#FF0000>clg</FONT>
<P>
<FONT COLOR=#FF0000>[x1,num_y1]=ode23('g1',2,4,0.5);</FONT>
<P>
<FONT COLOR=#FF0000>anl_y1=x1.^3-7.5;</FONT>
<P>
<FONT COLOR=#FF0000>error_1=abs(anl_y1-num_y1)./abs(anl_y1); %
ode23 的计算误差</FONT>
<P>
<FONT COLOR=#FF0000>[x2,num_y2]=ode45('g1',2,4,0.5);</FONT>
<P>
<FONT COLOR=#FF0000>anl_y2=x2.^3-7.5; % 注意 x2 个数与 x1 不一定相同</FONT>
<P>
<FONT COLOR=#FF0000>error_2=abs(anl_y2-num_y2)./abs(anl_y2); %
ode45 的计算误差<BR>
</FONT>
<P>
<FONT COLOR=#FF0000>hold on</FONT>
<P>
<FONT COLOR=#FF0000>subplot(2,2,1)</FONT>
<P>
<FONT COLOR=#FF0000>plot(x1,num_y1,x1,anl_y1,'o')</FONT>
<P>
<FONT COLOR=#FF0000>title('ODE23 solution'), ylabel('y')</FONT>
<P>
<FONT COLOR=#FF0000>subplot(2,2,2)</FONT>
<P>
<FONT COLOR=#FF0000>plot(x1,error_y1) % 注意二种方法的误差</FONT>
<P>
<FONT COLOR=#FF0000>title('ODE23 error'), ylabel('y') % ode23
的误差的量级为 1.e-16</FONT>
<P>
<FONT COLOR=#FF0000>subplot(2,2,3)</FONT>
<P>
<FONT COLOR=#FF0000>plot(x2,num_y2,x2,anl_y2,'o')</FONT>
<P>
<FONT COLOR=#FF0000>title('ODE45 solution'), ylabel('y')</FONT>
<P>
<FONT COLOR=#FF0000>subplot(2,2,4)</FONT>
<P>
<FONT COLOR=#FF0000>plot(x1,error_y2)</FONT>
<P>
<FONT COLOR=#FF0000>title('ODE45 error'), ylabel('y') % ode45
的解没有误差</FONT>
<P>
<FONT COLOR=#FF0000>hold off<BR>
</FONT>
<P>
另一个例子：
<P>
<FONT COLOR=#FF0000>% m-function, g5.m</FONT>
<P>
<FONT COLOR=#FF0000>function dy=g5(x,y)</FONT>
<P>
<FONT COLOR=#FF0000>dy=-y+2*cos(x);<BR>
</FONT>
<P>
<FONT COLOR=#FF0000>% m-file, odes1.m</FONT>
<P>
<FONT COLOR=#FF0000>% Solve an ode using ode23 and ode45</FONT>
<P>
<FONT COLOR=#FF0000>clg</FONT>
<P>
<FONT COLOR=#FF0000>[x1,num_y1]=ode23('g5',0,5,1);</FONT>
<P>
<FONT COLOR=#FF0000>anl_y1=sin(x1)+cos(x1);</FONT>
<P>
<FONT COLOR=#FF0000>error_1=abs(anl_y1-num_y1)./abs(anl_y1);</FONT>
<P>
<FONT COLOR=#FF0000>[x2,num_y2]=ode45('g5',0,5,1);</FONT>
<P>
<FONT COLOR=#FF0000>anl_y2=sin(x2)+cos(x2);</FONT>
<P>
<FONT COLOR=#FF0000>error_2=abs(anl_y2-num_y2)./abs(anl_y2);
<BR>
</FONT>
<P>
<FONT COLOR=#FF0000>hold on</FONT>
<P>
<FONT COLOR=#FF0000>subplot(2,2,1)</FONT>
<P>
<FONT COLOR=#FF0000>plot(x1,num_y1,x1,anl_y1,'o')</FONT>
<P>
<FONT COLOR=#FF0000>title('ODE23 solution'), ylabel('y')</FONT>
<P>
<FONT COLOR=#FF0000>subplot(2,2,2)</FONT>
<P>
<FONT COLOR=#FF0000>plot(x1,error_y1) % 注意二种方法的误差</FONT>
<P>
<FONT COLOR=#FF0000>title('ODE23 error'), ylabel('y') % ode23
的误差的量级为 1.e-4</FONT>
<P>
<FONT COLOR=#FF0000>subplot(2,2,3)</FONT>
<P>
<FONT COLOR=#FF0000>plot(x2,num_y2,x2,anl_y2,'o')</FONT>
<P>
<FONT COLOR=#FF0000>title('ODE45 solution'), ylabel('y')</FONT>
<P>
<FONT COLOR=#FF0000>subplot(2,2,4)</FONT>
<P>
<FONT COLOR=#FF0000>plot(x1,error_y2)</FONT>
<P>
<FONT COLOR=#FF0000>title('ODE45 error'), ylabel('y') % ode45
的误差的量级为 1.e-6</FONT>
<P>
<FONT COLOR=#FF0000>hold off<BR>
</FONT>
<HR>
<P>
<A HREF="ch10_1.htm" tppabs="http://166.111.167.223/computer/cai/matlabjc/ch10_1.htm"><IMG SRC="lastpage.gif" tppabs="http://166.111.167.223/computer/cai/matlabjc/img/lastpage.gif" BORDER=0></A>
<A HREF="ch10_3.htm" tppabs="http://166.111.167.223/computer/cai/matlabjc/ch10_3.htm"><IMG SRC="nextpage-1.gif" tppabs="http://166.111.167.223/computer/cai/matlabjc/img/nextpage.gif" BORDER=0 HSPACE=10></A>
<A HREF="index.html" tppabs="http://166.111.167.223/computer/cai/matlabjc/index.html"><IMG SRC="outline-1.gif" tppabs="http://166.111.167.223/computer/cai/matlabjc/img/outline.gif" BORDER=0 HSPACE=6></A>
<BR>
<FONT SIZE=2 COLOR=#FF00FF>上一页 下一页 讲义大纲 </FONT>
</BODY>
</HTML>
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -