数学实验“线性多步法(数值积分法,Taylor展开法)”实验报告(内含matlab程序)

西京学院数学软件实验任务书

实验二十五实验报告

一、实验名称：线性多步法（数值积分法，Taylor 展开法）。 二、实验目的：进一步熟悉线性多步法（数值积分法，Taylor 展开法）。

三、实验要求：运用Matlab/C/C++/Java/Maple/Mathematica 等其中一种语言完成程序设计。 四、实验原理：

1．数值积分法： 常微分方程初值问题

000

(,),

()()dy

f x y x x dx

y x y ⎧=⎪≤⎨⎪=⎩ (1) 数值解法中，某一步的公式不仅与前一步解的值有关，而且与前若干步解的值有关，利用前面多步的信息预测下一步的值，这就是多步法的基本思想，可以期望获得较高的精度。

将(1)中的方程在区间[]1n n x x +，上积分，可以得到：

()()()()11n n

x n n x y x y x f x y x dx ++=+⎰

，。

用等距节点的插值多项式来替代被积函数，再对插值多项式积分，这样就得到一系列求积公式。

用梯形方法计算积分项

()()()()()()1

112

n n

x n n n n x h f x y x dx f x y x f x y x +++⎡⎤≈

+⎣⎦⎰

，，， 代入(1)中得：

()()()()()()1112

n n n n n n h y x y x f x y x f x y x +++⎡⎤≈+

+⎣⎦，， 设由1r +个数据点()()()11n n n n n r n r x f x f x f ---- ，，，，，， 构造插值多项式()r P x ，这里()0k k k k f f x y x x kh ==+，，，运用插值公式有：

()()()00r

r

n k

r n j j j j k k j n j n k

x x P x f l x l x x x --==≠---==

-∑∏，， 得到下列计算公式：

10r

n n rj n j j y y h f α+-==+∑ (2)

其中，()1100101n n r x rj j x k k j

t k

l x dx dt j r h k j α+=≠+===-∏⎰⎰ ，，，，

由此可得(2)中的系数。公式(2)是一个r+1步的显式公式，称为Adams 显式公式。

2．Taylor 展开法：

基于数值积分可以构造出一系列求解常微分方程的计算公式，下面介绍基于 Taylor 展开的待定系数法，它可灵活地构造出线性多步法。对固定的系数，可以选取待定系数使线性

多步法的阶尽可能高。还可以根据需要，确定显式还是隐式。

设构造如下具有 p 阶精度的线性多步公式：

()1011110n n n r n r n n r n r y y y y h f f f αααβββ+---+-=+++++++

(4)

它们的局部截断误差为:

()()()()1101r r

n n k n k k n k n k k k T y x y x h f x y x αβ++---==-⎡⎤

=-+⎢⎥⎣⎦∑∑，

即：()()()1101r r

n n k n k k n k k k T y x y x h y x αβ++--==-⎡⎤

'=-+⎢⎥⎣⎦

∑∑ (5)

对(5)式的右端各项在n x 点处作Taylor 展开有:

()()

()()()()()()()()()()()()()()()()1

120

11

11!1!1!

!j p p

j p p n k n n j j p p

j p p n k n n j kh kh y x y x y x O h j p kh kh y x y x y x O h j p +++-=+++-=--=+++--'=++-∑∑，

。

代入(5)中得：

()()

()()()1011

1

1

111

121

!

(1)[1()! ()][1()1! (1)()]()()

j

p

r

r n k n k

k j k j p p r

r

j k n k k k r

p p p k n k h T y x k j h j k y x k p p k y x O h ααβαβ+===-++==++==-+----+--+-+-+∑∑∑∑∑∑

使()2p n y x h h h ，，，，的系数为零，得到关于k α和k β的线性方程组：

()()0

11

11 112r

k k j j r r

k k k k k j k j p ααβ=-==-⎧=⎪⎪⎨⎪-+-==⎪⎩∑∑∑ ，，，，。 而且得到线性多步法的局部截断误差：

()()()()()()()1

1

11

121

[11! 1]p p r

n k

k p

r

p p k n k h T k p p k y x O h αβ+++=++=-=--+-+-+∑∑。

五、实验内容：

%数值积分法

function [k,X,Y,wucha,P]= Adams4x(funfcn,x0,b,y0,h) x=x0; y=y0;p=128; n=fix((b-x0)/h); if n<5

return ; end

X=zeros(p,1);

Y=zeros(p,length(y)); f=zeros(p,1); k=1; X(k)=x; Y(k,:)=y'; for k=2:4

c1=1/6;c2=2/6;c3=2/6; c4=1/6;a2=1/2; a3=1/2; a4=1;b21=1/2;b31=0;

b32=1/2; b41=0;b42=0;b43=1; x1=x+a2*h; x2=x+a3*h;

x3=x+a4*h; k1=feval(funfcn,x,y); y1=y+b21*h*k1; x=x+h; k2=feval(funfcn,x1,y1); y2=y+b31*h*k1+b32*h*k2; k3=feval(funfcn,x2,y2);

y3=y+b41*h*k1+b42*h*k2+b43*h*k3; k4=feval(funfcn,x3,y3);

y=y+h*(c1*k1+c2*k2+c3*k3+c4*k4); X(k)=x; Y(k,:)=y;