常微分方程的数值解法(欧拉法、改进欧拉法、泰勒方法和龙格-库塔法)

［例1］用欧拉方法与改进的欧拉方法求初值问题

223(0)1dy x dx y y ⎧=⎪⎨

⎪=⎩

在区间[0,1]上取0.1h =的数值解。

［解］ 欧拉方法的计算公式为

102

2(

)1,0.13n

n n n

x y y h y h y +=+==

现用matlab 编程，程序如下

x0=0; y0=1; x(1)=0.1;

y(1)=y0+0.1*2*x0/(3*y0^2); for n=1:9

x(n+1)=0.1*(n+1);

y(n+1)=y(n)+0.1*2*x(n)/(3*y(n)^2); end ; x y 结果为 x =

Columns 1 through 8

0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000

Columns 9 through 10

0.9000 1.0000 y =

Columns 1 through 8

1.0000 1.0067 1.0198 1.0391 1.0638 1.0932 1.1267 1.1634

Columns 9 through 10

1.2028 1.2443 改进的欧拉方法其计算公式为

(0)

1

21

12(0)2

102()322[]233()1,0.1n n n n

n n n n n n x y y h y x x h y y y y y h ++++⎧=+⎪⎪⎪=++⎨⎪

⎪==⎪⎩

本题的精确解为()y x = 使用excel 表格进行运算，相应如下

现用matlab编程，程序如下

x0=0;

y0=1;

x(1)=0.1;

ya(1)=y0+0.1*2*x0/(3*y0^2);

y(1)=y0+0.05*(2*x0/(3*y0^2)+2*x0/(3*ya^2));

for n=1:9

x(n+1)=0.1*(n+1);

ya(n+1)=ya(n)+0.1*2*x(n)/(3*ya(n)^2);

y(n+1)=y(n)+0.05*(2*x(n)/(3*y(n)^2)+2*x(n+1)/(3*ya(n+1)^2));

end;

x

y

结果为

x =

Columns 1 through 8

0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000

Columns 9 through 10

0.9000 1.0000

y =

Columns 1 through 8

1.0000 1.0099 1.0261 1.0479 1.0748 1.1059 1.1407 1.1783

Columns 9 through 10

1.2183 1.2600

［例2］用泰勒方法解

223(0)1dy x dx y y ⎧=⎪⎨⎪=⎩

分别用二阶、四阶泰勒方法计算点n x ＝0.1, 0.2, …, 1.0处的数值解，并与精确解进行

比较。

解：二阶泰勒方法

()1

,2n n n n x y h f f y y h f f x y +⎡⎤

⎛⎫∂∂=+++⎢⎥ ⎪∂∂⎝⎭⎣

⎦

对于本题

2

23x

f y

=

223223224=33324133f f x x f x y y y y x y y ⎛⎫∂∂++- ⎪∂∂⎝⎭

⎛⎫

=- ⎪⎝⎭

故

()1

,2

2

3242133n n n n x y n n n n

n h f f y y h f f x y x h

y x h y y +⎡⎤

⎛⎫∂∂=+++⎢⎥

⎪∂∂⎝⎭⎣⎦⎡⎤⎛⎫=++-⎢⎥

⎪⎢⎥⎝⎭⎣⎦

使用excel 表格进行运算，相应结果如下

现用matlab编程，程序如下

x0=0;

y0=1;

x(1)=0.1;

y(1)=y0+0.1/(3*y0^2)*(2*x0+0.1*(1-4*x0^2/(3*y0^3)));

for n=1:9

x(n+1)=0.1*(n+1);

y(n+1)=y(n)+0.1/(3*y(n)^2)*(2*x(n)+0.1*(1-4*x(n)^2/(3*y(n)^3)));

end;

x

y

结果为

x =

Columns 1 through 9

0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000

Column 10

1.0000

y =

Columns 1 through 9

1.0033 1.0132 1.0293 1.0510 1.0776 1.1084 1.1427 1.1799 1.2193

Column 10

1.2606

四阶泰勒方法

()()()()223232

2

25

253625362

358

2358692,324,33242333283916440,93944402939348032748020640,39327x f x y y x f x y y y y x x y y y x y y x x f x y y y y y x x x y y y y x x y y x x x f x y y y y y y =⎛⎫''

=+- ⎪⎝⎭⎛⎫=+- ⎪⎝⎭=-⎛⎫'''

=-+-+ ⎪⎝⎭⎛⎫=-+-+ ⎪⎝⎭=-+⎛⎫

''''=-++- ⎪⎝⎭23

58692

24581148020640239327344012803381x x x x y y y y y x x y y y ⎧

⎪⎪⎪⎪⎪

⎪⎪⎪

⎪⎪⎪⎪

⎪

⎪⎪

⎨

⎪⎪⎪

⎪

⎪⎪⎪⎪⎪

⎪⎛⎫⎪=-++- ⎪⎪⎝⎭⎪

⎪

=-+-⎪⎩

()()()()2

2

25

358

(4)

24

58112,3,82

39,480 327,40128043381n n n n n n n n n n n n n n n n

n n n n n n n n n n x y f x y y y f x y x y y y f x y x x y y y f x y x x y y y ⎧'==⎪⎪

⎪'''=⎪

⎪=-⎪⎪

⎪'''

''=⎨⎪⎪=-+⎪⎪'''=⎪⎪⎪=-+-⎪⎩