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

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

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

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

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

出师表

两汉：诸葛亮

先帝创业未半而中道崩殂，今天下三分，益州疲弊，此诚危急存亡之秋也。然侍卫之臣不懈于内，忠志之士忘身于外者，盖追先帝之殊遇，欲报之于陛下也。诚宜开张圣听，以光先帝遗德，恢弘志士之气，不宜妄自菲薄，引喻失义，以塞忠谏之路也。

宫中府中，俱为一体；陟罚臧否，不宜异同。若有作奸犯科及为忠善者，宜付有司论其刑赏，以昭陛下平明之理；不宜偏私，使内外异法也。

侍中、侍郎郭攸之、费祎、董允等，此皆良实，志虑忠纯，是以先帝简拔以遗陛下：愚以为宫中之事，事无大小，悉以咨之，然后施行，必能裨补阙漏，有所广益。

将军向宠，性行淑均，晓畅军事，试用于昔日，先帝称之曰“能”，是以众议举宠为督：愚以为营中之事，悉以咨之，必能使行阵和睦，优劣得所。