贵州师范大学数学与计算机科学学院学生实验报告
课程名称： 数值分析 班级： 实验日期： 年 月 日 学 号： 姓名： 指导教师： 实验成绩： 一、实验名称
实验六: 常微分方程初值问题数值解法 二、实验目的及要求
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 法求解初值
问题
⎩⎨
⎧=≤≤++-=1
)0()
10(2222'y t t t y y 要求:
1) 画出准确解(准确解22t e y t +=-)的曲线,近似解折线;
2) 把节点0.1和0.5上的精确解与近似解比较,观察误差变化情况.
2. 用 Euler 法,隐式Euler 法和经典4阶R-K 法取不同步长解初值问题
⎪⎩
⎪
⎨⎧=
∈-=21
)0(],1,0[,50'y x y y 并画出曲线观察稳定性. 注：题1必须写实验报告
五、算法描述及实验步骤
Euler 法：
输入 000),(,,,),,(y a x x h b a y x f = 输出 Euler 解y
步1 ),,2,1(;m n h n a x h
a
b m n =⨯+=-⇐
步2 对1,,2,1,0-=m n 执行),(1n n n n y x f h y y ⨯+⇐+ 步3 输出T m y y y y ),,,(21 = 经典4阶R-K 法：
输入 000),(,,,),,(y a x x h b a y x f = 输出 4阶R-K 解y 步1 ),,2,1(;m n h n a x h
a
b m n =⨯+=-⇐
步2 对1,,2,1,0-=m n 执行),(1n n y x f K ⇐，)5.0,(15.02hK y x f K n n +⇐+， )5.0,(25.03hK y x f K n n +⇐+，),(314hK y x f K n n +⇐+ )22(6
43211K K K K h
y y n n ++++⇐+ 步3 输出T m y y y y ),,,(21 =
六、调试过程及实验结果
>> 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('¾«È·½âÇúÏßÓë½üËƽâÕÛÏß','fonts ize',10,'fontweight','bold')
text(0.735,0.7,'\leftarrowy=Euler·¨½üËƽâÕÛÏß','fontsize',8)
text(0.15,0.775,'\leftarrowy=¾-µä4½×R unge-Kutta·¨½üËƽâÕÛÏß','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
出师表
两汉：诸葛亮
先帝创业未半而中道崩殂，今天下三分，益州疲弊，此诚危急存亡之秋也。然侍卫之臣不懈于内，忠志之士忘身于外者，盖追先帝之殊遇，欲报之于陛下也。诚宜开张圣听，以光先帝遗德，恢弘志士之气，不宜妄自菲薄，引喻失义，以塞忠谏之路也。
宫中府中，俱为一体；陟罚臧否，不宜异同。若有作奸犯科及为忠善者，宜付有司论其刑赏，以昭陛下平明之理；不宜偏私，使内外异法也。
侍中、侍郎郭攸之、费祎、董允等，此皆良实，志虑忠纯，是以先帝简拔以遗陛下：愚以为宫中之事，事无大小，悉以咨之，然后施行，必能裨补阙漏，有所广益。
将军向宠，性行淑均，晓畅军事，试用于昔日，先帝称之曰“能”，是以众议举宠为督：愚以为营中之事，悉以咨之，必能使行阵和睦，优劣得所。
