四阶龙格库塔和向前欧拉方法和隐式欧拉法以及改进欧拉法

实验九欧拉方法

信息与计算科学金融崔振威201002034031

一、实验目的：

1、掌握欧拉算法设计及程序实现

二、实验内容：

1、p364-9.2.4、p386-9.5.6

三、实验要求：

主程序：

欧拉方法（前项）：

function [x,y]=DEEuler(f,a,b,y0,n);

%f:一阶常微分方程的一般表达式的右端函数

%a:自变量的取值下限

%b:自变量的取值上限

%y0：函数的初值

%n:积分的步数

if nargin<5,n=50;

end

h=(b-a)/n;

x(1)=a;y(1)=y0;

for i=1:n

x(i+1)=x(i)+h;

y(i+1)=y(i)+h*feval(f,x(i),y(i));

end

format short

欧拉方法（后项）：

function [x,y]=BAEuler(f,a,b,y0,n);

%f:一阶常微分方程的一般表达式的右端函数

%a:自变量的取值下限

%b:自变量的取值上限

%y0：函数的初值

%n:积分的步数

if nargin<5,n=50;

end

h=(b-a)/n;

x(1)=a;y(1)=y0;

for i=1:n

x(i+1)=x(i)+h;

y1=y(i)+h*feval(f,x(i),y(i));

y(i+1)=y(i)+h*feval(f,x(i+1),y1);

end

format short

梯形算法：

function [I,step,h2] = CombineTraprl(f,a,b,eps)

%f 被积函数

%a,b 积分上下限

%eps 精度

%I 积分结果

%step 积分的子区间数

if(nargin ==3)

eps=1.0e-4;

end

n=1;

h=(b-a)/2;

I1=0;

I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h;

while abs(I2-I1)>eps

n=n+1;

h=(b-a)/n;

I1=I2;

I2=0;

for i=0:n-1

x=a+h*i;

x1=x+h;

I2=I2+(h/2)*(subs(sym(f),findsym(sym(f)),x)+subs(sym(f),findsym(sym(f)),x1));

end

end

I=I2;

step=n;

h2=(b-a)/n;

改进欧拉方法：

function [x,y]=MoEuler(f,a,b,y0,n);

%f:一阶常微分方程的一般表达式的右端函数

%a:自变量的取值下限

%b:自变量的取值上限

%y0：函数的初值

%n:积分的步数

if nargin<5,n=50;

end

h=(b-a)/n;

x(1)=a;y(1)=y0; for i=1:n

x(i+1)=x(i)+h;

y1=y(i)+h*feval(f,x(i),y(i)); y2=y(i)+h*feval(f,x(i+1),y1); y(i+1)=(y1+y2)/2; end

format short

四阶龙格-库塔法：

function y = DELGKT4_lungkuta(f, h,a,b,y0,varvec) %f:一阶常微分方程的一般表达式的右端函数 %h:积分步长

%a ：自变量的取值下限 %b:自变量的取值上限 %varvec ：常微分方程的变量组 format long; N = (b-a)/h; y = zeros(N+1,1); y(1) = y0;

x = a:h:b;

var = findsym(f); for i=2:N+1

K1 = Funval(f,varvec,[x(i-1) y(i-1)]);%Funval 为程序所需要的函数 K2 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]); K3 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K2*h/2]); K4 = Funval(f,varvec,[x(i-1)+h y(i-1)+h*K3]); y(i) = y(i-1)+h*(K1+2*K2+2*K3+K4)/6; end format short;

p364-9.2.4

欧拉方法（前项）：

1、22)(,1)0(,2

2

'

+-+-==-=-t t e t y y y t y t

解:

执行20步时：

编写函数文件doty.m 如下： function f=doty(x,y); f=x^2-y;

在Matlab 命令窗口输入： >> [x1,y1]=DEEuler('doty',0,2,1,20)

可得到结果：

x1 =

0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000

0.8000 0.9000 1.0000 1.1000 1.2000 1.3000 1.4000 1.5000

1.6000 1.7000 1.8000 1.9000

2.0000

y1 =

1.0000 0.9000 0.8110 0.7339 0.6695 0.6186 0.5817 0.5595

0.5526 0.5613 0.5862 0.6276 0.6858 0.7612 0.8541 0.9647

1.0932 1.2399 1.4049 1.5884 1.7906

在Matlab 命令窗口输入：

>> y3=-exp(-x1)+x1.^2-2*x1+2

求得解析解：

y3 =

1.0000 0.9052 0.8213 0.7492 0.6897 0.6435 0.6112 0.5934

0.5907 0.6034 0.6321 0.6771 0.7388 0.8175 0.9134 1.0269

1.1581 1.3073 1.4747 1.6604 1.8647

输入：

>> plot(x1,y1,'o');hold on

>> plot(x1,y3,'*');hold on

可得近似值与解析解的图像比较：

执行40步时：

在Matlab 命令窗口输入：

>> [x1,y1]=DEEuler('doty',0,2,1,40)

可得到结果：

x1 =

0 0.0500 0.1000 0.1500 0.2000 0.2500 0.3000 0.3500 0.4000 0.4500 0.5000 0.5500 0.6000 0.6500 0.7000 0.7500

0.8000 0.8500 0.9000 0.9500 1.0000 1.0500 1.1000 1.1500

1.2000 1.2500 1.3000 1.3500 1.4000 1.4500 1.5000 1.5500

1.6000 1.6500 1.7000 1.7500 1.8000 1.8500 1.9000 1.9500

2.0000

y1 =

1.0000 0.9500 0.9026 0.8580 0.8162 0.7774 0.7417 0.7091 0.6798 0.6538 0.6312 0.6121 0.5967 0.5848 0.5767 0.5724 0.5719 0.5753 0.5826 0.5940 0.6094 0.6290 0.6526 0.6805

0.7126 0.7490 0.7897 0.8347 0.8841 0.9379 0.9961 1.0588

1.1260 1.1977 1.2739 1.3547 1.4401 1.5301 1.6247 1.7240 1.8279

在Matlab 命令窗口输入：

>> y3=-exp(-x1)+x1.^2-2*x1+2

求得解析解：

y3 =

1.0000 0.9513 0.9052 0.8618 0.8213 0.7837 0.7492 0.7178 0.6897 0.6649 0.6435 0.6256 0.6112 0.6005 0.5934 0.5901 0.5907 0.5951 0.6034 0.6158 0.6321 0.6526 0.6771 0.7059

0.7388 0.7760 0.8175 0.8633 0.9134 0.9679 1.0269 1.0903

1.1581 1.2305 1.3073 1.3887 1.4747 1.5653 1.6604 1.7602 1.8647

输入：

>> plot(x1,y1,'o');hold on

>> plot(x1,y3,'*');hold on

可得近似值与解析解的图像比较：

从上面结果可以看出，执行40步时近似解的值要接近于解析解，误差更小，结果更精确。

5、)1/(1)(,1)0(,222't t y y ty y -=== 解:

执行20步时：

编写函数文件doty.m 如下： function f=doty(x,y); f=2*x*y^2;

在Matlab 命令窗口输入： >> [x1,y1]=DEEuler('doty',0,2,1,20) 可得到结果：

x1 =

0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 1.1000 1.2000 1.3000 1.4000 1.5000 1.6000 1.7000 1.8000 1.9000 2.0000 y1 =

1.0e+176 *

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 6.4573

在Matlab 命令窗口输入： >> y3=1./(1-x1.^2)

求得解析解：

y3 =

1.0e+015 *

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 4.5036 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000

输入：

>> plot(x1,y1,'o');hold on

>> plot(x1,y3,'*');hold on

可得近似值与解析解的图像比较：

执行40步时：

在Matlab 命令窗口输入：

>> [x1,y1]=DEEuler('doty',0,2,1,40)

可得到结果：

x1 =

0.0000 0.0500 0.1000 0.1500 0.2000 0.2500 0.3000 0.3500 0.4000 0.4500 0.5000 0.5500 0.6000 0.6500 0.7000 0.7500

0.8000 0.8500 0.9000 0.9500 1.0000 1.0500

1.1000 1.1500 1.2000 1.2500 1.3000 1.3500 1.4000 1.4500 1.5000 1.5500 1.6000 1.6500 1.7000 1.7500 1.8000 1.8500 1.9000 1.9500

2.0000

y1 =

1.0e+224 *

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 3.6893 Inf Inf Inf Inf Inf Inf Inf Inf Inf 在Matlab 命令窗口输入： >> y3=1./(1-x1.^2) 求得解析解： y3 =

1.0e+015 *

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -2.2518 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 输入：

>> plot(x1,y1,'o');hold on >> plot(x1,y3,'*');hold on

可得近似值与解析解的图像比较：

欧拉方法（后项）：

1、22)(,1)0(,2

2

'

+-+-==-=-t t e

t y y y t y t

解:

执行20步时：

编写函数文件doty.m 如下：

function f=doty(x,y);

f=x^2-y;

在Matlab 命令窗口输入：

>> [x1,y1]=BAEuler('doty',0,2,1,20)

可得到结果：

x1 =

0.0000 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000

0.8000 0.9000 1.0000 1.1000 1.2000 1.3000 1.4000 1.5000

1.6000 1.7000 1.8000 1.9000

2.0000

y1 =

1.0000 0.9110 0.8329 0.7665 0.7127 0.6719 0.6449 0.6323

0.6345 0.6520 0.6852 0.7345 0.8003 0.8829 0.9825 1.0995

1.2341 1.3864 1.5567 1.7452 1.9520

在Matlab 命令窗口输入：

>> y3=-exp(-x1)+x1.^2-2*x1+2

求得解析解：

y3 =

1.0000 0.9052 0.8213 0.7492 0.6897 0.6435 0.6112 0.5934

0.5907 0.6034 0.6321 0.6771 0.7388 0.8175 0.9134 1.0269

1.1581 1.3073 1.4747 1.6604 1.8647

输入：

>> plot(x1,y1,'o');hold on

>> plot(x1,y3,'*');hold on

可得近似值与解析解的图像比较：

执行40步时：

在Matlab 命令窗口输入：

>> [x1,y1]=BAEuler('doty',0,2,1,40)

可得到结果：

x1 =

0 0.0500 0.1000 0.1500 0.2000 0.2500 0.3000 0.3500 0.4000 0.4500 0.5000 0.5500 0.6000 0.6500 0.7000 0.7500

0.8000 0.8500 0.9000 0.9500 1.0000 1.0500 1.1000 1.1500

1.2000 1.2500 1.3000 1.3500 1.4000 1.4500 1.5000 1.5500

1.6000 1.6500 1.7000 1.7500 1.8000 1.8500 1.9000 1.9500

2.0000

y1 =

1.0000 0.9526 0.9079 0.8658 0.8267 0.7904 0.7572 0.7272 0.7003 0.6768 0.6566 0.6399 0.6268 0.6172 0.6114 0.6092 0.6109 0.6164 0.6258 0.6392 0.6566 0.6780 0.7035 0.7332

0.7671 0.8052 0.8475 0.8942 0.9451 1.0005 1.0602 1.1243

1.1929 1.2660 1.3435 1.4256 1.5122 1.6034 1.6992 1.7996 1.9046

在Matlab 命令窗口输入：

>> y3=-exp(-x1)+x1.^2-2*x1+2

求得解析解：

y3 =

1.0000 0.9513 0.9052 0.8618 0.8213 0.7837 0.7492 0.7178

0.6897 0.6649 0.6435 0.6256 0.6112 0.6005 0.5934 0.5901 0.5907 0.5951 0.6034 0.6158 0.6321 0.6526 0.6771 0.7059 0.7388 0.7760 0.8175 0.8633 0.9134 0.9679 1.0269 1.0903 1.1581 1.2305 1.3073 1.3887 1.4747 1.5653 1.6604 1.7602 1.864 输入：

>> plot(x1,y1,'o');hold on

>> plot(x1,y3,'*');hold on

可得近似值与解析解的图像比较：

从上面结果可以看出，执行40步时近似解的值要接近于解析解，误差更小，结果更精确。但是后项欧拉要比前项接近解析解。

5、)1/(1)(,1)0(,22

2

'

t t y y ty y -=== 解:

执行20步时：

编写函数文件doty.m 如下： function f=doty(x,y); f=2*x*y^2;

在Matlab 命令窗口输入： >> [x1,y1]=BAEuler('doty',0,2,1,20) 可得到结果：

x1 =

0.0000 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000

0.8000 0.9000 1.0000 1.1000 1.2000 1.3000 1.4000 1.5000

1.6000 1.7000 1.8000 1.9000

2.0000

y1 =

1.0e+119 *

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 4.0544 Inf Inf Inf Inf Inf Inf Inf

在Matlab 命令窗口输入：

>> y3=1./(1-x1.^2)

求得解析解：

y3 =

1.0e+015 *

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 4.5036 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000

输入：

>> plot(x1,y1,'o');hold on

>> plot(x1,y3,'*');hold on

可得近似值与解析解的图像比较：

执行40步时：

在Matlab 命令窗口输入：

>> [x1,y1]=BAEuler('doty',0,2,1,40)

可得到结果：

x1 =

0 0.0500 0.1000 0.1500 0.2000 0.2500 0.3000 0.3500

0.4000 0.4500 0.5000 0.5500 0.6000 0.6500 0.7000 0.7500

0.8000 0.8500 0.9000 0.9500 1.0000 1.0500 1.1000 1.1500

1.2000 1.2500 1.3000 1.3500 1.4000 1.4500 1.5000 1.5500

1.6000 1.6500 1.7000 1.7500 1.8000 1.8500 1.9000 1.9500

2.0000

y1 =

1.0e+195 *

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.0834

Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf

在Matlab 命令窗口输入：

>> y3=1./(1-x1.^2)

求得解析解：

y3 =

1.0e+015 *

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

0.0000 0.0000 0.0000 0.0000 -2.2518 -0.0000 -0.0000 -0.0000

-0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000

-0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000

-0.0000

输入：

>> plot(x1,y1,'o');hold on

>> plot(x1,y3,'*');hold on

可得近似值与解析解的图像比较：

梯形算法：

1、22)(,1)0(,22'+-+-==-=-t t e t y y y t y t 解：

在matlab 窗口中输入：

>> [I,step,h2] = CombineTraprl('-exp(-x1)+x1^2-2*x1+2',0,2) 可得结果： I = 1.8032 step = 29 h2 =

0.0690

由此可以看出，梯形算法求得的结果比执行20步的欧拉算法要接近解析解，但执行40步的近似解要比梯形算法精确。

5、)1/(1)(,1)0(,22

2

'

t t y y ty y -=== 解：

在matlab 窗口中输入：

>> [I,step,h2] = CombineTraprl('1/(1-x^2)',0,2)

无法得出结果，程序报 Warning: Divide by zero. 错误

改进欧拉方法：

1、22)(,1)0(,22'+-+-==-=-t t e t y y y t y t 解:

执行20步时：

编写函数文件doty.m 如下： function f=doty(x,y); f=x^2-y;

在Matlab 命令窗口输入： >> [x1,y1]=MoEuler('doty',0,2,1,20) 可得到结果：

x1 =

0.0000 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 1.1000 1.2000 1.3000 1.4000 1.5000 1.6000 1.7000 1.8000 1.9000 2.0000 y1 =

1.0000 0.9055 0.8219 0.7501 0.6909 0.6450 0.6130 0.5954 0.5929 0.6059 0.6348 0.6800 0.7418 0.8207 0.9167 1.0304 1.1617 1.3111 1.4786 1.6644 1.8687 在Matlab 命令窗口输入： >> y3=-exp(-x1)+x1.^2-2*x1+2 求得解析解：

y3 =

1.0000 0.9052 0.8213 0.7492 0.6897 0.6435 0.6112 0.5934 0.5907 0.6034 0.6321 0.6771 0.7388 0.8175 0.9134 1.0269 1.1581 1.3073 1.4747 1.6604 1.8647 输入：

>> plot(x1,y1,'o');hold on >> plot(x1,y3,'*');hold on

可得近似值与解析解的图像比较：

执行40步时：

在Matlab 命令窗口输入：

>> [x1,y1]=MoEuler('doty',0,2,1,40)

可得到结果：

x1 =

0.0000 0.0500 0.1000 0.1500 0.2000 0.2500 0.3000 0.3500 0.4000 0.4500 0.5000 0.5500 0.6000 0.6500 0.7000 0.7500

0.8000 0.8500 0.9000 0.9500 1.0000 1.0500 1.1000 1.1500

1.2000 1.2500 1.3000 1.3500 1.4000 1.4500 1.5000 1.5500

1.6000 1.6500 1.7000 1.7500 1.8000 1.8500 1.9000 1.9500

2.0000

y1 =

1.0000 0.9513 0.9052 0.8619 0.8214 0.7839 0.7494 0.7181 0.6900 0.6652 0.6438 0.6260 0.6116 0.6009 0.5939 0.5907 0.5912 0.5957 0.6040 0.6164 0.6328 0.6532 0.6778 0.7066

0.7395 0.7768 0.8182 0.8641 0.9142 0.9688 1.0277 1.0911

1.1590 1.2313 1.3082 1.3897

1.4756 1.5662 1.6614 1.7612 1.8657

在Matlab 命令窗口输入：

>> y3=-exp(-x1)+x1.^2-2*x1+2

求得解析解：

y3 =

1.0000 0.9513 0.9052 0.8618 0.8213 0.7837 0.7492 0.7178 0.6897 0.6649 0.6435 0.6256 0.6112 0.6005 0.5934 0.5901

0.5907 0.5951 0.6034 0.6158 0.6321 0.6526 0.6771 0.7059 0.7388 0.7760 0.8175 0.8633 0.9134 0.9679 1.0269 1.0903 1.1581 1.2305 1.3073 1.3887 1.4747 1.5653 1.6604 1.7602 1.8647 输入：

>> plot(x1,y1,'o');hold on

>> plot(x1,y3,'*');hold on

可得近似值与解析解的图像比较：

从上面结果可以看出，改进的欧拉法更接近于解析解，并且执行40步时近似解的值要接近于解析解，误差更小，结果更精确。

5、)1/(1)(,1)0(,22

2

'

t t y y ty y -=== 解:

执行20步时：

编写函数文件doty.m 如下： function f=doty(x,y); f=2*x*y^2;

在Matlab 命令窗口输入： >> [x1,y1]=MoEuler('doty',0,2,1,20) 可得到结果： x1 =

0.0000 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000

0.8000 0.9000 1.0000 1.1000 1.2000 1.3000 1.4000 1.5000

1.6000 1.7000 1.8000 1.9000

2.0000

y1 =

1.0e+113 *

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 3.6181 Inf Inf Inf Inf Inf Inf

在Matlab 命令窗口输入：

>> y3=1./(1-x1.^2)

求得解析解：

y3 =

1.0e+015 *

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 4.5036 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000

输入：

>> plot(x1,y1,'o');hold on

>> plot(x1,y3,'*');hold on

可得近似值与解析解的图像比较：

执行40步时：

在Matlab 命令窗口输入：

>> [x1,y1]=MoEuler('doty',0,2,1,40)

可得到结果：

x1 =

0.0000 0.0500 0.1000 0.1500 0.2000 0.2500 0.3000 0.3500

0.4000 0.4500 0.5000 0.5500 0.6000 0.6500 0.7000 0.7500

0.8000 0.8500 0.9000 0.9500 1.0000 1.0500 1.1000 1.1500

1.2000 1.2500 1.3000 1.3500 1.4000 1.4500 1.5000 1.5500

1.6000 1.6500 1.7000 1.7500 1.8000 1.8500 1.9000 1.9500

2.0000

y1 =

1.0e+102 *

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

2.1796 Inf Inf Inf Inf Inf Inf Inf

Inf Inf Inf Inf Inf Inf Inf Inf Inf 在Matlab 命令窗口输入：

>> y3=1./(1-x1.^2)

求得解析解：

y3 =

1.0e+015 *

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

0.0000 0.0000 0.0000 0.0000 -2.2518 -0.0000 -0.0000 -0.0000

-0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000

-0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000

-0.0000

输入：

>> plot(x1,y1,'o');hold on

>> plot(x1,y3,'*');hold on

可得近似值与解析解的图像比较：

四阶龙格-库塔法

1、22)(,1)0(,22'+-+-==-=-t t e t y y y t y t 解：

执行20步时：

在matlab 窗口中输入： >> syms x y; >> z=x^2-y

>> y= DELGKT4_lungkuta(z,0.1,0,2,1,[x y])

可以得出下面结果： y = 1.0000 0.9052 0.8213 0.7492 0.6897 0.6435 0.6112 0.5934 0.5907 0.6034

初值问题

《计算机数学基础(2)》辅导六 第14章常微分方程的数值解法 一、重点内容 1.欧拉公式： (k＝0，1，2，…，n－1) 局部截断误差是O(h2)。 2. 改进欧拉公式： 或表示成： 平均形式： 局部截断误差是O(h3)。 3. 四阶龙格――库塔法公式： 其中κ1＝f(x k，y k)；κ2＝f(x k+ 0.5h，y k+ 0.5 hκ1)；κ3＝f(x k+ 0.5 h，y k+ 0.5 hκ2)； κ4＝f(x k+h，y k+hκ3) 局部截断误差是O(h5)。

二、实例 例1用欧拉法解初值问题 取步长h＝0.2。计算过程保留4位小数。 解h＝0.2，f(x，y)＝－y－xy2。首先建立欧拉迭代格式 ＝0.2y k(4－x k y k) (k＝0，1，2) 当k＝0，x1＝0.2时，已知x0＝0，y0＝1，有 y(0.2)≈y1＝0.2×1(4－0×1)＝0.8 当k＝1，x2＝0.4时，已知x1＝0.2，y1＝0.8，有 y(0.4)≈y2＝0.2×0.8×(4－0.2×0.8)＝0.6144 当k＝2，x3＝0.6时，已知x2＝0.4，y2＝0.6144，有 y(0.6)≈y3＝0.2×0.6144×(4－0.4×0.6144)＝0.4613 例2 用欧拉预报－校正公式求解初值问题 取步长h＝0.2，计算y(1.2)，y(1.4)的近似值，小数点后至少保留5位。 解步长h＝0.2，此时f(x，y)＝－y－y2sin x 欧拉预报－校正公式为： 有迭代格式：

当k＝0，x0＝1，y0＝1时，x1＝1.2，有 ＝y0(0.8－0.2y0sin x0)＝1×(0.8－0.2×1sin1)＝0.63171 y(1.2)≈y1 ＝1×(0.9－0.1×1×sin1)－0.1(0.63171＋0.631712sin1.2)＝0.71549 当k＝1，x1＝1.2，y1＝0.71549时，x2＝1.4，有 ＝y1(0.8－0.2y1sin x1)＝0.71549×(0.8－0.2×0.71549sin1.2) ＝0.47697 y(1.4)≈y2 ＝0.71549×(0.9－0.1×0.71549×sin1.2) －0.1(0.47697＋0.476972sin1.4) ＝0.52611 例3写出用四阶龙格――库塔法求解初值问题 的计算公式，取步长h＝0.2计算y(0.4)的近似值。至少保留四位小数。 解此处f(x，y)＝8－3y，四阶龙格――库塔法公式为 其中κ1＝f(x k，y k)；κ2＝f(x k+ 0.5h，y k+ 0.5 hκ1)；κ3＝f(x k+ 0.5 h，y k+ 0.5 hκ2)；

龙格库塔方法matlab实现

龙格库塔方法matlab实现~ function ff=rk(yy,x0,y0,h,a,b)%yy为y的导函数,x0,y0,为初值,h为步长,a,b为区间 c=(b-a)/h+1;i1=1; %c为迭代步数；i1为迭代步数累加值 y=y0;z=zeros(c,6); %z生成c行，5列的零矩阵存放结果； %每行存放c次迭代结果，每列分别存放k1~k4及y的结果 for x=a:h:b if i1<=c k1=feval(yy,x,y); k2=feval(yy,x+h/2,y+(h*k1)/2); k3=feval(yy,x+h/2,y+(h*k2)/2); k4=feval(yy,x+h,y+h*k3); y=y+(h/6)*(k1+2*k2+2*k3+k4); z(i1,1)=x;z(i1,2)=k1;z(i1,3)=k2;z(i1,4)=k3;z(i1,5)=k4;z(i1,6)=y; i1=i1+1; end end fprintf(‘结果矩阵，第一列为x(n),第二列~第五列为k1~k4，第六列为y（n+1）的结果') z %在命令框输入下列语句 %yy=inline('x+y'); %>> rk(yy,0,1,0.2,0,1) %将得到结果 %结果矩阵，第一列为x(n),第二列~第五列为k1~k4第六列为y（n+1）的结果 %z = % 0 1.0000 1.2000 1.2200 1.4440 1.2428 % 0.2000 1.4428 1.6871 1.7115 1.9851 1.5836 % 0.4000 1.9836 2.2820 2.3118 2.6460 2.0442 % 0.6000 2.6442 3.0086 3.0451 3.4532 2.6510 % 0.8000 3.4510 3.8961 3.9407 4.4392 3.4365 % 1.0000 4.4365 4.9802 5.0345 5.6434 4.4401

四阶龙格库塔和向前欧拉方法和隐式欧拉法以及改进欧拉法教材

实验九欧拉方法 信息与计算科学金融崔振威201002034031 一、实验目的： 1、掌握欧拉算法设计及程序实现 二、实验内容： 1、p364-9.2.4、p386-9.5.6 三、实验要求： 主程序： 欧拉方法（前项）： function [x,y]=DEEuler(f,a,b,y0,n); %f:一阶常微分方程的一般表达式的右端函数 %a:自变量的取值下限 %b:自变量的取值上限 %y0：函数的初值 %n:积分的步数 if nargin<5,n=50; end h=(b-a)/n; x(1)=a;y(1)=y0; for i=1:n x(i+1)=x(i)+h; y(i+1)=y(i)+h*feval(f,x(i),y(i)); end format short 欧拉方法（后项）： function [x,y]=BAEuler(f,a,b,y0,n); %f:一阶常微分方程的一般表达式的右端函数 %a:自变量的取值下限 %b:自变量的取值上限 %y0：函数的初值 %n:积分的步数 if nargin<5,n=50; end h=(b-a)/n; x(1)=a;y(1)=y0; for i=1:n x(i+1)=x(i)+h; y1=y(i)+h*feval(f,x(i),y(i));

y(i+1)=y(i)+h*feval(f,x(i+1),y1); end format short 梯形算法： function [I,step,h2] = CombineTraprl(f,a,b,eps) %f 被积函数 %a,b 积分上下限 %eps 精度 %I 积分结果 %step 积分的子区间数 if(nargin ==3) eps=1.0e-4; end n=1; h=(b-a)/2; I1=0; I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h; while abs(I2-I1)>eps n=n+1; h=(b-a)/n; I1=I2; I2=0; for i=0:n-1 x=a+h*i; x1=x+h; I2=I2+(h/2)*(subs(sym(f),findsym(sym(f)),x)+subs(sym(f),findsym(sym(f)),x1)); end end I=I2; step=n; h2=(b-a)/n; 改进欧拉方法： function [x,y]=MoEuler(f,a,b,y0,n); %f:一阶常微分方程的一般表达式的右端函数 %a:自变量的取值下限 %b:自变量的取值上限 %y0：函数的初值 %n:积分的步数 if nargin<5,n=50; end

龙格库塔积分算法

龙格库塔法 龙格库塔法是常用于模拟常微分方程的解的重要的一类隐式或显式迭代法。这些技术由数学家C. Runge和M.W. Kutta于1900年左右发明。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。 龙格库塔法是一种在工程上应用广泛的高精度单步算法，可以应用在物理、工程、控制、动力学中，如模糊控制、弹道分析以及分析光纤特性等，在系统仿真中得到广泛应用。 龙格库塔法源自于相应的泰勒级数方法，在每一插值节点用泰勒级数展开,其截断误差阶数也是，根据可省略更高阶的导数计算, 这种方法可构造任意阶数的龙格库塔法。其中4 阶龙格库塔法是最常用的一种方法。因为它相当精确、稳定、容易编程。在计算中一般不必使用高阶方法, 因为附加的计算误差可由增加精度来弥补。如果需要较高的精度, 可采取减小步长的方法即可。4 阶龙格库塔法的精度类似4 阶泰勒级数法的精度。 1、初值问题 对于一阶常微分方程的初值问题 根据常微分方程的理论可知，此初值问题的解在区间[a，b]上存在，且唯一。 2、离散化

取步长h=(b-a)/n，将区间[a , b]分成n个子区间： a=<=b 在其中任意两点的曲线段上，根据积分中值定理，一段光滑曲 线上至少有一点，它的斜率与整段曲线的平均斜率相同， 得=y’() (0<<1) 其中，= 可以将上式改写成y()=y()+h*K (2.1) 其中K为平均斜率，K=f() 公式（2.1）表明，如果能够确定平均斜率K，就可以根据（2.1）式得到y()的值。 欧拉法和龙格库塔法就是用不同方法确定不同精度的平均斜率K，从而求得y()的近似值。 3、Euler法 欧拉法虽然精度低，但它是最简单的一种显式单步法，也是龙 格库塔法的基础。 首先，令、为y() 及y()的近似值，并且令平均斜 率K=f()，即以点的斜率作为平均斜率K，便得到欧拉公式=+h* f() （3.1） 4、改进的欧拉法 此种方法是取、两点的斜率的平均值作为平均斜率K， 即K= ,其中、均为y()以及y()的近似值，就得到 改进后的欧拉公式（4.1）

Matlab中龙格-库塔(Runge-Kutta)方法原理及实现

函数功能编辑本段回目录 ode是专门用于解微分方程的功能函数，他有ode23,ode45,ode23s等等，采用的是Runge-Kutta算法。ode45表示采用四阶，五阶runge-kutta单步算法,截断误差为(Δx)3。解决的是Nonstiff(非刚性)的常微分方程.是解决数值解问题的首选方法,若长时间没结果,应该就是刚性的,换用ode23来解. 使用方法编辑本段回目录 [T,Y] = ode45(odefun,tspan,y0) odefun 是函数句柄,可以是函数文件名,匿名函数句柄或内联函数名 tspan 是区间[t0 tf] 或者一系列散点[t0,t1,...,tf] y0 是初始值向量 T 返回列向量的时间点 Y 返回对应T的求解列向量 [T,Y] = ode45(odefun,tspan,y0,options) options 是求解参数设置,可以用odeset在计算前设定误差,输出参数,事件等 [T,Y,TE,YE,IE] =ode45(odefun,tspan,y0,options) 在设置了事件参数后的对应输出 TE 事件发生时间 YE 事件解决时间 IE 事件消失时间 sol =ode45(odefun,[t0 tf],y0...) sol 结构体输出结果 应用举例编辑本段回目录 1 求解一阶常微分方程

程序: 一阶常微分方程 odefun=@(t,y) (y+3*t)/t^2; %定义函数 tspan=[1 4]; %求解区间 y0=-2; %初值 [t,y]=ode45(odefun,tspan,y0); plot(t,y) %作图 title('t^2y''=y+3t,y(1)=-2,1

数值计算方法复习题9

习题九 1. 取步长h = 0.1，分别用欧拉法与改进的欧拉法解下列初值问题 （1）；（2） 准确解：（1）；（2）； 欧拉法：，，， 改进的欧拉法：，，， 2. 用四阶标准龙格—库塔法解第1题中的初值问题，比较各法解的精度。，，， 3. 用欧拉法计算下列积分在点处的近似值。 0.5000,1.1420,2.5011,7.2450 4. 求下列差分格式局部截断误差的首项，并指出其阶数。 （1），2 （2），3； （3），4 （4），4 5.用Euler法解初值问题取步长h=0.1，计算到x=0.3(保留到小数点后4位).

解: 直接将Eulerr法应用于本题，得到 由于，直接代入计算，得到 6.用改进Euler法和梯形法解初值问题取步长 h=0.1,计算到x=0.5，并与准确解相比较. 解：用改进Euler法求解公式，得 计算结果见下表 用梯形法求解公式，得 解得 精确解为 7.证明中点公式(7.3.9)是二阶的，并求其局部截断误差主项. 证明根据局部截断误差定义，得 将右端Taylor展开，得

故方法是二阶的，且局部截断误差主项是上式右端含h3的项。 8.用四阶R-K方法求解初值问题取步长 h=0.2. 解直接用四阶R－K方法 其中 计算结果如表所示： 9.对于初值问题 解因f'(y)=-100，故由绝对稳定区间要求(1)用Euler法解时， (2)用梯形法解时，绝对稳定区间为，由因f 对y是线性的，故不用迭代，对h仍无限制。(3)用四阶R-K方法时， 10. (1) 用Euler法求解，步长h应取在什么范围内计算才稳定?(2) 若用梯形法求解，对步长h有无限制? (3) 若用四阶R-K方法求解，步长h如何选取?

龙格-库塔法

龙格－库塔法（Runge-Kutta） 数值分析中，龙格－库塔法（Runge-Kutta）是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。这些技术由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。 经典四阶龙格库塔法 龙格库塔法的家族中的一个成员如此常用，以至于经常被称为“RK4”或者就是“龙格库塔法”。令初值问题表述如下。 则，对于该问题的RK4由如下方程给出： 其中 这样，下一个值(yn+1)由现在的值(yn)加上时间间隔(h)和一个估算的斜率的乘积决定。该斜率是以下斜率的加权平均： k1是时间段开始时的斜率； k2是时间段中点的斜率，通过欧拉法采用斜率k1来决定y在点tn + h/2的值； k3也是中点的斜率，但是这次采用斜率k2决定y值； k4是时间段终点的斜率，其y值用k3决定。 当四个斜率取平均时，中点的斜率有更大的权值： RK4法是四阶方法，也就是说每步的误差是h5阶，而总积累误差为h4阶。 注意上述公式对于标量或者向量函数(y可以是向量)都适用。 显式龙格库塔法 显示龙格－库塔法是上述RK4法的一个推广。它由下式给出 其中

(注意：上述方程在不同著述中由不同但却等价的定义)。 要给定一个特定的方法，必须提供整数s (阶段数)，以及系数aij (对于1 ≤ j < i ≤ s), bi (对于i = 1, 2, ..., s)和ci (对于i = 2, 3, ..., s)。这些数据通常排列在一个助记工具中，称为龙格库塔表： c2 a21 c3 a31 a32 cs as1 as2 as,s ? 1 b1 b2 bs ? 1 bs 龙格库塔法是自洽的，如果 如果要求方法有精度p则还有相应的条件，也就是要求舍入误差为O(hp+1)时的条件。这些可以从舍入误差本身的定义中导出。例如，一个2阶精度的2段方法要求b1 + b2 = 1, b2c2 = 1/2, 以及b2a21 = 1/2。

计算方法及答案汇总

《计算方法》练习题一 一、填空题 1． 14159.3=π的近似值3.1428，准确数位是（ ）。 2．满足d b f c a f ==)(,)(的插值余项=)(x R （ ）。 3．设)}({x P k 为勒让德多项式，则=))(),((22x P x P （ ）。 4．乘幂法是求实方阵（ ）特征值与特征向量的迭代法。 5．欧拉法的绝对稳定实区间是（ ）。 6． 71828.2=e 具有3位有效数字的近似值是（ ）。 7．用辛卜生公式计算积分 ?≈+1 01x dx （ ） 。 8．设)()1()1(--=k ij k a A 第k 列主元为) 1(-k pk a ，则=-)1(k pk a （ ）。 9．已知?? ? ? ??=2415A ，则=1A （ ）。 10．已知迭代法：),1,0(),(1 ==+n x x n n ? 收敛，则)(x ?'满足条件（ ）。 二、单选题 1．已知近似数,,b a 的误差限)(),(b a εε，则=)(ab ε（ ）。 A ．)()(b a εε Ｂ．)()(b a εε+ Ｃ．)()(b b a a εε+ Ｄ．)()(a b b a εε+ 2．设x x x f +=2 )(，则=]3,2,1[f （ ）。 Ａ．１ Ｂ．２ Ｃ．３ Ｄ．４ 3．设Ａ＝?? ? ? ??3113，则化Ａ为对角阵的平面旋转=θ（ ）． Ａ． 2π Ｂ．3π Ｃ．4π Ｄ．6 π 4．若双点弦法收敛，则双点弦法具有（ ）敛速． Ａ．线性 Ｂ．超线性 Ｃ．平方 Ｄ．三次 5．改进欧拉法的局部截断误差阶是（ ）. A ．)(h o Ｂ．)(2h o Ｃ．)(3h o Ｄ．)(4 h o 6．近似数2 1047820.0?=a 的误差限是（ ）。 Ａ． 51021-? Ｂ．41021-? Ｃ．31021-? Ｄ．2102 1 -? 7．矩阵Ａ满足（ ）,则存在三角分解A=LR 。 A ．0det ≠A Ｂ． )1(0det n k A k <≤≠ Ｃ．0det >A Ｄ．0det

龙格-库塔法

龙格－库塔法 格－库塔法（Runge-Kutta）数值分析中，龙格－库塔法（Runge-Kutta）是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。这些技术由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。 经典四阶龙格库塔法 龙格库塔法的家族中的一个成员如此常用，以至于经常被称为“RK4”或者就是“龙格库塔法”。 令初值问题表述如下。 则，对于该问题的RK4由如下方程给出： 其中 这样，下一个值(y n+1)由现在的值(y n)加上时间间隔(h)和一个估算的斜率的乘积决定。该斜率是以下斜率的加权平均： ?k1是时间段开始时的斜率； ?k2是时间段中点的斜率，通过欧拉法采用斜率k1来决定y在点t n + h/2的值； ?k3也是中点的斜率，但是这次采用斜率k2决定y值； ?k4是时间段终点的斜率，其y值用k3决定。

当四个斜率取平均时，中点的斜率有更大的权值： RK4法是四阶方法，也就是说每步的误差是h5阶，而总积累误差为h4阶。 注意上述公式对于标量或者向量函数(y可以是向量)都适用。 显式龙格库塔法 显示龙格－库塔法是上述RK4法的一个推广。它由下式给出 其中 (注意：上述方程在不同著述中由不同但却等价的定义)。 要给定一个特定的方法，必须提供整数s (阶段数)，以及系数 a ij (对于1 ≤ j < i ≤ s), b i (对于i = 1, 2, ..., s)和c i (对于i = 2, 3, ..., s)。这些数据通常排列在一个助记工具中，称为龙格库塔表： 0 c2 a21 c3 a31 a32c s a s1 a s2a s,s ? 1 b1 b2b s ? 1 b s 龙格库塔法是自洽的，如果

欧拉法与龙格库塔法比较分析

解微分方程的欧拉法，龙格-库塔法简单实例比较 欧拉方法（Euler method)用以对给定初值的常微分方程(即初值问题)求解分为前EULER 法、后退EULER 法、改进的EULER 法。 缺点： 欧拉法简单地取切线的端点作为下一步的起点进行计算，当步数增多时，误差会因积累而越来越大。因此欧拉格式一般不用于实际计算。 改进欧拉格式（向前欧拉公式）： 为提高精度，需要在欧拉格式的基础上进行改进。采用区间两端的斜率的平均值作为直线方程的斜率。改进欧拉法的精度为二阶。 算法： 微分方程的本质特征是方程中含有导数项，数值解法的第一步就是设法消除其导数值。对于常微分方程： (,)dy f x y dx = [,]x a b ∈ 0()y a y = 可以将区间[,]a b 分成n 段，那么方程在第i x 点有'()(,())i i i y x f x y x =，再用向前差商近似代替导数则为： ((1)())(,())i i i i y x y x f x y x h +-= 在这里，h 是步长，即相邻两个结点间的距离。因此可以根据i x 点和i y 的数值计算出1i y +来： 1(,)i i i i y y h f x y +=+?0,1,2,i L =L 这就是向前欧拉公式。 改进的欧拉公式：

将向前欧拉公式中的导数(,)i i f x y 改为微元两端导数的平均，即上式便是梯形的欧拉公式。 可见，上式是隐式格式，需要迭代求解。为了便于求解，使用改进的欧拉公式： 数值分析中，龙格－库塔法（Runge-Kutta ）是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。实际上，龙格-库塔法是欧拉方法的一种推广，向前欧拉公式将导数项简单取为(,)n n f x y ,而改进的欧拉公式将导数项取为两端导数的平均。 龙格-库塔方法的基本思想： 在区间1[,]n n x x +内多取几个点，将他们的斜率加权平均，作为导数的近似。 龙格库塔法的家族中的一个成员如此常用，以至于经常被称为“RK4”或者就是“龙格库塔法”。 令初值问题表述如下。 '(,)y f t y =00()y t y = 则，对于该问题的RK4由如下方程给出： 11234(22)6 n n h y y k k k k +=++++ 其中 1(,)n n k f t y = 21(,)22 n n h h k f t y k =++ 32(,)22 n n h h k f t y k =++ 43(,)n n k f t h y hk =++

Matlab中龙格-库塔(Runge-Kutta)方法原理及实现

Matlab中龙格-库塔(Runge-Kutta)方法原理及实现 龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。龙格库塔方法的理论基础来源于泰勒公式和使用斜率近似表达微分，它在积分区间多预计算出几个点的斜率，然后进行加权平均，用做下一点的依据，从而构造出了精度更高的数值积分计算方法。如果预先求两个点的斜率就是二阶龙格库塔法，如果预先取四个点就是四阶龙格库塔法。一阶常微分方程可以写作：y'=f(x,y),使用差分概念。 (Yn+1-Yn)/h= f(Xn,Yn)推出（近似等于，极限为Yn'） Yn+1=Yn+h*f(Xn,Yn) 另外根据微分中值定理，存在0

所以，为了更好更准确地把握时间关系，应自己在理解龙格库塔原理的基础上，编写定步长的龙格库塔函数，经过学习其原理，已经完成了一维的龙格库塔函数。 仔细思考之后，发现其实如果是需要解多个微分方程组，可以想象成多个微分方程并行进行求解，时间，步长都是共同的，首先把预定的初始值给每个微分方程的第一步，然后每走一步，对多个微分方程共同求解。想通之后发现，整个过程其实很直观，只是不停的逼近计算罢了。编写的定步长的龙格库塔计算函数： function [x,y]=runge_kutta1(ufunc,y0,h,a,b)%参数表顺序依次是微分方程组的函数名称，初始值向量，步长，时间起点，时间终点（参数形式参考了ode45函数） n=floor((b-a)/h);%求步数 x(1)=a;%时间起点 y(:,1)=y0;%赋初值，可以是向量，但是要注意维数 for ii=1:n x(ii+1)=x(ii)+h; k1=ufunc(x(ii),y(:,ii)); k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2); k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2); k4=ufunc(x(ii)+h,y(:,ii)+h*k3); y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6; %按照龙格库塔方法进行数值求解

证明隐式Euler方法稳定性

第六章数值积分 6.1数值积分基本概念 6.1.1引言 在区间上求定积分 (6.1.1) 是一个具有广泛应用的古典问题，从理论上讲，计算定积分可用Newton-Leibniz公式 (6.1.2) 其中F(x)是被积函数f(x)的原函数.但实际上有很多被积函数找不到用解析式子表达的原函数，例如等等，表面看它们并不复杂， 但却无法求得F(x).此外，有的积分即使能找到F(x)表达式，但式子非常复杂，计算也很困难.还有的被积函数是列表函数，也无法用(6.1.2)的公式计算.而数值积分则只需计算f(x) 在节点xi(i=0，1,…,n)上的值，计算方便且适合于在计算机上机械地实现. 本章将介绍常用的数值积分公式及其误差估计、求积公式的代数精确度、收敛性和稳定性以及Romberg求积法与外推原理等. 6.1.2插值求积公式 根据定积分定义，对及都有 (极限存在)若不取极限，则积分I(f)可近似表示为 (6.1.3) 这里称为求积节点，与f无关，称为求积系数，(6.1.3)称为机械求积公式. 为了得到形如(6.1.3)的求积公式，可在上用Lagrange插值多项式 ，则得

其中 (6.1.4) 这里求积系数由插值基函数积分得到，它与f(x)无关.如果求积公式(6.1.3)中的系数由(6.1.4)给出，则称(6.1.3)为插值求积公式.此时可由插值余项得到 (6.1.5) 这里ξ∈，(6.1.5)称为插值求积公式余项. 当n=1时，，此时 由(6.1.4)可得 于是 (6.1.6) 称为梯形公式.从几何上看它是梯形A bB(见图6-1)的面积近似曲线y=f(x)下的曲边梯形面积，公式(6.1.6)的余项为 (6.1.7)

第8章 常微分方程数值解法 本章主要内容： 1.欧拉法、改进欧拉法 2

第8章 常微分方程数值解法 本章主要内容： 1．欧拉法、改进欧拉法. 2．龙格-库塔法。 3．单步法的收敛性与稳定性。 重点、难点 一、微分方程的数值解法 在工程技术或自然科学中，我们会遇到的许多微分方程的问题，而我们只能对其中具有较简单形式的微分方程才能够求出它们的精确解。对于大量的微分方程问题我们需要考虑求它们的满足一定精度要求的近似解的方法，称为微分方程的数值解法。本章我们主要 讨论常微分方程初值问题?????==00 )() ,(y x y y x f dx dy 的数值解法。 数值解法的基本思想是：在常微分方程初值问题解的存在区间[a,b]内，取n+1个节点a=x 0＜x 1＜…＜x N =b （其中差h n = x n –x n-1称为步长，一般取h 为常数，即等步长），在这些节点上把常微分方程的初值问题离散化为差分方程的相应问题，再求出这些点的上的差分方程值作为相应的微分方程的近似值（满足精度要求）。 二、欧拉法与改进欧拉法 欧拉法与改进欧拉法是用数值积分方法对微分方程进行离散化的一种方法。 将常微分方程),(y x f y ='变为() *+=?++1 1))(,()()(n x n x n n dt t y t f x y x y 1．欧拉法（欧拉折线法） 欧拉法是求解常微分方程初值问题的一种最简单的数值解法。 欧拉法的基本思想：用左矩阵公式计算（＊）式右端积分，则得欧拉法的计算公式为：N a b h N n y x hf y y n n n n -= -=+=+)1,...,1,0(),(1 欧拉法局部截断误差 11121 )(2 ++++≤≤''=n n n n n x x y h R ξξ或简记为O （h 2）。

5.2龙格库塔法

第五章 常微分方程的数值解法 5.2 龙格-库塔法 一、教学目标及基本要求 通过对本节课的学习，使学生掌握常微分方程、常微分方程方程组的数值解法。 二、教学内容及学时分配 本节课主要介绍常微分方程的数值解法。具体内容如下： 讲授内容：龙格库塔方法、收敛性与稳定性 三、教学重点难点 1．教学重点：龙格库塔方法、收敛性与稳定性。 2. 教学难点：收敛性与稳定性。 四、教学中应注意的问题 多媒体课堂教学为主。适当提问，加深学生对概念的理解。 五、正文 龙格－库塔方法 引言 龙格－库塔方法的基本思想 '11()() ()()()(,())n n n n y x y x y y x y x hf y h ξξξ++-=?=+ 令*(,())K f y ξξ=称为区间1[,]n n x x +上的平均斜率，只要对*K 提供一种算法，即可推导出一种计算公式。 欧拉公式只是取点n x 的斜率作为区间1[,]n n x x +的平均斜率*K ，精度自然很低。 考察改进的欧拉公式 11211 ()22 n n y y h k k +=++1(,)n n k f x y =21(,)n n k f x h y hk =++ 它利用n x ，1n x +两点的斜率取算术平均，1n x +处斜率通过已知信息n y 用欧拉公式预报得到。

可以考虑设法在1[,]n n x x +多取几个点的斜率值，将它们加权平均作为区间 1[,]n n x x +的平均斜率*K 。这就是龙格－库塔方法的基本思想。 1、二阶龙格－库塔方法 考察1[,]n n x x +内一点,01n p n x x hp p +=+<≤，希望通过,n n p x x +两个点的斜率 12,K K 加权平均得到*K ，即令112((1))n n y y h K K λλ+=+-+ 取1(,)n n K f x y =，如何预报n p x +处斜率2K ？ 仿照改进的欧拉公式，先用欧拉公式预测()n p y x +的值n p y +： 1n p n y y phK +=+ 然后用n p y +计算2(,)n p n p K f x y ++=，从而得 112121[(1)](,)(,) n n n n n p n y y h K K K f x y K f x y phK λλ++=+-+==+ 适当选取,p λ，使上述公式具有较高得精度。假定()n n y y x =，分别将12 ,K K 泰勒展开： '1(,)()n n n K f x y y x == 212 ' '' 2 (,) (,)[(,)(,)(,)]()()()() n p n n n x n n n n y n n n n K f x y pkK f x y ph f x y f x y f x y O h y x phy x O h +=+=+++=++ 代入得 '2''31()()()()n n n n y y x hy x ph y x O h λ+=+++ 按泰勒展开2' '' 31()()()()2 n n n n h y y x hy x y x O h +=+++ 比较得，只要1 2 p λ= ，公式截断误差为3()O h 特别，当1,1/2p λ==，就是改进的欧拉公式， 改取1,1/2p λ==， 12n n y y hk +=+，1(,)n n k f x y =，21(,)22 n n h h k f x y k =++

数值分析--6微分方程数值解习题课

微分方程 初值问题数值解 习题课

一、应用向前欧拉法和改进欧拉法求由如下积分 2 x t y e dt -=? 所确定的函数y 在点x =0.5,1.0,1.5的近似值。 解：该积分问题等价于常微分方程初值问题 2 '(0)0x y e y -?=??=?? 其中h=0.5。其向前欧拉格式为 2 ()100ih i i y y he y -+?=+?? =?? 改进欧拉格式为 22()2(1)10()20 ih i h i i h y y e e y --++? =++???=? 将两种计算格式所得结果列于下表

二、应用4阶4步阿达姆斯显格式求解初值问题 '1(0)1y x y y =-+??=? 00.6x ≤≤ 取步长h=0.1. 解：4步显式法必须有4个起步值，0y 已知，其他3个123,,y y y 用4阶龙格库塔方法求出。 本题的信息有： 步长h=0.1；结点0.1(0,1, ,6)i x ih i i ===； 0(,)1,(0)1f x y x y y y =-+== 经典的4阶龙格库塔公式为

11234(22)6 i i h y y k k k k +=++++ 1(,)1i i i i k f x y x y ==-+ 121(,)0.05 1.0522 i i i i hk h k f x y x y k =++=--+ 232(,)0.05 1.0522 i i i i hk h k f x y x y k =++=--+ 433(,)0.1 1.1i i i i k f x h y hk x y k =++=--+ 算得1 1.0048375y =,2 1.0187309y =,3 1.0408184y = 4阶4步阿达姆斯显格式 1123(5559379) 24i i i i i i h y y f f f f +---=+-+- 1231 (18.5 5.9 3.70.90.24 3.24)24 i i i i i y y y y y i ---=+-+++ 由此算出 4561.0703231, 1.1065356, 1.1488186y y y === 三、用Euler 方法求 ()'1,0101 x y e y x x y =-++≤≤= 问步长h 应该如何选取，才能保证算法的稳定性？ 解：本题 (),1x f x y e y x =-++ (),0,01x y f x y e x λ'==-<≤≤

龙格库塔方法的Miline-Hamming预测-校正算法实验报告

2011-2012学年第2学期实验报告 实验名称：微分方程数值解实验学院：****** 专业：**************班级：********** 班内序号：** 学号：******** 姓名：****** 任课教师：****** 北京邮电大学 时间：****年**月**日

实验目标 用多环节Miline-Hamming预测-校正算法求下列方程的解 y‘=y?2x , y0=1, 0≤x≤4 其中解析解为 y x=1+2x 实验原理 计算龙格库塔显示公式计算预测值，然后用隐式公式进行校正。Miline-Hamming预测-校正公式为 p n+1=u n?3+4 h(2f n?f n?1+f n?2) m n+1=p n+1+112 121 (c n?p n) c n+1=1 9u n?u n?2+ 3 h[f t n+1,m n+1+2f n?f n?1] u n+1=c n+1? 9 (c n+1?p n+1) 其对应的算法流程为 1）输入a，b，f(t,u)，N，u0 2）置h=(b-a)/N，t0=a，n=1 3）计算f n-1=f(t n-1,u n-1) K1=hf n-1 K2=hf(t n-1+h/2, u n-1+K1/2) K3=hf(t n-1+h/2, u n-1+K2/2) K4=hf(t n-1+h, u n-1+K3) u n= u n-1+1/6(K1 +2K2 +2K3 +K4) t n=a+nh 4）输出(t n ,u n) 5）若n<3，置n+1→n，返回3； 否则，置n+1→n，0→p0，0→c0，转6. 6）计算 f3=f(t3,u3) t= t3+h p=u0+4/3(2 f3–f2 +2f1) m=p+112/121(c0-p0) c=1/8(9u3- u1)+3/8h[f(t,m)+ 2 f3–f2] u=c-9/121(c-p)

数值分析习 题 六 解 答

习 题 六 解 答 1、在区间［0，1］上用欧拉法求解下列的初值问题，取步长h=0.1。 (1)210(1)(0)2y y y '?=--?=?(2)sin (0)0 x y x e y -'?=+?=? 解：（1）取h=0.1，本初值问题的欧拉公式具体形式为 21(1)(0,1,2,)n n n y y y n +=--= 由初值y 0=y(0)=2出发计算，所得数值结果如下： x 0=0，y 0=2； x 1=0.1，2100(1)211y y y =--=-= x 2=0.2，2211(1)101y y y =--=-= 指出： 可以看出，实际上求出的所有数值解都是1。 （2）取h=0.1，本初值问题的欧拉公式具体形式为 21(sin )(0,1,2,)n x n n n y y h x e n -+=++= 由初值y 0=y(0)=0出发计算，所得数值结果如下： x 0=0，y 0=0； x 1=0.1， 02 1000 (sin )00.1(sin 0)00.1(01)0.1x y y h x e e -=++=+?+=+?+= x 2=0.2， 12 2110.1 (sin )0.10.1(sin 0.1)0.10.1(0.10.9)0.2 x y y h x e e --=++=+?+=+?+= 指出： 本小题的求解过程中，函数值计算需要用到计算器。 2、用欧拉法和改进的欧拉法(预测－校正法)求解初值问题，取步长h=0.1。 22(00.5) (0)1 y x y x y '?=-≤≤? =? 解：(1) 取h=0.1，本初值问题的欧拉公式具体形式为 2 1(2)(0,1,2,)n n n n y y h x y n +=+-= 由初值y 0=y(0)=1出发计算，所得数值结果如下： x 0=0，y 0=1； x 1=0.1，2 21000(2)10.1(021)0.8y y h x y =+-=+?-?= x 2=0.2，222111(2)0.80.1(0.120.8)0.641y y h x y =+-=+?-?=

龙格库塔方法及其matlab实现

龙格-库塔方法及其matlab实现 摘要：本文的目的数值求解微分方程精确解，通过龙格-库塔法，加以利用matlab为工具 达到求解目的。龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法，用于数值求解微分方程。MatLab软件是由美国Mathworks公司推出的用于数值计算和图形 处理的科学计算系统环境。MatLab是英文MATrix LABoratory（矩阵实验室）的缩写。在MratLab环境下，用户可以集成地进行程序设计、数值计算、图形绘制、输入输出、文件 管理等各项操作。 关键词：龙格-库塔 matlab 微分方程 1.前言 1.1：知识背景 龙格－库塔法（Runge-Kutta）是用于非线性常微分方程的解的重要的一类隐式或显式迭代法。这些技术由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。通常所说的龙格库塔方法是相对四阶龙格库塔而言的，成为经典四阶龙格库塔法。该方法具有精度高，收敛，稳定，计算过程中可以改变步长不需要计算高阶导数等优点，但是仍需计算在 一些点上的值，比如四阶龙格-库塔法没计算一步需要计算四步，在实际运用中是有一定复杂性的。 Matlab是在20世纪七十年代后期的事：时任美国新墨西哥大学计算机科学系主任的Cleve Moler教授出于减轻学生编程负担的动机，为学生设计了一组调用LINPACK和EISPACK库 程序的“通俗易用”的接口，此即用FORTRAN编写的萌芽状态的MATLAB。 经几年的校际流传，在Little的推动下，由Little、Moler、Steve Bangert合作，于1984年成立了MathWorks公司，并把MATLAB正式推向市场。从这时起，MATLAB的内核 采用C语言编写，而且除原有的数值计算能力外，还新增了数据图视功能。 MATLAB以商品形式出现后，仅短短几年，就以其良好的开放性和运行的可靠性， 使原先控制领域里的封闭式软件包（如英国的UMIST，瑞典的LUND和SIMNON，德国的KEDDC）纷纷淘汰，而改以MATLAB为平台加以重建。在时间进入20世纪九十年代的时候，MATLAB已经成为国际控制界公认的标准计算软件。 到九十年代初期，在国际上30几个数学类科技应用软件中，MATLAB在数值计算方面独占 鳌头，而Mathematica和Maple则分居符号计算软件的前两名。Mathcad因其提供计算、 图形、文字处理的统一环境而深受中学生欢迎。 1.2研究的意义 精确求解数值微分方程，对龙格库塔的深入了解与正确运用，主要是在已知方程导数和初 值信息，利用计算机仿真时应用，省去求解微分方程的复杂过程。利用matlab强大的数值计算功能，省去认为计算的过程，达到快速精确求解数值微分方程。在实际生活中可以利 用龙格库塔方法和matlab的完美配合解决问题。 1.3研究的方法 对实例的研究对比，实现精度的要求，龙格库塔是并不是一个固定的公式，所以只是对典 型进行分析

数值计算方法复习题9

习题九 1. 取步长h = ，分别用欧拉法与改进的欧拉法解下列初值问题（1）；（2） 准确解：（1）；（2）； 欧拉法：，，， 改进的欧拉法：，，， 2. 用四阶标准龙格—库塔法解第1题中的初值问题，比较各法解的精度。，，， 3. 用欧拉法计算下列积分在点处的近似值。 ,,, 4. 求下列差分格式局部截断误差的首项，并指出其阶数。（1），2 （2），3；（3），4 （4），4 5.用Euler法解初值问题取步长h=，计算到x=(保留到小数点后4位).

解:直接将Eulerr法应用于本题，得到 由于，直接代入计算，得到 6.用改进Euler法和梯形法解初值问题取步长h=,计算到x=，并与准确解相比较. 解：用改进Euler法求解公式，得 计算结果见下表 用梯形法求解公式，得 解得 精确解为 7.证明中点公式(7.3.9)是二阶的，并求其局部截断误差主项. 证明根据局部截断误差定义，得 将右端Taylor展开，得

故方法是二阶的，且局部截断误差主项是上式右端含h3的项。 8.用四阶R-K方法求解初值问题取步长h=.解直接用四阶R－K方法 其中 计算结果如表所示： 9.对于初值问题 解因f'(y)=-100，故由绝对稳定区间要求(1)用Euler法解时， (2)用梯形法解时，绝对稳定区间为，由因f对y是线性的，故不用迭代，对h仍无限制。(3)用四阶R-K方法时， 10. (1) 用Euler法求解，步长h应取在什么范围内计算才稳定(2) 若用梯形法求解，对步长h有无限制(3) 若用四阶R-K方法求解，步长h如何选取 解：用四阶显式Adams公式先要算出，而，其余3点可用四阶R-K方法计算。由，得