贵州师范大学数学与计算机科学学院学生实验报告
课程名称： 数值分析 班级： 实验日期： 年 月 日 学 号： 姓名： 指导教师： 实验成绩： 一、实验名称
实验六: 常微分方程初值问题数值解法 二、实验目的及要求
1. 让学生掌握用Euler法, Runge-Kutta法求解常微分方程初值问题. 2. 培养Matlab编程与上机调试能力. 三、实验环境
每人一台计算机，要求安装Windows XP操作系统，Microsoft office2003、MATLAB6.5(或7.0). 四、实验内容
1. 取步长h=0.1,0.05,0.01, ? ,用Euler法及经典4阶Runge-Kutta 法求解初值
问题
?y'??2y?2t2?2t(0?t?1) ??y(0)?1要求:
1) 画出准确解(准确解y?e?2t?t2)的曲线,近似解折线;
2) 把节点0.1和0.5上的精确解与近似解比较,观察误差变化情况.
2. 用 Euler法,隐式Euler法和经典4阶R-K法取不同步长解初值问题
?y'??50y,x?[0,1],? ?1?y(0)?2?并画出曲线观察稳定性. 注：题1必须写实验报告
五、算法描述及实验步骤
Euler法：
输入 f(x,y),a,b,h,x0(x0?a),y0 输出 Euler解y 步1 m?b?a;xn?a?n?h(n?1,2,?,m) h步2 对n?0,1,2,?,m?1执行yn?1?yn?h?f(xn,yn)
步3 输出y?(y1,y2,?,ym)T 经典4阶R-K法：
输入 f(x,y),a,b,h,x0(x0?a),y0 输出 4阶R-K解y 步1 m?b?a;xn?a?n?h(n?1,2,?,m) h步2 对n?0,1,2,?,m?1执行K1?f(xn,yn)，K2?f(xn?0.5,yn?0.5hK1)， K3?f(xn?0.5,yn?0.5hK2)，K4?f(xn?1,yn?hK3) yn?1?yn?步3 输出y?(y1,y2,?,ym)T
h(K1?2K2?2K3?K4) 6
六、调试过程及实验结果
>> shiyan6 Y1 =
0.8000 0.6620 0.5776 0.5401 0.5441 0.5853 0.6602 0.7662 0.9009 1.0627 Y2 =
0.8287 0.7103 0.6388 0.6093 0.6179 0.6612 0.7366 0.8419 0.9753 1.1353
e1 = 0.0287 e2 =
4.2469e-006
e1 = 0.0738 e2 =
1.1609e-005
注：至于h=0.05、0.01的情况将程序中的h值作相应的改动即可得。
七、总结
在1阶的Euler法解初值问题时，若步长过大的话，误差将会较大，其解不可靠，只有
控制步长，在一定误差范围内才可用。
就精度阶而言经典4阶Runge-Kutta法最可取，但并不是阶高的方法就一定可用，还必须考虑初值问题的性态、数值方法的计算量稳定性等因素。因此在实际计算中，应根据问题的具体情况来选择合适的方法。
八、附录（源程序清单）
图像、计算、误差比较程序：
x=0:0.001:1; y=exp(-2*x)+x.^2; plot(x,y,'r') axis([0,1,0.5,1.2]) hold on
a=0;b=1;y0=1;h=0.1;d=y0; Y1=Euler('fun',a,b,y0,h) u1=0:0.001:h; v1=Y1(1)+0*u1; plot(u1,v1,'g--') hold on
u2=0:0.001:5*h; v2=Y1(5)+0*u2; plot(u2,v2,'g--') hold on Y1=[d,Y1]; t=0:h:1;
scatter(t,Y1,'r') hold on plot(t,Y1) hold on
Y2=RK('fun',a,b,y0,h) u3=0:0.001:h; v3=Y2(1)+0*u3; plot(u3,v3,'y--') hold on
u4=0:0.001:5*h;
v4=Y2(5)+0*u4; plot(u4,v4,'y--') hold on
v5=0:0.001:Y2(1); u5=h+0*v5;
plot(u5,v5,'k--') hold on
v6=0:0.001:Y2(5); u6=5*h+0*v6; plot(u6,v6,'k--') hold on Y2=[d,Y2];
scatter(t,Y2,'r') hold on plot(t,Y2)
title('??è·?a?ú??ó??ü???a????','fontsize',10,'fontweight','bold')
text(0.735,0.7,'\\leftarrowy=Euler·¨?ü???a????','fontsize',8)
text(0.15,0.775,'\\leftarrowy=?-μ?4?×Runge-Kutta·¨?ü???a????','fontsize',8) x=[0.1,0.5]; y=exp(-2*x)+x.^2; e1=abs(y(1)-Y1(2)) e2=abs(y(1)-Y2(2)) e1=abs(y(2)-Y1(6)) e2=abs(y(2)-Y2(6))
Euler法程序：
function Y=Euler(f,a,b,y0,h) m=(b-a)/h;Y=zeros(1,m);x=a;d=y0; for n=1:m
K=feval(f,x,y0); x=x+h; y0=y0+h*K; Y(n)=y0; end
定义函数：
function z=fun(t,y) z=-2*y+2*t^2+2*t;
经典4阶Runge—Kutta法程序：
function Y=RK(f,a,b,y0,h) m=(b-a)/h;Y=zeros(1,m);x=a;d=y0; for n=1:m
K1=feval(f,x,y0); x=x+0.5*h;y1=y0+0.5*h*K1; K2=feval(f,x,y1); y2=y0+0.5*h*K2; K3=feval(f,x,y2); x=x+0.5*h;y3=y0+h*K3; K4=feval(f,x,y3);
y0=y0+h*(K1+2*K2+2*K3+K4)/6; Y(n)=y0; end