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

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

h 的数值解。

在区间[0,1]上取0.1

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

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

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

本题的精确解为()

y x=

x0=0;

y0=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］用泰勒方法解

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

解：二阶泰勒方法

对于本题

故

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

x0=0;

y0=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

四阶泰勒方法

x0=0;

y0=1;

ya0=2*x0/(3*y0^2);%%一阶导数

yb0=2/(3*y0^2)-8*x0^2/(9*y0^5);%%二阶导数

yc0=-4*x0/(3*y0^5)-80*x0^3/(27*y0^8);%%三阶导数

yd0=-4/(3*y0^5)+40*x0^2/(3*y0^8)-1280*x0^4/(81*y0^11);%%四阶导数

x(1)=0.1;

y(1)=y0+0.1*ya0+0.01/2*yb0+0.001/6*yc0+0.0001/24*yd0; ya(1)=2*x(1)/(3*y(1)^2);%%一阶导数

yb(1)=2/(3*y(1)^2)-8*x(1)^2/(9*y(1)^5);%%二阶导数

yc(1)=-4*x(1)/(3*y(1)^5)-80*x(1)^3/(27*y(1)^8);%%三阶导数

yd(1)=-4/(3*y(1)^5)+40*x(1)^2/(3*y(1)^8)-1280*x(1)^4/(81*y(1)^11);%%四阶导数 for n=1:9

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

y(n+1)=y(n)+0.1*ya(n)+0.01/2*yb(n)+0.001/6*yc(n)+0.0001/24*yd(n); ya(n+1)=2*x(n+1)/(3*y(n+1)^2);%%一阶导数

yb(n+1)=2/(3*y(n+1)^2)-8*x(n+1)^2/(9*y(n+1)^5);%%二阶导数

yc(n+1)=-4*x(n+1)/(3*y(n+1)^5)-80*x(n+1)^3/(27*y(n+1)^8);%%三阶导数

yd(n+1)=-4/(3*y(n+1)^5)+40*x(n+1)^2/(3*y(n+1)^8)-1280*x(n+1)^4/(81*y(n+1)^11);%%四阶导数 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.0033 1.0132 1.0292 1.0508 1.0773 1.1081 1.1423 1.1793 Columns 9 through 10

1.2187 1.2598 ［例3］用标准四阶R －K 方法求

在区间[0, 1]上，

0.1h =的数值解以及在区间[1, 10]上，1=h 的数值解，并与精确解

进行比较。 解：对于本题

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