常微分方程数值解

第八章 常微分方程数值解

姓名 学号 班级

习题主要考察点：欧拉方法的构造，单步法的收敛性和稳定性的讨论，线性多步法中亚当姆斯方法的构造和讨论。

1 用改进的欧拉公式，求以下微分方程

]1,0[1)0(2∈⎪⎩⎪⎨⎧=-='x y y x y y

的数值解（取步长2.0=h ），并与精确解作比较。（改进的尤拉公式的应用）

解：原方程可转化为 x y y y 22

-='，令22

y z =，有x z dx dz 22-=- 解此一阶线性微分方程，可得 12+=

x y 。

利用以下公式 )4,3,2,1,0()(21)2(2.0)2(2.01=⎪⎪⎪⎩

⎪⎪⎪⎨⎧+=-⋅+=-⋅+=+i y y y y x y y y y x y y y c p i p i p i c i i i i p 求在节点)5,4,3,2,1(2.0=⋅=i i x i 处的数值解i y ，其中，初值为1,000==y x 。 MATLAB 程序如下：

x(1)=0;%初值节点

y(1)=1;%初值

fprintf('x(%d)=%f,y(%d)=%f,yy(%d)=%f\n',1,x(1),1,y(1),1,y(1)); for i=1:5

yp=y(i)+0.2*(y(i)-2*x(i)/y(i));%预报值

yc=y(i)+0.2*(yp-2*x(i)/yp);%校正值

y(i+1)=(yp+yc)/2;%改进值

x(i+1)=x(i)+0.2;%节点值

yy(i+1)=sqrt(2*x(i+1)+1);%精确解

fprintf('x(%d)=%f,y(%d)=%f,yy(%d)=%f\n',i+1,x(i+1),i+1,y(i+1),i+1,yy(i+1));

end

程序运行的结果如下：

x(1)=0.000000, y(1)=1.000000, yy(1)=1.000000

x(2)=0.200000, y(2)=1.220000, yy(2)=1.183216

x(3)=0.400000, y(3)=1.420452, yy(3)=1.341641

x(4)=0.600000, y(4)=1.615113, yy(4)=1.483240

x(5)=0.800000, y(5)=1.814224, yy(5)=1.612452

x(6)=1.000000, y(6)=2.027550, yy(6)=1.732051

2用四阶龙格－库塔法求解初值问题⎩⎨⎧==+'0

)0(1y y y ，取2.0=h , 求4.0,2.0=x 时的数值解.

要求写出由n n y x h ,,直接计算1+n y 的迭代公式，计算过程保留3位小数。（龙格－库塔方法的应用）

解：四阶龙格-库塔经典公式为

11234(22)6

n n h y y k k k k +=++++ 1(,)n n k f x y =

2111(,)22

n n k f x h y hk =++ 3211(,)22

n n k f x h y hk =++ 43(,)n n k f x h y hk =++

由于y y x f -=1),(，在各点的斜率预报值分别为：

n y k -=11

)2

1)(1()1(21)2(112h y y h y k h y k n n n n --=---=+-= )]2

1(21)[1()21)(1(21)2(123h h y h y h y k h y k n n n n ---=----=+-= ))]2

1(21(1)[1()]21(21)[1(1)(134h h h y h h y h y hk y k n n n n ----=-----=+-= 四阶经典公式可改写成以下直接的形式：

)4

36)(1(6321h h h y h y y n n n -+--+=+ 在2.01==x x 处，有

1813.0)4

)2.0()2.0(2.036)(01(62.003

21=-+⨯--+=y 在4.02==x x 处，有

3297.0)4

)2.0()2.0(2.036)(1813.01(62.01813.03

22=-+⨯--+=y 注：这两个近似值与精确解x e

y --=1在这两点的精确值十分接近。

3 用梯形方法解初值问题

⎩⎨⎧==+'1

)0(0y y y 证明其近似解为

n n h h y ⎪⎭

⎫ ⎝⎛+-=22 并证明当0→h 时，它收敛于原初值问题的准确解x e y -=。

解：显然，x e y -=是原初值问题的准确解。

求解一般微分方程初值问题的梯形公式的形式为

)],(),([2

111+++++=n n n n n n y x f y x f h y y 对于该初值问题，其梯形公式的具体形式为

)(211++--+

=n n n n y y h y y ，n n y h y h )21()21(1-=++，n n y h h y ))

22(1+-=+ 于是： 101121222222))22(++-+⎪⎭⎫ ⎝⎛+-=⎪⎭

⎫ ⎝⎛+-==⎪⎭⎫ ⎝⎛+-=+-=n n n n n h h y h h y h h y h h y 亦即：n n h h y ⎪⎭

⎫ ⎝⎛+-=22 注意到：nh nh x n =+=0，h x n n =，令h h t +-=22，2111--=t h 有 22)1()1()

1(221n n n n n

x t x x t x h x n t t t h h y ----++=+=⎪⎭⎫ ⎝⎛+-= 从而 n n

n

x x t t x t n h e t t y --→-→→=+⋅+=2000)1(lim )1(lim lim

即：当0→h 时，n y 收敛于原初值问题的准确解n x n e x y -=)(。

4对于初值问题⎩⎨⎧=-='1

)0(10y y y ，证明当2.0

证明：显式的欧拉公式为n n n n n y h y x hf y y )101(),(1-=+=+

从而n n e h e )101(1-=+，由于2.00<