常微分方程初值问题的数值解法

贵州师范大学数学与计算机科学学院学生实验报告

课程名称： 数值分析 班级： 实验日期： 年 月 日 学 号： 姓名： 指导教师： 实验成绩： 一、实验名称

实验六: 常微分方程初值问题数值解法 二、实验目的及要求

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))