数值微分的计算方法
数值微分的计算方法
内容摘要 求解数值微分问题,就是通过测量函数在一些离散点上的值,求得函数的近似导数。本文就所学知识，归纳性地介绍了几种常用的数值微分计算方法。并举例说明计算，实验结果表明了方法的有效性。
关键词 数值微分 Taylor 展开式 Lagrange 插值 三对角矩阵
引言：数值微分即根据函数在一些离散点的函数值，推算它在某点的导数或高阶导数的近似值的方法。常见的可以用一个能够近似代替该函数的较简单的可微函数（如多项式或样条函数等）的相应导数作为能求导数的近似值，由此也可导出多点数值微分计算公式。当函数可微性不太好时，利用样条插值进行数值微分要比多项式插值更适宜。
1.Taylor 展开式方法
理论基础：Taylor 展开式
()()()()
()()
()()()00000002
2!
!
n
n
x x x x f x f x x x f x f x f x n --'''=+-+
++
+
我们借助Taylor 展开式，可以构造函数()f x 在点0x x =的一阶导数和二阶导数的数值微分公式。取步长0h >则
),()
(2
)()()(0011'
'20'
00h x x f h x hf x f h x f +∈++=+ξξ （1）
所以
),()
(2
)()()(0011'
'000'h x x f h h x f h x f x f +∈--+=
ξξ （2）
同理
),()
(2
)()()(0022'
'20'
00x h x f h x hf x f h x f -∈+-=-ξξ （3） ),()
(2
)()()(0022'
'000'x h x f h h h x f x f x f -∈+--=
ξξ （4）
式(2)和式(4)是计算()'
0f x 的数值微分公式，其截断误差为()O h ，为提高精度，将
Taylor 展开式多写几项
),()
(24
)(6)(2)()()(0011)
4(40'''30''20'
00h x x f h x f h x f h x hf x f h x f +∈++++=+ξξ ),()
(24
)(6)(2)()()(0022)
4(40'''30''20'
00x h x f h x f h x f h x hf x f h x f -∈+-+-=-ξξ
两式相减得
)()(6
2)()()(40'
''2000'
h O x f h h h x f h x f x f +---+= （5）
上式为计算)(0'x f 的微分公式，其截断误差为O(h 2
),比式（2）和(4)精度高。
两式相加，如果],[)(00)
4(h x h x C x f
+-∈，则有
2''
(4)
00002()2()()()()12f x h f x f x h h f x f h ξ+-+-=- ,00(,)
x h x h ξ∈-+ (6)
式(6)是计算)(0''x f 的数值微分公式,其截断误差为2()O h 。
例1．设函数()0ln ,2,0.1f x x x h ===，试用数值微分公式计算()'2f 的值。 解 由式(2)、式(4)和式(5)分别计算结果为
4879.01.02
ln 1.2ln )2('≈-≈f
5129.01
.09
.1ln 2ln )2('≈-≈f
5004.02
.09
.1ln 1.2ln )2('≈-≈
f
与真值()'
20.5f
=相比，式(5)计算的结果精度较高。
2.数值微分的Lagrange 插值方法
设函数()y f x =具有1n +个实验数据：()()(),0,1,2,,i i x f x i n =，我们希望估计
()'f x 的值，特别i x x =时，估计()'i f x 的值。
基于插值方法的数值微分做法是，由已知()()(),0,1,,i i x f x i n =建立Lagrange 插值多项式或Newton 插值多项式（这里以Lagrange 插值方法为例），即
()()()n n f x L x R x =+
于是
()()()'''n n f x L x R x =+
当i x x =时，有
'''
()()()(0,1,2,...,)n i n i i f x L x R x i n =+=
其中
1(1)'
'
()()()(1)!
n n n i i f R x x n ξω++=+
略去误差项有
''()()i n i f x L x ≈
实际运用中，等距节点更为常见。设
(),0,1,2,,,i b a
h x a ih i n n
-=
=+= x a th =+， 于是有
)(!
)1)...(1(...!2)1(!1)(00200x R y n n t t t y t t y t y x f n n
+?+--++?-+?+=
所以
'000(1)0
(1)...(1)(){...}1!!
()()
(1)!n n
i t i n n n
j j i
d t t t t n f x y y y dt n h f i j n ξ=+=≠--+=
+?++?+-+∏ 3.多点数值微分公式
由于高阶插值的不稳定性，实际应用时多采用n=1,2,4的两点、三点和五点等多点插值型求导公式。
(1) 两点公式(n=1)
??
???=-=-
=-=)
(2
)()
(1
)()(2)()(1)(''1'011''
'0'010'
ξξf h
x R y y h x f f h x R y y h x f (7) （2）三点公式(n=2)
???
?
??
???=+-=-=+-==-+-=)(3)()34(21)()(6)()(21)()(3)()43(21)()3(20'2
102')3(21'
201')3(20'
2100'ξξξf h x R y y y h x f f h x R y y h x f f h x R y y y h x f （8） (3)五点公式（n =4）
?????
?????
?
??=+-+-=-=++-+-=-=-+-=-=+-+--==-+-+-=)
(5
)()
254836163(121)()
(20)()310186(121)()
(30)()88(121
)()
(20
)()
618103(121)()
(5
)()316364825(121)()5(40'
432104')5(43'
432103'
)5(42'
43102')5(41'
4
32101')
5(40'
432100'ξξξξξf h x R y y y y y h x f f h x R y y y y y h x f f h x R y y y y h x f f h x R y y y y y h x f f h x R y y y y y h x f （9） 例2．设()ln f x x =，取0.05h =，分别用三点公式和五点公式计算()'2f 的近似值。
解:由式(8)有
499802861.0))10.2()05.2(4)2(3(05
.021
)2('=-+-?=
f f f f 500104205.0))05.2()95.1((05
.021
)2('=+-?=f f f 499779376.0))2(3)95.1(4)90.1((05
.021
)2('=+-?=
f f f f 由式(9)有
499999843.0))10.2()05.2(6)95.1(8)90.1((05
.0121
)2('=-+-?=
f f f f f
与真值()'20.5f =相比，三点公式已有相当满意精度，而五点公式的结果是十分满意的。 4.数值微分的隐式格式
前述的Taylor 展式法和数值微分格式均称为显式格式，即直接由已知的
()(0,1,2,,)i f x i n =，经过适当的算术四则运算，立即可得()'i f x 的近似值。显式格式
优点是计算方便，工作量小，缺点是数值不稳定。为克服后一缺点，隐式格式常常具有数值稳定性。
数值微分的隐式格式建立方法常用通过Taylor 展开式方法或数值积分方法等不同途径。
首先我们用Taylor 展开式方法来推导数值微分的隐式格式。
由式(5)和式(6)，我们用k x 代替0x 得：
)()(6
)]()([21)(4'
''2'h O x f h h x f h x f h x f k k k k +---+=
)
(12
)]()(2)([1)()4(22'
'ξf h h x f x f h x f h x f k k k k --+-+=
所以也有
)
(12
)]()(2)([1)()5(2'
''2'
''ξf h h x f x f h x f h x f k k k k --+-+= 将最后()'"
k f
x 表达式代入()'k f x 表达式可得
)()]()(2)([6
1
)]()([21)(4''''h O h x f x f h x f h x f h x f h x f k k k k k k +-+-+---+=
略去误差项()4O h ，并且k m 表示()'k f x 的近似值，则有
),...,2,1()]()([3
41111n k x f x f h
m m m k k k k k =-=
++-++- 同样，用数值积分方法也可推导数值微分的隐式格式。
),...,2,1()(')(1
1
1
1
n k dx
x f dx x f k k k k x x k x x ==?
?
+-+-
根据simpson 公式)]()2
(4)([6)(b f b
a f a f a
b dx x f b
a +++-=
?有
),...,2,1(],4[62)()(1111n k m m m h
x f x f k k k k k =++=-+--+
),...,2,1()]()([3
41111n k x f x f h
m m m k k k k k =-=++-++-
上式是关于1n +个未知量01,,n m m m 的1n -个方程，如果已知0m =
()()''0,n n f x m f x =，记)]()([3
11-+-=k k k x f x f h
d ，故有方程组：
???
???
???
?
??????????--=??????????????
???
???????????????????????----n n n n n m d d d d
m d m m m m m 123201123214114114114114 （10） 这就是求121,,,n m m m - 的线性方程组。由于系数矩阵是严格对角占优的三对角矩阵，因此非奇异，解存在且唯一，可由追赶法求解，且数值稳定。
例3．设函数()ln f x x =在节点 1.5,1.6,1.7,1.8,1.9,x =上的函数值，及
()()''1.5, 2.0f f 如表1所示。试用数值微分隐式格式式(10)和式(11)求出相应节点上一阶
导数的近似值。
表1 例3节点数据
解：由式(10)有
?
????
?
??????--=????????????????????????500000000.03.1608155433676905.353349105.3666666667.075489429.341141141144321m m m m 解方程组得：
1m =0.62499828611483
2m =0.58823447854067 3m =0.55555484972249
4m =0.52631517256938
与()'
1
f
x x
=
在各相应节点上数值相比，上述结果约有五位有效数字，精度较高。h 越小，精度越高。 附追赶法程序： %追赶法
%定义三对角矩阵A的各组成单元。方程为Ax=d
% a为主对角线下面的次对角线，b为主对角线，c为主对角线上面的次对角线，d为等式右端向量。
% A=[4 1 0 0
% 1 4 1 0
% 0 1 4 1
% 0 0 1 4]
function x=zhuiganfa(a,b,c,d)
a=[0 1 1 1];c=[1 1 1];b=[4 4 4 4];
d=[3.75489429-0.666666667,3.53349105,3.33676905,3.16081554-0.5];
n=length(b);
u0=0;y0=0;a(1)=0;
L(1)=b(1)-a(1)*u0;
y(1)=(d(1)-y0*a(1))/L(1);
u(1)=c(1)/L(1);
for i=2:(n-1)
L(i)=b(i)-a(i)*u(i-1);
y(i)=(d(i)-y(i-1)*a(i))/L(i);
u(i)=c(i)/L(i);
end
L(n)=b(n)-a(n)*u(n-1);
y(n)=(d(n)-y(n-1)*a(n))/L(n);
x(n)=y(n);
for i=(n-1):-1:1
x(i)=y(i)-u(i)*x(i+1);
end
运行结果：
ans =
0.62499828611483 0.58823447854067 0.55555484972249 0.52631517256938
>> x=1.5:0.1:2.0;
y=log(x);
dy=[1/1.5 0.62499828611483 0.58823447854067 0.55555484972249 0.52631517256938 1/2.0];
plot(x,y,'b-')
hold on
plot(x,1./x,'r-o')
hold on
plot(x,dy,'g-*')
作图为：
后序：除了本文介绍的常用方法，还有其它许多比较成熟的方法计算数值微分，可以提高稳定近似导数的收敛率。
参考文献
[1]施吉林,刘淑珍,陈桂芝.计算机数值方法(第三版)[M]. 高等教育出版社.2009
[2]吴天毅. 数值积分中点公式的改进[J]. 天津理工大学学报, 2007, (03) .
[3]宫浩. 数值微分算法研究及应用[D]黑龙江大学, 2009
[4]尹秀玲,闫立梅. 一种数值微分方法[J]德州学院学报, 2007,(02)
数值微分的计算方法
数值微分的计算方法 内容摘要 求解数值微分问题,就是通过测量函数在一些离散点上的值,求得函数的近似导数。本文就所学知识，归纳性地介绍了几种常用的数值微分计算方法。并举例说明计算，实验结果表明了方法的有效性。 关键词 数值微分 Taylor 展开式 Lagrange 插值 三对角矩阵 引言：数值微分即根据函数在一些离散点的函数值，推算它在某点的导数或高阶导数的近似值的方法。常见的可以用一个能够近似代替该函数的较简单的可微函数（如多项式或样条函数等）的相应导数作为能求导数的近似值，由此也可导出多点数值微分计算公式。当函数可微性不太好时，利用样条插值进行数值微分要比多项式插值更适宜。 1.Taylor 展开式方法 理论基础：Taylor 展开式 ()()()() ()() ()()()00000002 2! ! n n x x x x f x f x x x f x f x f x n --'''=+-+ ++ + 我们借助Taylor 展开式，可以构造函数()f x 在点0x x =的一阶导数和二阶导数的数值微分公式。取步长0h >则 ),() (2 )()()(0011' '20' 00h x x f h x hf x f h x f +∈++=+ξξ （1） 所以 ),() (2 )()()(0011' '000'h x x f h h x f h x f x f +∈--+= ξξ （2） 同理 ),() (2 )()()(0022' '20' 00x h x f h x hf x f h x f -∈+-=-ξξ （3） ),() (2 )()()(0022' '000'x h x f h h h x f x f x f -∈+--= ξξ （4） 式(2)和式(4)是计算()' 0f x 的数值微分公式，其截断误差为()O h ，为提高精度，将 Taylor 展开式多写几项 ),() (24 )(6)(2)()()(0011) 4(40'''30''20' 00h x x f h x f h x f h x hf x f h x f +∈++++=+ξξ ),() (24 )(6)(2)()()(0022) 4(40'''30''20' 00x h x f h x f h x f h x hf x f h x f -∈+-+-=-ξξ 两式相减得 )()(6 2)()()(40' ''2000' h O x f h h h x f h x f x f +---+= （5） 上式为计算)(0'x f 的微分公式，其截断误差为O(h 2 ),比式（2）和(4)精度高。 两式相加，如果],[)(00) 4(h x h x C x f +-∈，则有
081数值计算方法—常微分方程(组)
科学计算—理论、方法 及其基于MATLAB 的程序实现与分析 微分方程（组）数值解法 §1 常微分方程初值问题的数值解法 微分方程（组）是科学研究和工程应用中最常用的数学模型之一。如揭示质点运动规律的Newton 第二定律： ()()()?????'='==0 00022x t x x t x t F dt x d m （１） 和刻画回路电流或电压变化规律的基尔霍夫回路定律等，但是，只有一些简单的和特殊的常微分方程及常微分方程组，可以求得用公式给出的所谓“解析解”或“公式解”，如一阶线性微分方程的初值问题： () ()0 0y y t f ay dt dy =+= （２） 的解为： ()()()τττd f e y e t y t t a at ?-+=00 （３） 但是，绝大多数在实际中遇到的常微分方程和常微分方程组得不到“解析解”，因此，基于如下的事实：
１、绝大多数的常微分方程和常微分方程组得不到（有限形式的）解析解； ２、实际应用中往往只需要知道常微分方程（组）的解在（人们所关心的）某些点处的函数值（可以是满足一定精度要求的近似值）； 如果只需要常微分方程（组）的解在某些点处的函数值，则没有必要非得通过求得公式解，然后再计算出函数值不可，事实上，我们可以采用下面将介绍的常微分方程（组）的初值问题的数值解法，就可以达到这一目的。 一般的一阶常微分方程（组）的初值问题是指如下的一阶常微分方程（组）的定解问题： ()()0 00,y t y t t t y t F dt dy f =≤≤= （７） 其中 ()()()()???? ?? ? ??=t y t y t y t y n 21 （８） ()()()()???? ?? ? ??=y t f y t f y t f y t F n ,,,,21 （９） 常微分方程（组）的初值问题通常是对一动态过程（动态系统、动力系统）演化规律的描述，求解常微分方程（组）的初值问题就是要了解和掌握动态过程演化规律。 §１.1 常微分方程（组）的Cauch 问题数值解法概论
数值微分
数值微分 数值微分(numerical differentiation) 根据函数在一些离散点的函数值，推算它在某点的导数或高阶导数的近似值的方法。通常用差商代替微商，或者用一个能够近似代替该函数的较简单的可微函数（如多项式或样条函数等）的相应导数作为能求导数的近似值。例如一些常用的数值微分公式(如两点公式、三点公式等)就是在等距步长情形下用插值多项式的导数作为近似值的。此外，还可以采用待定系数法建立各阶导数的数值微分公式，并且用外推技术来提高所求近似值的精确度。当函数可微性不太好时，利用样条插值进行数值微分要比多项式插值更适宜。如果离散点上的数据有不容忽视的随机误差，应该用曲线拟合代替函数插值，然后用拟合曲线的导数作为所求导数的近似值，这种做法可以起到减少随机误差的作用。数值微分公式还是微分方程数值解法的重要依据。 7.1 数值微分 7.1.1 差商与数值微分 当函数是以离散点列给出时，当函数的表达式过于复杂时，常用数值微分近似计算 的导数。在微积分中，导数表示函数在某点上的瞬时变化率，它是平均变化率的极限；在几何上可解释为曲线的斜率；在物理上可解释为物体变化的速率。 以下是导数的三种定义形式： （7.1） 在微积分中，用差商的极限定义导数；在数值计算中返璞归真，导数取用差商（平均变化率）作为其近似值。 最简单的计算数值微分的方法是用函数的差商近似函数的导数，即取极限的近似值。下面是与式（7.1）相应的三种差商形式的数值微分公式以及相应的截断误差。 向前差商 用向前差商（平均变化率）近似导数有： （7.2）
其中的位置在的前面，因此称为向前差商。同理可得向后差商、中心差商的定义。 由泰勒展开 得向前差商的截断误差： 向后差商 用向后差商近似导数有：（7.3） 与计算向前差商的方法类似，由泰勒展开得向后差商的截断误差： 中心差商 用中心差商（平均变化率）近似导数有： （7.4） 由泰勒展开 得中心差商的截断误差：
微积分的基本运算
第4章微积分的基本运算 本章学习的主要目的： 1．复习高等数学中有关函数极限、导数、不定积分、定积分、二重积分、级数、方程近似求解、常微分方程求解的相关知识. 2．通过作图和计算加深对数学概念：极限、导数、积分的理解. 3．学会用MatLab软件进行有关函数极限、导数、不定积分、级数、常微分方程求解的符号运算； 4．了解数值积分理论,学会用MatLab软件进行数值积分;会用级数进行近似计算. 1 有关函数极限计算的MatLab命令 (1)limit(F,x,a) 执行后返回函数F在符号变量x趋于a的极限 (2)limit(F,a) 执行后返回函数F在符号变量findsym(F)趋于a的极限 (3)limit(F) 执行后返回函数F在符号变量findsym(F)趋于0的极限 52
53 (4)limit(F,x,a,’left’) 执行后返回函数F 在符号变量x 趋于a 的左极限 (5)limit(F,x,a,’right’) 执行后返回函数F 在符号变量x 趋于a 的右极限 注:使用命令limit 前,要用syms 做相应符号变量说明. 例7 求下列极限 (1)42 20 x cos lim x e x x -→- 在MatLab 的命令窗口输入: syms x limit((cos(x)-exp(-x^2/2))/x^4,x,0) 运行结果为 ans =-1/12 理论上用洛必达法则或泰勒公式计算该极限: 方法1 =-+-=---=-- - →- →-→2 2 222 20 x 3 22 x 4 2 20 x 12cos lim 4) (sin lim cos lim x x e e x x x e x x e x x x x x 12112112)2(2 lim 1211cos lim 222 220x 2 2 22220 x -=--+=--++-- →- - →x x x e x x x x x e e x 方法2 4 42 224420x 4 2 20 x ))(2) 2()2(1()(!421lim cos lim x x o x x x o x x x e x x +-+---++-=-→- →
数值计算_第7章 数值微分和数值积分
第7章数值微分和数值积分 7.1 数值微分 7.1.1 差商与数值微分 当函数是以离散点列给出时，当函数的表达式过于复杂时，常用数值微分近似计算的导数。在微积分中，导数表示函数在某点上的瞬时变化率，它是平均变化率的极限；在几何上可解释为曲线的斜率；在物理上可解释为物体变化的速率。 以下是导数的三种定义形式： （7.1） 在微积分中，用差商的极限定义导数；在数值计算中返璞归真，导数取用差商（平均变化率）作为其近似值。 最简单的计算数值微分的方法是用函数的差商近似函数的导数，即取极限的近似值。下面是与式（7.1）相应的三种差商形式的数值微分公式以及相应的截断误差。 向前差商 用向前差商（平均变化率）近似导数有： （7.2） 其中的位置在的前面，因此称为向前差商。同理可得向后差商、中心差商的定义。 由泰勒展开
得向前差商的截断误差： 向后差商 用向后差商近似导数有： （7.3） 与计算向前差商的方法类似，由泰勒展开得向后差商的截断误差： 中心差商 用中心差商（平均变化率）近似导数有： （7.4） 由泰勒展开
得中心差商的截断误差： 差商的几何意义 微积分中的极限定义，表示在处切线的斜率，即图7.1中直线的斜率；差商表示过和 两点直线的斜率，是一条过的割线。可见数值微分是用近似值内接弦的斜率代替准确值切线的斜率。 图7.1 微商与差商示意图 例7.1给出下列数据，计算，
解：（5.07－5.06）/（0.04－0.02）= 0.5 （5.05－5.07）/（0.08－0.04）= －0.5 （5.05－5.055）/（0.08－0.10）= 0.25 （(0.10) －(0.06)）/（0.10－0.06）= 18.75 设定最佳步长 在计算数值导数时，它的误差由截断误差和舍入差两部分组成。用差商或插值公式近似导数产生截断误差，由原始值的数值近似产生舍入误差。在差商计算中，从截断误差的逼近值的角度看，越小，则误差也越小；但是太小的会带来较大的舍入误差。怎样 选择最佳步长，使截断误差与舍入误差之和最小呢？ 一般对计算导数的近似公式进行分析可得到误差的表示式，以中心差商为例，截断误差不超过 而舍入误差可用量估计（证明略），其中是函数的原始值的绝对误差限，总误 差为 当时，总误差达到最小值，即 （*） 可以看到用误差的表达式确定步长，难度较大，难以实际操作。
计算物理数值微分与数值积分编程
%离散拉普拉斯算符（del2） [x,y]=meshgrid(-2:2,-3:3); U=x.*x+y.*y V=4*del2(U) %梯形积分 a=[1, 5, 7, 2, 3]; trapz(a) b =[0.95013,0.48598;0.23114,0.8913;0.60684,0.7621]; trapz(b) %默认对列向量积分 trapz(b,2) %对指定的行向量积分 %符号变量积分 Q=int('sin(x)/x',0,pi) %正弦积分函数 eval(Q) %用于简化结果 %函数积分 f=inline('exp(x.^(-3))','x'); format long quad(f,0.5,0.6) quadl(f,0.5,0.6) %精度更高的积分 f=@(x)exp(x.^(-3)) %匿名函数 quad(f,0.5,0.6) function f=sinc(x) %函数文件 if x==0 f=1; else f=sin(x)./x; end >>quad(@sinc,0,pi) %含有参数的函数积分应该用函数文件； %方法1 function f=betaf(t,z,w) %函数文件中变量必须放在前面，参数放在后面f=t.^(z-1).*(1-t).^(w-1); >>z=8/3;w=10/3;tol=1e-6; >>quad(@betaf,0,1,tol,[],z,w) %方法2 F=inline('t.^(z-1).*(1-t).^(w-1)','t','z','w'); z=8/3;w=10/3;tol=1e-6;
Q=quad(F,0,1,tol,[],z,w) %方法3 function chans z=8/3;w=10/3;tol=1e-6; quad(@(t)betaf(t,z,w),0,1) function f=betaf(t,z,w) f=t.^(z-1).*(1-t).^(w-1); %函数的三重积分 Q=triplequad(@(x,y,z)(y*sin(x)+z*cos(x)),0,pi,0,1,-1,1) function f=integrnd(x,y,z) %也可编写函数文件 f=y*sin(x)+z*cos(x); >>Q=triplequad(@integrnd,0,pi,0,1,-1,1) %下面总结下： %(1)quad：采用自适应变步长simpson方法，速度和精度都是最差的，建议不要使用 %(2)quad8：使用8阶Newton-Cotes算法，精度和速度均优于quad，但在目前版本下已被取消 %(3)quadl：采用lobbato算法，精度和速度均较好，建议全部使用该函数 %(4)quadg：NIT(数值积分)工具箱函数，效率最高，但该工具箱需要另外下载 %(5)quadv：quad的矢量化函数，可以同时计算多个积分 %(6)quadgk：很有用的函数，功能在Matlab中最强大 %(7)quad2dggen：一般区域二重积分，效率很好，需要NIT支持 %(8)dblquad：长方形区域二重积分 %(9)triplequadL：长方体区域三重积分 %(10)quadndg：超维长方体区域积分，需要NIT支持 例1：磁场的二维图 clear all a=0.35; the=0:pi/20:2*pi; y= -1:0.04:1; z= -1:0.04:1; [Y,Z,T]=meshgrid(y,z,the); r=sqrt((a*cos(T)).^2+(Y-a*sin(T)).^2+Z.^2); r3=r.^3; dby =a*Z.*sin(T)./r3; by=pi/40*trapz(dby,3); %B的y分量 dbz =a*(a-Y.*sin(T))./r3; bz=pi/40*trapz(dbz,3); %B的X分量 figure(1) [bSY,bSZ]=meshgrid([0:0.05:0.2],0); h1=streamline(Y(:,:,1),Z(:,:,1),by,bz,bSY,bSZ,[0.1,1000]);
微积分计算公式
§3-6 常用积分公式表·例题和点评 ⑴ d k x kx c =+? (k 为常数) ⑵1 1 d (1)1 x x x c μ μμμ+≠-= ++? 特别， 2 1 1d x c x x =- +?, 3 223 x x c = +? , x c =? ⑶ 1 d ln ||x x c x =+? ⑷d ln x x a a x c a = +?, 特别， e d e x x x c =+? ⑸sin d cos x x x c =-+? ⑹cos d sin x x x c =+? ⑺ 2 2 1 d csc d cot sin x x x x c x ==-+?? ⑻ 2 2 1 d sec d tan cos x x x x c x ==+?? ⑼arcsin (0)x x c a a =+>?,特别，arcsin x x c =+? ⑽2 2 1 1d arctan (0)x x c a a a a x = +>+?,特别， 21 d arctan 1x x c x =++? ⑾2 2 1 1d ln (0)2a x x c a a a x a x += +>--? 或 2 2 1 1d ln (0)2x a x c a a x a x a -= +>+-? ⑿ tan d ln cos x x x c =-+? ⒀cot d ln sin x x x c =+? ⒁ln csc cot 1csc d d ln tan sin 2x x c x x x x c x ?-+? = =?+?? ? ? ⒂πln sec tan 1 sec d d ln tan cos 24x x c x x x x c x ?++?= =?? ?++ ?????? ?
郑州大学研究生课程数值分析复习---第八章 常微分方程数值解法
郑州大学研究生课程（2012-2013学年第一学期）数值分析 Numerical Analysis 习题课 第八章常微分方程数值解法
待求解的问题：一阶常微分方程的初值问题/* Initial-Value Problem */: ?????=∈=0 )(] ,[),(y a y b a x y x f dx dy 解的存在唯一性（“常微分方程”理论）：只要f (x , y ) 在[a , b ] ×R 1 上连续，且关于y 满足Lipschitz 条件，即存在与x , y 无关的常数L 使 对任意定义在[a , b ] 上的y 1(x ) 和y 2(x ) 都成立，则上述IVP 存在唯一解。 1212|(,)(,)||| f x y f x y L y y ?≤?一、要点回顾
§8.2 欧拉(Euler)法 通常取(常数),则Euler 法的计算格式 h h x x i i i ==?+1?? ?=+=+) (),(001x y y y x hf y y i i i i i =0,1,…,n ( 8.2 )
§8.2 欧拉(Euler)法(1) 用差商近似导数 )) (,()()()()(1n n n n n n x y x hf x y x y h x y x y +=′+≈+?? ?=+=+) (),(01a y y y x hf y y n n n n 差分方程初值问题向前Euler 方法h x y x y x y n n n ) ()()(1?≈ ′+)) (,() ()(1n n n n x y x f h x y x y ≈?+))(,()(n n n x y x f x y =′
微积分的数值计算方法
第七章 微积分的数值计算方法 7.1 微积分计算存在的问题/数值积分的基本概念 1. 微分计算问题 求函数的导数(微分)，原则上没有问题。当然，这是指所求函数为连续形式且导数存在的情形。但如果函数一表格形式给出，要求函数在某点的导数值；或者是希望某点的导数值只用其附近离散点上的函数值近似地表示，这就是新问题了，它称为微分的数值计算，或称为数值微分。 2．定积分计算问题 计算函数f 在],[b a 上的定积分 dx x f I b a ?= )( 当被积函数f 的原函数能用有限形式)(x F 给出时，可用积分基本公式来计算： )()()(a F b F dx x f I b a -==? 然而，问题在于：① f 的原函数或者很难找到，或者根本不存在；②f 可能给出一个函数表；③仅仅知道f 是某个无穷级数的和或某个微分方程的解等等。这就迫使人们不得不寻求定积分的近似计算，也称数值积分。 3．数值积分的基本形式 数值积分的基本做法是构造形式如下的近似公式 ∑?=≈n k k k b a x f A dx x f 0 )()( (7.1.1) 或记成 ∑?=+=n k n k k b a f R x f A dx x f 0 ][)()( (7.1.2) ∑==n k k k x f A I 0 * )( 和 ][f R n 分别成为],[b a 上的f 的数值求积公式及其 余项(截断误差)，k x 和k A ),,1,0(n k =分别称为求积节点和求积系数(求积系数与被积函数无关)。 这种求积公式的特点是把求积过(极限过程)程转化为乘法与加法的代数运算。构造这种求积公式需要做的工作是：确定节点k x 及系数 k A ),,1,0(n k =，估计余项][f R n 以及讨论* I 的算法设计及其数值稳定 性。 4．插值型求积公式 如何构造求积公式呢？基本的技术是用被积函数f 的Lagrange 插值多项式 )(x L n 近似代替f ，也即对],[b a 上指定的1+n 个节点
数值分析_第五章_常微分方程数值解法
图５畅２ 令珔h ＝h λ，则y n ＋１＝１＋珔 h ＋１２珔h ２ ＋１６珔h ３＋１２４ 珔 h ４y n ．由此可知，绝对稳定性区域在珔h ＝h λ复平面上满足 ｜１＋珔 h ＋１２珔h ２＋１６珔h ３＋１２４珔h ４ ｜≤１的区域，也就是由曲线 １＋珔h ＋ １２珔h ２＋１６珔h ３＋１２４ 珔h ４＝e i θ 所围成的区域．如图５畅２所示． 例22 用Euler 法求解 y ′＝－５y ＋x ，y （x ０）＝y ０， x ０≤x ≤X ． 从绝对稳定性考虑，对步长h 有何限制？ 解 对于模型方程y ′＝λy （λ＜０为实数）这里λ＝抄f 抄y ＝－５．由 ｜１＋h λ｜＝｜１－５h ｜＜１ 得到对h 的限制为：０＜h ＜０畅４． 四、习题 １畅取步长h ＝０畅２，用Euler 法解初值问题 y ′＝－y －x y ２ ， y （０）＝ １． （０≤x ≤０畅６）， ２畅用梯形公式解初值问题 y ′＝８－３y ， （１≤x ≤２），
取步长h＝０畅２，小数点后至少保留５位． ３畅用改进的Euler公式计算初值问题 y′＝１x y－１x y２， y（１）＝０畅５， １＜x＜１畅５， 取步长h＝０畅１，并与精确解y（x）＝ x １＋x比较． ４畅写出用梯形格式的迭代算法求解初值问题 y′＋y＝０， y（０）＝１ 的计算公式，取步长h＝０畅１，并求y（０畅２）的近似值，要求迭代误差不超过１０－５． ５畅写出用四阶经典Runge唱Kutta法求解初值问题 y′＝８－３y， y（０）＝２ 的计算公式，取步长h＝０畅２，并计算y（０畅４）的近似值，小数点后至少保留４位． ６畅证明公式 y n＋１＝y n＋h９（２K１＋３K２＋４K３）． K１＝f（x n，y n）， K２＝f x n＋h２，y n＋h２K１， K３＝f x n＋３４h，y n＋３４h K２， 至少是三阶方法． ７畅试构造形如 y n＋１＝α（y n＋y n－１）＋h（β０f n＋β１f n－１）
数值积分与微分方法
数值积分与微分 摘要 本文首先列举了一些常用的数值求积方法，一是插值型求积公式，以N e w t o n C o t e s -公式为代表，并分析了复合型的Newton Cotes -公式；另一个是Gauss Ledendre -求积公式，并给出几个常用的Gauss Ledendre -求积公式。其次，本文对数值微分方法进行分析，主要是差分型数值微分和插值型数值微分，都给出了几种常用的微分方法。然后，本文比较了数值积分与微分的关系，发现数值积分与微分都与插值或拟合密不可分。 本文在每个方法时都分析了误差余项，并且在最后都给出了MATLAB 的调用程序。 关键词：插值型积分Gauss Ledendre -差分数值微分插值型数值微分 MATLAB
一、常用的积分方法 计算积分时，根据Newton Leibniz -公式， ()()()b a f x dx F b F a =-? 但如果碰到以下几种情况： 1）被积函数以一组数据形式表示； 2）被积函数过于特殊或者原函数无法用初等函数表示 3）原函数十分复杂难以计算 这些现象中，Newton Leibniz -公式很难发挥作用，只能建立积分的近似计算方法，数值积分是常用的近似计算的方法。 1.1 插值型积分公式 积分中的一个常用方法是利用插值多项式来构造数值求积公式，具体的步骤如下： 在积分区间上[,]a b 上取一组节点：01201,,,,()n n x x x x a x x x b ≤<<≤ 。已知()k x f 的函数值，作()x f 的n 次插值多项式，则 (1) ()10()()()()() (1)!n n x n n k k n k f f L x R x f x l x w x n ++==+=++∑ 其中，()k l x 为n 次插值基函数，则得 (1)+10()(()())1 =[()]()[()](1)!b b n n a a n b b n k k n a a k f x dx L x R x dx l x dx f x f x w x dx n ξ+==+++? ?∑??（） 公式写成一般形式： ()()[]n b k k n a k f x dx A f x R f ==+∑? 其中， 01100110 ()()()() ()()()()()b b k k k k a a k k k k k k x x x x x x x x A l x dx dx x x x x x x x x -+-+----==----?? (1)+11 [][()](1)!b n n n a R f f x w x dx n ξ+=+?（） 显然，当被积函数f 为次数小于等于n 的多项式时，其相应的插值型求积公式为准确公式，即： ()() n b k k a k f x dx A f x ==∑? 1.1.1 求积公式的代数精度 定义：求积公式对于任何次数不大于m 的代数多项式()f x 均精确成立，而对于 1()m f x x +=不精确成立，则称求积公式具有m 次代数精度。 定理：含有1n +个节点(0,1,,)k x k n = 的插值型求积公式的代数精度至少为n 。
偏微分方程数值解法
《偏微分方程数值解法》 课程设计 题目: 六点对称差分格式解热传导方程的初边 值问题 姓名: 王晓霜 学院: 理学院 专业: 信息与计算科学 班级: ０91１01２ 学号: ０９１10１218 指导老师：翟方曼 ２０12年12月14
日 一、题目 用六点对称差分格式计算如下热传导方程的初边值问题 222122,01,01(,0),01 (0,),(1,),01x t t u u x t t x u x e x u t e u t e t +???=<<<≤?????=≤≤??==≤≤??? 已知其精确解为 2(,)x t u x t e += 二、理论 1.考虑的问题 考虑一维模型热传导方程 （1.１) )(22x f x u a t u +??=??,T t ≤<0 其中a 为常数。)(x f 是给定的连续函数。（1.１)的定解问题分两类: 第一,初值问题（Cauch ｙ 问题)：求足够光滑的函数()t x u ,,满足方程（１.1）和初始条件： (1.２) ()()x x u ?=0,， ∞<<∞-x 第二,初边值问题（也称混合问题）:求足够光滑的函数()t x u ,，满足方程(１.１)和初始条件： ()13.1 ()()x x u ?=0,， l x l <<- 及边值条件 ()23.1 ()()0,,0==t l u t u ， T t ≤≤0 假定()x f 和()x ?在相应的区域光滑,并且于()0,0，()0,l 两点满足相容条件,则上述问题有唯一的充分光滑的解。 现在考虑边值问题（１.1),（1.３)的差分逼近 取 N l h = 为空间步长,M T =τ为时间步长,其中N ,M 是自然数,
微积分计算方法
学号 1330101009 毕业论文 对概率积分解法的研究和讨论 院（系）名称：书信学院 专业名称：数学教育 学生姓名：李建鹏 指导教师：杜争光 二○一五年
摘要：文章给出了计算概率积分 2 x e dx ∞- -∞ ?的几种简便的计算方法；对以 后概率积分的研究和应用具有较好的帮助。 关键词：格林公式；奥高公式；重积分；含参变量 概率积分 2 x e dx ∞- -∞ ?是重要的积分之一，再数理方程、概率论等方面经 常用到，且有广泛的应用。而关于这个积分值的计算问题，有不少人讨论过，大多数方法要用到较深的预备知识，本文给出了几种所需预备知识而又简便的计算方法。
目录 方法一：二重积分法 (1) 方法二：三重积分法 (1) 方法三：线积分法 (2) 方法四:面积分法 (3) 方法五:含参变量的无穷积分法 (4) 方法六：二重积分证明法 (6) 参考文献： (8) 致谢: (9)
对概率积分2 x e dx ∞ --∞ ? 解法的研究和讨论 概率积分 2 x e dx ∞ --∞ ? 是重要的积分之一，再数理方程、概率论等方面经常用 到，且有广泛的应用。而关于这个积分值的计算问题，有不少人讨论过，大多数方法要用到较深的预备知识，本文给出了几种所需预备知识而又简便的计算方法。 方法一：二重积分法 现有连续函数 22() (,)x y f x y e -+=在正方形区域:(;)D a x a a y a -≤≤-≤≤； 圆域2 2 2 1:()R x y a +≤；圆域：2 222 :(2)R x y a +≤上的二重积分分别为12,,I I I ， 即： 22 22 2 () () 2 ()a a a x y x y x a a a D I e d x d y d x e d y e d x -+-+----===????? 22 22 1 2() 10 .(1) a x y r a R I e d x d y d r e d r e πθπ-+--===-???? 2222 2 22() 220 .(1) a x y r a R I e dxdy d r e dr e πθπ-+--===-???? （用极坐标） 同时又因：1 2I I I ≤≤,故有 12 lim lim lim a a a I I I →∞ →∞ →∞ ≤≤，即有2 2 lim()a t a a e dt π--→∞ =? ，从而 2 x e dx π ∞ --∞ =? [] 4 方法二：三重积分法 首先我们把旋转体的体积概念推广到积分限无穷的情况。再设XOZ 平面上的曲线2 x Z e -=绕Z 轴旋转一周得到的曲面22() x y Z e -+=与平面XOY 围成 的体V 。显然，一方面，该体的体积 22() 2 2 () x y e x v V dxdydz dx dy dz e dx -+∞ ∞ ∞ --∞ -∞ -∞ ===?????? ? 另一方面，根据旋转体的体积公式有：
常微分方程数值解
第四章常微分方程数值解 [课时安排]6学时 [教学课型]理论课 [教学目的和要求] 了解常微分方程初值问题数值解法的一些基本概念，如单步法和多步法，显式和隐式，方法的阶数，整体截断误差和局部截断误差的区别和关系等；掌握一阶常微分方程初值问题的一些常用的数值计算方法，例如欧拉（Euler）方法、改进的欧拉方法、龙贝－库塔（Runge-Kutta）方法、阿达姆斯（Adams）方法等，要注意各方法的特点及有关的理论分析；掌握构造常微分方程数值解的数值积分的构造方法和泰勒展开的构造方法的基本思想，并能具体应用它们导出一些常用的数值计算公式及评估截断误差；熟练掌握龙格－库塔（Ｒ－Ｋ）方法的基本思想，公式的推导，Ｒ－Ｋ公式中系数的确定，特别是能应用“标准四阶Ｒ－Ｋ公式”解题；掌握数值方法的收敛性和稳定性的概念，并能确定给定方法的绝对稳定性区域。[教学重点与难点] 重点：欧拉方法，改进的欧拉方法，龙贝－库塔方法。 难点：R—K方法，预估-校正公式。 [教学内容与过程] 4.1 引言 本章讨论常微分方程初值问题 (4.1.1) 的数值解法，这也是科学与工程计算经常遇到的问题，由于只有很特殊的方程能用解析方法求解，而用计算机求解常微分方程的初值问题都要采用数值方法.通常我们假定(4.1.1)中 f(x,y)对y满足Lipschitz条件，即存在常数L＞0，使对，有 (4.1.2) 则初值问题(4.1.1)的解存在唯一. 假定(4.1.1)的精确解为，求它的数值解就是要在区间上的一组离散点 上求的近似.通常取 ，h称为步长，求(4.1.1)的数值解是按节点的顺序逐步 推进求得.首先，要对方程做离散逼近，求出数值解的公式，再研究公式的局部截
计算方法6_微分方程
习题6 6.1 试用三种方法导出线性二步方法 122+++=n n n hf y y 6.2 用Taylor 展开法求三步四阶方法类，并确定三步四阶显式方法. 6.3 形如 ∑=++=k i k n k j n j f h y 0βα 的k 阶方法称为Gear 方法，试确定一个三步Gear 方法，并给出其截断误差主项。 6.4 试用显式Euler 法及改进的Euler 法 )],(),([2 11n n n n n n n hf y t f y t f h y y +++=++ 6.5 给出线性多步法 ])13()3[(4 )1(212n n n n n f f h y y y +++=--++++αααα 为零稳定的条件，并证明该方法为零稳定时是二阶收敛的. 6.6 给出题(6.5)题中1=α时的公式的绝对稳定域. 6.7 指出Heun 方法 0 0 0 0 1/3 1/3 0 0 2/3 0 2/3 0 1/4 0 3/4 的相容阶，并给出由该方法以步长h 计算初值问题(6.45)的步骤. 6.8 试述刚性问题的基本特征，并给出s 级Runge-Kutta 方法为A -稳定的条件. 6.9 设有???=='00 )(),(y x y y x f y ，试构造形如 )()(11011--++++=n n n n n f f h y y y ββα 的二阶方法，并推导其局部截断误差首项。
6.10设有常微分方程初值问题???=='00 )(),(y x y y x f y 的单步法)],(2),([3 111+++++=n n n n n n y x f y x f h y y ，证明该方法是无条件稳定的。
数值积分与数值微分 编程计算
解：卫星轨道的示意图如右上图所示，,a b 分别是椭圆轨道的长半轴和短半轴，地球位于椭圆的一个焦点处，焦距为c ，地球半径为r ，近地点和远地点与地球表面的距离分别是1s 和2s . 由图中可知，上述数据存在如下关系： 12122,a s s r c a s r =++=-- 由椭圆性质 b =，将12,,s s r 的数据代入以上各式可得7782.5a km =，7721.5b km =. 椭圆的参数方程为： c o s ,s i n x a t y b t == ， (02)t π≤<
根据计算参数方程弧长的公式，椭圆长度可表为如下积分： /2 22221/20 4(sin cos )L a t b t dt π=+? 由于该积分无法求得解析解，下面我们编写MATLAB 程序对其进行数值求解。 s1=439;s2=2384;r=6371; a=(s1+s2)/2+r a = 7.7825e+003 >> c=a-s1-r; >> b=sqrt(a^2-c^2) b = 7.7215e+003 y=inline('sqrt(7782.5^2*sin(t).^2+7721.5^2*cos(t).^2)'); %建立积分内联函数 >> t=0:pi/10:pi/2; y1=y(t); format long >> L1=4*trapz(t,y1) %梯形积分 L1 = 4.870744099902405e+004 >> L2=4*quad(y,0,pi/2,1e-6) %辛普森积分 L2 = 4.870744099903280e+004 求解结果显示：两种方法计算求得的积分结果相当接近，轨道长度约为：4 4.8710km ?.
高数微分公式
初等数学基础知识 一、三角函数 1．公式 同角三角函数间的基本关系式： ·平方关系： sin^2(α)+cos^2(α)=1; tan^2(α)+1=sec^2(α);cot^2(α)+1=csc^2(α)·商的关系： tanα=sinα/cosαcotα=cosα/sinα ·倒数关系： tanα·cotα=1; sinα·cscα=1; cosα·secα=1 三角函数恒等变形公式： ·两角和与差的三角函数： cos(α+β)=cosα·cosβ-sinα·sinβ cos(α-β)=cosα·cosβ+sinα·sinβ sin(α±β)=sinα·cosβ±cosα·sinβ tan(α+β)=(tanα+tanβ)/(1-tanα·tanβ) tan(α-β)=(tanα-tanβ)/(1+tanα·tanβ) 倍角公式： sin(2α)=2sinα·cosα cos(2α)=cos^2(α)-sin^2(α)=2cos^2(α)-1=1-2sin^2(α) tan(2α)=2tanα/[1-tan^2(α)] ·半角公式： sin^2(α/2)=(1-cosα)/2 cos^2(α/2)=(1+cosα)/2 tan^2(α/2)=(1-cosα)/(1+cosα) tan(α/2)=sinα/(1+cosα)=(1-cosα)/sinα ·万能公式： sinα=2tan(α/2)/[1+tan^2(α/2)] cosα=[1-tan^2(α/2)]/[1+tan^2(α/2)] tanα=2tan(α/2)/[1-tan^2(α/2)]
·积化和差公式： sinα·cosβ=(1/2)[sin(α+β)+sin(α-β)] cosα·sinβ=(1/2)[sin(α+β)-sin(α-β)] cosα·cosβ=(1/2)[cos(α+β)+cos(α-β)] sinα·sinβ=-1/2[ cos(α-β)-cos(α+β)] ·和差化积公式： sinα+sinβ=2sin[(α+β)/2]cos[(α-β)/2] sinα-sinβ=2cos[(α+β)/2]sin[(α-β)/2] cosα+cosβ=2cos[(α+β)/2]cos[(α-β)/2] cosα-cosβ=-2sin[(α+β)/2]sin[(α-β)/2] 2 只需记住这两个特殊的直角三角形的边角关系，依照三角函数的定义即可推出上面的三角值。 3诱导公式： 记忆规律： 竖变横不变（奇变偶不变），符号看象限（一全，二正弦割，三切，四余弦割 即第一象限全是正的，第二象限正弦、正割是正的，第三象限正切是正的，第四象限余弦、余割是正的） 1 ο45 2 1 ο45 1 2 ο30 ο60 3
数值微分计算方法实验课件
数值微分计算方法实验 郑发进 2012042020022 【摘要】数值微分(numerical differentiation )根据函数在一些离散点的函数值，推算它在某点的导数或高阶导数的近似值的方法。通常用差商代替微商，或者用一个能够近似代替该函数的较简单的可微函数（如多项式或样条函数等）的相应导数作为能求导数的近似值。例如一些常用的数值微分公式(如两点公式、三点公式等)就是在等距步长情形下用插值多项式的导数作为近似值的。此外，还可以采用待定系数法建立各阶导数的数值微分公式，并且用外推技术来提高所求近似值的精确度。当函数可微性不太好时，利用样条插值进行数值微分要比多项式插值更适宜。如果离散点上的数据有不容忽视的随机误差，应该用曲线拟合代替函数插值，然后用拟合曲线的导数作为所求导数的近似值，这种做法可以起到减少随机误差的作用。数值微分公式还是微分方程数值解法的重要依据。 关 键 词 中心差分公式；理查森外推法；牛顿多项式微分； 一、实验目的 1.通过本次实验熟悉并掌握各种数值微分算法。 2.掌握如何通过程序设计实现数值微分算法，从而更好地解决实际中的问题。 二、实验原理 1. 精度为2()O h 的中心差分公式： 设[]3,f C a b ∈，且[],,,x h x x h a b -+∈，则 ()()() '2f x h f x h f x h +--≈ 而且存在数()[],c c x a b =∈，满足
()()() ()',2trunc f x h f x h f x E f h h +--≈ + 其中 ()()() ()322,6 trunc h f c E f h O h =-= 项(),trunc E f h 称为截断误差。 2. 精度为4()O h 的中心差分公式： 设[]5,f C a b ∈，且[]2,,,,2,x h x h x x h x h a b --++∈，则 ()()() '2f x h f x h f x h +--≈ 而且存在数()[],c c x a b =∈，满足 ()()()()()2882'12f x h f x h f x h f x h f x h -+++--+-≈ 其中 ()()() ()5 44,30 trunc h f c E f h O h == 项(),trunc E f h 称为截断误差。 3.理查森外推法： 利用低阶公式推出高阶求解数值微分的公式，定理如下： 设()0'f x 的两个精度为2()k O h 的近似值分别为()1k D h -和()12k D h -，而且 它们满足 ()()2220112'k k k f x D h c h c h +-=+++ 和 ()()21220112'244k k k k k f x D h c h c h ++-=+++ 这样可得到改进的近似值表达式
计算方法--常微分方程求解实验
实验五 常微分方程求解实验 一、 实验目的 通过本实验学会对给定初值我呢他，用欧拉法、改进欧拉法、四阶龙格-库塔法求数值解和误差，并比较优缺点.对给定刚性微分方程，求其数值解，并与精确解比较，分析计算结果. 二、 实验题目 1. 解初值问题各种方法比较 实验题目：给定初值问题 ?? ???=≤<+=,0)1(, 21,e d d y x x x y x y x 精确解为)e -e (x x y =，按 （1） 欧拉法，步长;1.0,025.0==h h （2） 改进欧拉法，步长;01.0,05.0==h h （3） 四阶标准龙格-库塔法，步长;1.0=h 求在节点)10,,2,1(1.01 =+=k k x k 处的数值解及误差，比较各方法的优缺点. 2. 刚性方程计算 实验题目：给定刚性微分方程 ?? ???=≤<-+-=-,2)0(, 50,600e 8.1199600d d 1.0y x y x y x 其精确解为12e e )(-0.1x -600x -+=x y .任选取一显示方法，取不同的步长求解，并分析计算 结果. 三、 实验原理 1.欧拉格式 由数值微分的向前差商公式可以解决初值问题（6.1）()()? ??=≤≤=0 0', ,,y x y b x a y x f y 中的导数的 数值计算问题： ()()() ()() ,'1h x y x y h x y h x y x y n n n n -= -+≈ + 由此可得 ()()().'1n n n x hy x y x y +≈+
（6.1）实际上给出 ()()()()()().,','n n n x y x f x y x y x f x y =?= 于是有 ()()()().,1n n n n x y x hf x y x y +≈+ 再由()()11,++≈≈n n n n x y y x y y 得 ().,1,0,,111 =+=+++n y x hf y y n n n n （6.2） 递推公式（6.2）称为欧拉格式。 2.改进欧拉格式 先对欧拉格式（6.2）对1+n y 进行计算，并将结果记为1+n y ，再代入（6.7） ()()[]111,,2 +++++ =n n n n n n y x f y x f h y y 可得“预报-校正”形式的差分格式： ()()()[] ?? ? ? ?++=+=++++.,,2,,1111n n n n n n n n n n y x f y x f h y y y x hf y y 公式（6.8）称为改进欧拉格式。 3.四阶经典龙格-库塔格式： ()()()?? ?? ? ? ??? ???? +=???? ??+=???? ??+==++++=++++.,,2,,2,,,,226 31422131212143211 hK y x f K K h y x f K K h y x f K y x f K K K K K h y y n n n n n n n n n n 四、 实验内容 1. 解初值问题各种方法比较 实验程序： function shiyan51 h=0.1; dyfun=inline('y./x+x*exp(x)'); [x,y1]=maeuler1(dyfun,[1,2],0,h); [x,y2]=maeuler(dyfun,[1,2],0,h); [x,y3]=marunge4(dyfun,[1,2],0,h); y=x.*(exp(x)-exp(1));