试验四：四阶 Runge-Kutta 方法 一、 试验目的 在计算机上实现用四阶 Runge-Kutta 求一阶常微分方程初值问题   ' xy    ay  ,   , xyxf 1 y     , bz  ,  的数值解，并利用最后绘制的图形直观分析近似解与准确解之间的误差。 二、 计算公式或算法 本算法是根据四阶 R-K 法而设计，它是求一阶常微分方程初值问题    ' y    ay ,   , xyxf y  0  , ba   , ; 的计算解，即求出解  xy 在节点 ix 上的计算值  iyi yxf （编写或调用计算   yxf , , ynba , , , 0  ,2,1 ,  n . , 的函数文件  yxF , ）， % x 0  , xa n  b  i   1 h 0 2  xf i  xf  xf i  xf 1  i 1  1  i 1  ,  y 1 i  , yh  1 , yh  1 , yh  Kh  6  1 i   1 Kh  1 Kh  1 3 hK  2 i 1  i 1  1  2 K 2  2 K 3  K 4 1．输入 2．   h ab   y  0 3．For i n  0 xy :1 n i x x  h h 1 K K K K 1       3 2 1 4 y i y   1 i  End 4．输出 , yy 1 2  ny . 三、 Matlab 程序 x=[a:h:b]; y(1)=y1; n=(b-a)/h+1; for i=2:n 

fk1=f(x(i-1),y(i-1)); fk2=f(x(i-1)+h/2,y(i-1)+fk1*h/2); fk3=f(x(i-1)+h/2,y(i-1)+fk2*h/2); fk4=f(x(i-1)+h,y(i-1)+fk3*h); y(i)=y(i-1)+h*(fk1+2*fk2+2*fk3+fk4)/6; end y 四、 测试数据及结果 （2）用调试好的程序解决如下问题： 求初值问题 (  y 2 x  4 x  ),1 x  ],5.0,0[    ' )( xy  1 2 )0( 0 y  的解 y(x)在 xi  （h=0.05）处的近似值 iy ,并与问题的解析解 ih )( xy   e x 2  2 x  1 相比较。 步骤一：编写一函数子程序. 编写程序如下： %f.m function z=f(x,y) z=(-y+x^2+4*x-1)/2; end; 步骤二：执行上述 Runge-Kutta 算法，计算结果为 ix iy )iy x ( ix iy )iy x ( 0.05 0.10 0.15 0.20 0.25 0.022190087 -0.038770574 -0.049756511 -0.055162579 -0.055003093 -0.022190087 -0.038770575 -0.049756514 -0.055162582 -0.055003097 0.30 0.35 0.40 0.45 0.50 -0.049292019 -0.038042973 -0.021269240 0.0010612260 0.028800791 -0.04929024 -0.038042979 -0.021269274 0.0010162188 0.028800783 

（3）使用 Matlab 绘图函数“plot(x,y)”绘制（2）的数值解和解析解的图形。 数值解的图形： >>runge >> x=0:0.05:0.5; >> plot(x,y) 

解 析 解 >> x=0:0.05:0.5; >> y=exp(-x./2)+x.^2-1; >> plot(x,y) 0.03 0.02 0.01 0 -0.01 -0.02 -0.03 -0.04 -0.05 y -0.06 0 0.05 0.1 0.15 0.2 0.3 0.35 0.4 0.45 0.5 0.25 x （4）使用 Matlab 中的 ode45 求解（2），并绘图。 编写函数如下： %ode.m function dy=ode(x,y) dy=[0.5*(-y+x^2+4*x-1)]; >> [T,Y]=ode45('ode',[0 0.5],0); >>plot(T,Y) 运行结果如下： 

0.03 0.02 0.01 0 -0.01 -0.02 -0.03 -0.04 -0.05 -0.06 0 ode45求 解 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 五、 心得体会 学会用 Matlab 绘图来分析误差。