科学计算—理论、方法
及其基于MATLAB 的实现与分析
数值微积分
§1 数值微分
对于给定的函数()x f y =，如果 １、()x f y =的函数关系式比较复杂时；
２、()x f y =未知，而仅仅知道()x f y =在1+n 个相异点k x ，n k ,,1,0 =处 的函数值k y ；
则希望能用相对简单的计算方法，求得()x f y =导数的（近似）值。
基于上述考虑，选择的方法之一是利用函数()x f y =的插值多项式的导数作为函数()x f y =导数的近似值，例如Lagrange 插值多项式，由于
()()()()∑==≈n
k k k n x f x l x L x f 0 （１）
因而有
()()x L x f n
'≈' （２） 这里需要说明一点的是，尽管()x f 和()x L n 的函数值可能相差不多，但是它们的导数有可能相差很大，如下面的例子
例1： 考虑函数()2
2511
x x f +=
在区间[-1，1]的插值问题，取区间
[-1，1]的11个点k x k ?+-=2.01，10,,1,0 =k ，作函数()2
2511
x x f +=的
10次插值多项式：
()18552.163597.1234338.3819095.4949417.220246810+-+-+-=x x x x x x L n
函数()x f 和插值多项式()x L n 的导数分别为 ()()
2
2
25150x x
dx x df +-=
()x x x x x dx
x dLn 7.334.4936.22883.39594.22093579-+-+-= 对函数()x f 和插值多项式()x L n 及其导数分别比较，结果如图所示：
Derivative_Runge
下面我们基于Taylor 公式
()()()()()()()()()()??
?????++++''+'+=++++1
1
12!
1!!2n n n n n h O h n f h n x f h x f h x f x f h x f ξ (3)
讨论导数的近似计算问题 1 一阶导数的近似计算
令k k x x h -=+1,可得一阶向前有限差商公式(First Forward Finite Divided Difference ):
()()()()()
()()()()()
()()()()?????=-≈
'?+''--='?+''+
'+=+++h O error h x f x f x f h O h x f h x f x f x f h O h x f h x f x f x f k k k k k k k k k k k 1213
2
1!2!
2 (4)
类似地,当1--=k k x x h 时,可得一阶向后有限差商公式(First Backward Finite Divided Difference ):
()()()()()
()()()()()
()()()()?????=-≈
'?+''+-='?+''+
'-=---h O error h x f x f x f h O h x f h x f x f x f h O h x f h x f x f x f k k k k k k k k k k k 1213
2
1!2!
2 (5)
由近似计算公式(4)和(5),可得可得一阶中心有限差商公式(First Centered Finite Divided Difference ):
()()()()()()
()()()()
?????=-≈'
?+-=
'?
+-+-+2112112254h O error h
x f x f x f h O h
x f x f x f k k k k k k (6)
2 二阶导数的近似计算
在Taylor 公式中，用h 2代替h 可得
()()()()()()
3222!
22h O h x f h x f x f x f k k k k +''+
'+=+ (7) 式(7)与式(4)结合可得二阶向前有限差商公式(Second Forward Finite Divided Difference ):
()()()()()
32122h O h x f x f x f x f k k k k ++''+-=-++ (8)
()()()()()
?????
=+-≈
''?++h O error h x f x f x f x f k k k k 2
122 ()()()()()()
?
????
=---≈''?+++h O error h h x f x f h x f x f x f k k k k k 112 (9) 类似地利用
()()()()()()
3222!
22h O h x f h x f x f x f k k k k +''+'-=- ()()()()()
32
1!
2h O h x f h x f x f x f k k k k +''+
'-=- 可得二阶向后有限差商公式(Second Backward Finite Divided Difference ):
()()()()()32122h O h x f x f x f x f k k k k ++''+-=---
()()()()()
?????
=+-≈''?-+h O error h x f x f x f x f k k k k 2
212 (10) 进一步底，由公式
()()()()()()
()()()()()()
4
3
214
3
21!3!2!3!2h O h x f h x f h x f x f x f h O h x f h x f h x f x f x f k k k k k k k k k k +'''-''+'-=+'''+''+
'+=-+ (11) 可得二阶中心有限差商公式(Second Centered Finite Divided Difference ):
()()()()()
?????
=+-≈''-+22
112h
O error h x f x f x f x f k k k k (12)
§2 数值积分
§2.1 Newton-Cotes 数值积分公式
基于数值微分的同样原因，当函数()x f y =在],[b a 上的定积分不易计算时，利用函数()x f y =的Lagrange 插值多项式()()()∑==n
k k k n x f x l x L 0，
有
()()()()k n
k b
a k b
a n b
a
x f dx x l dx x L dx x f I ∑???
=?????
?=≈=0 (13) 由于Lagrange 插值多项式()()()∑==n
k k k n x f x l x L 0
的插值基函数()x l k 只依赖
插值节点k x ，n k ,,1,0 =，所以当k x ，n k ,,1,0 =取定后，()x l k 就是完全确定的多项式函数，令
()?=b
a
k k dx x l A ， n k ,,1,0 = (14)
则由式(13)得到Newton-Cotes 求积公式：
()()()k n
k k b
a
n b a x f A dx x L dx x f I ∑??==≈=0
(15)
特别地，当取插值节点为
kh a x k +=，n k ,,1,0 =，n
a
b h -=
(16) 时有
１） 两点公式（Trapezoidal Rule ）：
?=1n 2
210a
b h A A -=
=
= ()()()()()[]b f a f a
b x f A dx x L dx x f I k n
k k b
a
b
a
+-=
=≈=∑??=2
1 (17) 2） 三点公式（Simpson's 1/3 Rule ）：
3
4,322120h
A h A A a b h n ===?-=
?=， ()()()()()()??
????+???
??++=??
????+??? ??++-=
≈=??b f b a f a f h b f b a f a f a b dx
x L dx x f I b a
b
a
2432462 (18) 3） 四点公式（Simpson's 3/8 Rule ）：
8
9,83332130h
A A h A A a b h n ====?-=
?= ()()()()()()[]()()()()[]b f h a f h a f a f h
b f h a f h a f a f a
b dx
x L dx x f I b a
b
a
+++++=+++++-=
≈=??2338323383 (19) 4） 误差估计：
利用Lagrange 插值的误差公式：
()()()()()()()∏=+-+=-=n
k k n n n x x n f x L x f x R 0
1!1ξ (20) 容易得到Trapezoidal Rule 的误差估计：
()()()
()()()()()
()()()3
2321
121222a b M
a b f dx ab x b a x f dx b x a x f dx
x L dx x f b a b
a b
a
b a
-≤-''=++-''≈--''=
-????ξξξ (21) Simpson's 1/3 Rule 以及Simpson's 3/8 Rule 的误差估计分别为：
()()()()5
445229012901??
?
??-≤??? ??--≈-??
a b M f a b dx x L dx x f b
a b
a
ξ (22) ()()()()5
445338033803??
? ??-≤??? ??--≈-??
a b M f a b dx x L dx x f b
a b
a
ξ (23) 其中()(){}x f M k b a k ]
,[max =.
注１ 数值积分公式(13)表明，数值积分方法最大的优点是将复杂的函数积分转化为相对简单的加、减、乘、除运算；
注２ 从几何的角度看，上述三种求积公式（17）、（18）和（19）（在积分区间],[b a 上）分别用直线
()()b f a
b a
x a f b a b x y --+--=
（24） 抛物线
()()()()()()()()()b f h h a x a x h a f h b x a x a f h h a x h a x y 2
2
2
2---++-------=
（25）
三次立方曲线
()()()()
()()()()()()()()()()()()h a f h h a x h a x a x h a f h h a x h a x a x h a f h h a x h a x a x a f h h a x h a x h a x y 3323323233
3
3
+-----+
+------
+-----+
-------
= （26）
代替曲线()x f y =:
例1 计算积分?=10dx e I x
，并估计误差. 解 1） 用梯形公式计算，
8591409.1)(2011
0≈+-=
≈e e T I .
由于x e x f =)(，所以x
e x
f ='')(，于是，梯形公式的误差
23.012)(max 12110≈=''≤
-=≤≤e
x f T I R x T .
2） 用辛普森公式计算，
7188612.1)4(60112/10
≈++-=
≈e e e S I .
由于
x ）
（
e x
f =)(4，于是，辛普森公式的误差
00094.02880)(max 28801
41
0≈=
≤
-=≤≤e
x f S I R ）
（
x S .
【注】 （1） 当0)()2(=x f 时，梯形求积公式准确成立；当)()
4(x f =0
时，辛普森公式准确成立.即梯形公式对一次多项式准确成立，而辛普森公式对三次多项式准确成立.
（2） 一般地，由余项公式(7.9)知，当)(x f 是n 次多项式时，积分余项为零，从而牛顿-柯特斯求积公式准确成立.
（3） 数值求积方法是一种近似方法，因此，要求求积公式对尽可能多的被积函数)(x f 能准确计算积分值.作为衡量公式逼近好坏的标准之一，下面给出代数精度(the Degree of Precision of the Quadrature Formula)的概念.
§2.1.2 代数精度
【定义1】 如果某求积公式对于次数小于等于m 的多项式能准确求出积分值，而对某个1+m 次多项式就不能准确求出积分，则称该求积公式具有m 次代数精度.
由定义1，具有m 次代数精度的求积公式（7.3）对m
x
x x f ,,,1)( =时精确成立，而对1
)(+=m x x f 不能精确成立.为了构造形如公式（7.3）
的求积公式，通过解方程组
???
????????+-==-==-==∑?∑?∑?=++==n k b a m m m m
k k n k b a k k n
k b a k m a b dx x x A a b xdx x A a b dx A 0110
2
20
,1,2,1
可得求积节点k x 与求积系数k A ，由此求得具有m 次代数精度的求积公式.
不难验证，梯形公式和矩形公式均具有一次代数精度，而辛普森公式具有三次代数精度.对牛顿-柯特斯求积公式，有下面结论.
【定理1】 牛顿-柯特斯求积公式（7.4）至少具有n 次代数精度，而当n 为偶数时，至少具有1+n 次代数精度.
证明 当)(x f 为任何次数不高于n 的多项式时，0)()1(≡+x f n ，所以
0)(=?
b
a
n dx x R ,显然结论成立.
当n 为偶数时，只须对1)(+=n x x f 时的结论验证.因为)!1()()1(+=+n x f n ，由截断误差公式
???∏∏=+=-=-=b
a
b a
n
n
j n n
j j
n
dt
j t h dx x x dx x R 00
20
)()()(，
若令2
n u t +=，则有
?
?
∏-
=+-+=b
a
n n n j n n du
j n u h
dx x R 220
2
)2()(.
注意到∏=-+n
j j n u 0
)2(是u 的奇函数（n 为偶数），因此?=b a n dx x R 0)(，结论
成立.
例2 对于求积公式2
0011
()(1)(2)'(2)
f x dx A f B f B f -=-++?，试确定
100,,B B A ，使此求积公式有最高的代数精度，是几阶的？
§2.2 复化求积公式
由前面分析知道，在整个积分区间],[b a 上用直线或用抛物线代替原曲线()x f y =做积分，尽管计算简单，但精度较差，结果不好，因此，一个自然的想法是：在尽可能保持公式(17)、和(18)计算简单、
易于计算机实现的优点的前提下，寻找提高误差精度的途径和方法。下面给出的“复化求积公式”就是有效的方法之一
受“利用低次多项式分段插值”能够提高逼近精度且计算简单的事实启发，将上述的“在整个积分区间],[b a 上用一个插值多项式代替被积函数”的做法变为“用分段插值多项式代替被积函数”，方法具体如下：
首先，用1+n 个点b x x x a n ==,,,10 将区间],[b a 分成n 个小区间
],[1+k k x x ，1,,1,0-=n k ，然后在每个小区间],[1+k k x x 上用插值多项式代
替被积函数做积分，也即在整个积分区间],[b a 上用分段插值多项式代替被积函数做积分。
如果所取的节点
kh a x k +=，n
a
b h -=
，n k ,,1,0 = 那么在每个小区间],[1+k k x x 上做线性插值时，就得到
１） 复化的梯形求积公式
()()()()()n n k k n k x x b
a
T b f x f a f h dx x f dx x f I k k
=??
?
???++≈==∑∑?
?-=-=+1
11
221 (27)
２） 复化的Simpson's 1/3 求积公式
()()()()()n n k k n k k n k x x b
a
S b f h x f x f a f h dx
x f dx x f I k k
=??
????+??? ??
+++≈==∑∑∑?
?-=-=-=+1
0111
24261 (28) 3） 复化的Simpson's 3/8 求积公式
()()()()()n
n k k n k k n k k n k x x b
a
S b f h x f h x f x f a f h dx
x f dx x f I k k
=??
????+??? ??++??? ??
+++≈==∑∑∑∑?
?-=-=-=-=+101
0111
323332831
(29)
复化求积公式的截断误差
设?=b
a dx x f I )(，则复化梯形公式的截断误差
∑-=''--=''-=-1
2
3)(12)](12[n k k n f h a b f h T I ηη，
其中],[],,[1+∈∈k k k x x b a ηη.
证明 因为在k I 上梯形公式的截断误差为],[),(12
13
+∈''-k k k k x x f h ηη，所
以
∑-=''-=-1
3
)](12[n k k n f h T I η.
)(x f ''在],[b a 上连续，由连续函数性质知，在],[b a 中存在点η，使
)()(11
ηηf f n n k k ''=''∑-=， 因此
∑-=''-=-1
03)](12[n k k n f h T I η∑-=''?-?-=10
2)(1)(12n k k f n a b h η)(122
ηf h a b ''?--=.证毕.
类似于复化梯形公式截断误差的推导，可得复化辛普森公式的截断误差
∑-=-
=-1
0)4(4)]()2(180[n k k n f h h S I η],[),()2
(180)4(4b a f h a b ∈??--=ηη. 复化的Simpson's 3/8 求积公式的截断误差
],[),()4
(945)(2)
6(6b a f h a b C I n ∈--
=-ηη.
【注】 复化梯形公式、复化辛普森公式和复化的Simpson's 3/8 公式都是有效的求积计算公式.
（1） 收敛性.由截断误差公式（可知，复化梯形公式、复化辛普森公式和复化的Simpson's 3/8 公式的误差阶分别为2h 、4h 、6h ，收敛性都是显然的.实际上，只要],[)(b a C x f ∈，则可得到收敛性，即
dx x f T b a
n n ?=∞
→)(lim ，dx x f S b a
n n ?=∞
→)(lim ，dx x f C b
a
n n ?=∞
→)(lim .
（2） 稳定性.由于n T 、n S 和n C 的求积系数均为正数，所以由定理2可知复化梯形公式、复化辛普森公式和复化的Simpson's 3/8 公式均是稳定的求积公式.
例2 计算积分?+=1
2
14
dx x
π. 解 1） 将积分区间[0，1]八等分，分点及分点处的函数值见表7-3，用复化梯形公式计算，得
)]}87()43()85()21()83()41()81([2)1()0({1618f f f f f f f f f T ++++++++=
≈π 138988.3≈.
2） 将积分区间[0，1]四等分，用复化辛普森公式计算，得
)85()83()81([4)]43()21()41([2)1()0({6414f f f f f f f f S +++++++?=
≈π 141593.3)]}8
7
(≈+f .
两种复化方法都用到表7-3中九个点上的函数值，它们的计算工作量基本上相同，但所得结果与积分真值 14159265.3=π相比较，复化辛普森所得近似值4S 远比复化梯形所得近似值8T 要精确.因此，在实际计算时，较多地应用复化辛普森公式.
例3 分别用复化梯形公式和复化辛普森公式计算?π
0sin xdx ，使误差不超过5102-?，问各取多少个节点？
解 由复化梯形公式的截断误差，有
5
02
2102sin 12)(12max -<
?
?
??≤''?--=-x n f h a b T I x n π
ππη， 所以，53
2
1024
?≥
πn ，360≥n ，因此，复化梯形公式至少应取361个节点.
由复化辛普森公式的截断误差，有
5
02
)4(4102sin 280)(180max -<
?
?
??≤?--=-x n f h a b S I x n π
ππη，
所以，55
4
105760
?≥
πn ，9≥n ，因此，复化辛普森公式至少应取10个节
点.
4） 变步长的梯形求积公式
（1） 一般复化求积公式的不足：
数值积分是积分近似计算的一类方法，在满足精度要求的前下， 人们希望计算工作量尽可能的小，而计算量大小的主要标志就是积分区间分割数n ，对于给定的精度要求，由于没有确定适当的区间分割数
n 的准则，所以，可能会出现：
▲ 区间分割数n 太小，达不到精度要求； ▲ 区间分割数n 太大，造成不必要的浪费。
(２) 改进的办法
以复化的梯形求积公式为例，将区间],[b a 按m n 2=分成小区间
],[1+k k x x ，1,,1,0-=n k ，则有m a
b h 2
-=
()()()()()??
?
???++≈==∑∑?
?-==+b f x f a f h dx x f dx x f T n k k k x x b
a
m
k k
m 1
120
2221
(30)
()()()()()∑∑∑∑?
?-=-=-==??
? ??++=??????++???
??++?≈==+++120212112020
222222222111
1m m m m m k k
m k k k k k k k x x b
a
h x f h T b f x f h x f a f h dx
x f dx x f T (31)
即有
∑-=+??
? ??+-+=+1201222221
m m m k k m h x f a b T T ， ,1,0=m (32)
例5 利用递推公式（32）重新计算积分?+=1
2
14
dx x π. 解 1） 首先对区间[0，1]使用梯形公式，得
3)24(2
1
)]1()0([2011=+=+-=
f f T . 2） 将区间二等分，新增分点2/1=x ，由递推公式（32）得
1.3)2
1
(212112=+=f T T .
3） 再将各小区间二等分，新增分点4/3,4/1=x ，由递推公式（32）得
13117647.3)]4
3
()41([412124=++=f f T T .
4） 将区间八等分，新增分点8/7,8/5,8/3,8/1=x ，由递推公式（32）得
13898849.3)]8
7
()85()83()81([812148=++++=f f f f T T .
这样不断二分下去，计算结果见表7-4，其中k 代表二分的次数，区间等分数k n 2=.
表7-4 梯形法递推计算值
表7-4说明用复化梯形公式计算积分，将区间二分9次，即有分点513个，达到7位有效数字，计算量很大.
§2.3 Romberg 求积公式
梯形求积公式包括复化的梯形求积公式由于是利用分段线性插值多项式代替被积函数做积分，所以是低阶的方法，一般来说，为获得较高的精度，区间分割数n 要取得很大，基于这样的事实，人们考虑这样一个问题：在不增加区间分割数n 的前提下，适当地使用低阶的方法提高计算精度，为此，我们考察复化的梯形求积公式的误差估计：
()()12
212ξf h a b T dx x f n b
a
''??
?
??-=-?
()()22
222112ξf h a b T dx x f n b
a
''??
?
???-=-?
由此得到下面的估计和结果：
()()()()()1
44422112212222
12
2--≈?≈''??
? ???-''??? ??-=--???n
n b a n b a n b
a
T T dx x f f h a b f h a b T dx x f T dx x f ξξ （33） ()()()()()()()()()n
n n k k n k k n n k k n n k k n k k n n S x f h x f x f x f h x f x f x f h x f h x f x f x f h T T =???? ??+??? ??
+++=??
? ??++-???
? ??+??? ??
+++=--∑∑∑∑∑-=-=-=-=-=1
11101
101
111022426223122244311
44 (34)
即有
1
442--=
n
n n T T S （35） ()()
()ξ44
2180f h a b S dx x f n b
a
??
? ??-=-?
（36）
将上述公式用于变步长的复化梯形求积公式，
1
442221--=+m m m
T T S （37）
按上述的思路，还有
()()()???≈--=?≈--++b a b a
b
a dx x f S S C S
dx x f S dx x f m m m
m m
1
4442222
22
2211
(38)
以及
()()
62h O C dx x f m b
a
≈-?
(39)
特别地
1
44322321--=
+m
m m C C R (40)
()()
82h O R dx x f m b
a
≈-?
(41)
上述的递推公式(37)、(38)和(40)统称为（Romberg Quadrature ）求积法，公式(40)称为Romberg 求积公式。
§2.4 Gauss 求积法
１）代数精度：衡量数值积分方法优劣的重要标准 定义１：如果由计算公式
()()k n
k k b
a
x f A dx x f I ∑?=≈=0
（41）
给出的数值积分方法对次数不超过n 次的多项式函数是精确的，而对次数大于n 次的多项式函数不一定精确，那么称该数值积分方法具有n 阶的代数精度。
显然，梯形求积公式具有１阶的代数精度；Simpson 求积公式具有２阶的代数精度。
按着代数精度的标准，自然要产生一个想法：
（１） 能否在不增加插值节点的前提下，尽可能地提高数值积分方法
的代数精度？
（２） 在插值节点的个数固定的前提下，代数精度最高能达到多高？
我们先来研究第二个问题，由于讨论的前提是插值节点的个数固定，所以，能够做的只能是通过对插值节点的选择来提高计算的代数精度。
设用n 次插值多项式代替被积函数，插值节点为k x ，n k ,,1,0 =，那么当数值积分公式具有m 阶的代数精度时，就应该使下列1+m 个等式成立：
()
1
10
11++=-+==?
∑i i b
a
n
k i
k k i
a b i
x A dx x ， m i ,,1,0 = （42） 由于在等式（42）中，有k x 和k A ，n k ,,1,0 =，总共22+n 个参数可共选择，或者说方程组（42）中总共有22+n 个未知数，唯一确定这22+n 个未知数需要且仅需要22+n 个方程，因此，代数精度m 最高可达到
12+n 阶。
定义２ 在积分区间[]b a ,上用n 次插值多项式代替被积函数且具有12+n 阶代数精度的数值积分方法称为Gauss 求积法。其插值节点称为Gauss 点。
现在我们来讨论第一个问题，
方程组（42）是非线性方程组，一般求精确解很难。
当积分区间是特殊的区间[]1,1-时，方程组（42）就成为下面的情形：
()i x A i n
k i
k
k +--=+=∑1111
， 12,,1,0+=n i （43） 已经证明，在区间[]1,1-上Gauss 点是1+n 次Legendre 多项式
()()(
)
[]
12
1
1111!121+++++-+=n n n n n x dx
d n x L （44） 的1+n 个零点。
据此，非线性方程组（43）的求解可分两步来完成：
１） 求出是1+n 次Legendre 多项式的1+n 个零点k x ，n k ,,1,0 =； ２） 将k x ，n k ,,1,0 =代入（43）后解线性方程组。
()()()()∏∏∑∏∏≠====≠==-=-?-=-n k
j j j k x
x n
j j n
k n k
j j j n j j x x x x dx d x x x x dx d k
00000 n k ,,2,1,0 = （９）
()()∏=-=n
j j x x x w 0
所以拉格朗日插值多项式
()()()()()∑
∑=='-==n
k k
k k n
k k k y x w x x x w y x L x L 00 （１０）
因此得到两点高斯-勒让德求积公式
)31()31()(1
1
f f dx x f +-
≈?
-.
同理可得三点高斯-勒让德求积公式
)515
(95)0(98)515(95)(1
1
f f f dx x f ++-≈
?
-.
对于一般区间[]b a ,上的积分，可先利用变换
()a b x a
b t ---=
21
（46）
用递推公式计算定积分 实验目的： 1．充分理解不稳定的计算方法会造成误差的积累，在计算过程中会导致误差的迅速增加，从而使结果产生较大的误差。 2．在选择数值计算公式来进行近似计算时，应学会选用那些在计算过程中不会导致误差迅速增长的计算公式。 3．理解不稳定的计算公式造成误差积累的来源及具体过程； 4．掌握简单的matlab语言进行数值计算的方法。 实验题目： 对n=0,1,2,…,20，计算定积分： 实验原理： 由于y(n)= = – 在计算时有两种迭代方法，如下： 方法一： y(n)=– 5*y(n-1),n=1,2,3, (20) 取y(0)= = ln6-ln5 ≈ 0.182322 方法二： 利用递推公式：y(n-1)=-*y(n),n=20,19, (1) 而且，由 = * ≤≤* =
可取：y(20)≈*()≈0.008730. 实验容： 对算法一,程序代码如下： function [y,n]=funa() syms k n t; t=0.182322; n=0; y=zeros(1,20); y(1)=t; for k=2:20 y(k)=1/k-5*y(k-1); n=n+1; end y(1:6) y(7:11) 对算法二，程序代码如下： %计算定积分； %n--表示迭代次数； %y用来存储结果； function [y,n]=f(); syms k y_20;
y=zeros(21,1); n=1; y_20=(1/105+1/126)/2; y(21)=y_20; for k=21:-1:2 y(k-1)=1/(5*(k-1))-y(k)/5; n=n+1; end 实验结果： 由于计算过程中，前11个数字太小，后9个数字比较大，造成前面几个数字只显示0.0000的现象，所以先输出前6个，再输出7—11个，这样就能全部显示出来了。 算法一结果： [y,n]=funa %先显示一y(1)—y(6) ans = 0.1823 -0.4116 2.3914 -11.7069 58.7346
电子一班王申江 实验九数值微积分与方程数值求解 一、实验目的 1、掌握求数值导数和数值积分的方法 2、掌握代数方程数值求解的方法 3、掌握常微分方程数值求解的方法 二、实验内容 1、求函数在指定点的数值导数。 () 23 2 123,1,2,3 026 x x x f x x x x x == >>syms x >>f=[x x^2 x^3;1 2*x 3*x^2;0 2 6*x]; >>F=det(f) F=2*x^3 >>h=0.1 >>x=[0:h:4]; >>f=2*x^3; >>[dy,dx]=diff_ctr(f,h,1); >>y1=dy(dx==1) y1=6.0000 >>y2=dy(dx==2)
y2=24.0000 >>y3=dy(dx==3) y3=54.0000 2、用数值方法求定积分。 （1） 210I π =?的近似值 a=inline('sqrt(cos(t.^2)+4*sin((2*t).^2)+1)'); I=quadl(a,0,2*pi) I = 6.7992 + 3.1526i （2）()1 202ln 11x I dx x +=+? b=inline('log(1+x)./(1+x.^2)'); I=quadl(b,0,1) I = 0.2722 3、分别用3种不同的数值方法解线性方程组。 6525494133422139211 x y z u x y z u x y z u x y u +-+=-??-+-=??++-=??-+=? A=[6,5,-2,5;9,-1,4,-1;3,4,2,-2;3,-9,0,2]; b=[-4,13,1,11]'; x=A\b
1 单选(2分) 利用MATLAB求解科学计算问题的优势是（）。 得分/总分 ? A. 算法最优 ? B. 不需要编写程序 ? C. 程序执行效率高 ? D. 编程效率高 正确答案：D你没选择任何选项 2 单选(2分) 在MATLAB命令行窗口输入命令时，可使用续行符，其写法是（）。 得分/总分 ? A. 省略号（…） ? B. 分号（;） ? C. 三个小数点（…） ? D. 百分号（%） 正确答案：C你没选择任何选项 3
下列语句执行后，D的值为（）。 1.A=[1:3;4:6]; 2.D=sub2ind(size(A),[1,1],[2,3]) 得分/总分 ? A. 3 6 ? B. 2 5 ? C. 3 5 ? D. 4 5 正确答案：C你没选择任何选项 4 单选(2分) ceil(-2.1)+floor(-2.1)+fix(-2.1)的结果为（）。 得分/总分 ? A. -7 ? B. -6 ? C. -5 ? D. -9 正确答案：A你没选择任何选项 5
下列语句执行后，x的值是（）。 1.log=1:5; 2.x=log(1) 得分/总分 ? A. ? B. 1 ? C. 数学常数e ? D. 报错 正确答案：B你没选择任何选项 6 单选(2分) 下列语句执行后，c的值是（）。 1.ch=['abcdef';'123456']; 2.c=char(ch(2,4)-1) 得分/总分 ? A. '4' ? B. 4 ? C. '3' ? D. 3
7 单选(2分) 产生和A同样大小的全0矩阵的函数是（）。 得分/总分 ? A. zero(size(A)) ? B. zeros(size(A)) ? C. size(zero(A)) ? D. size(zeros(A)) 正确答案：B你没选择任何选项 8 单选(2分) 语句x=speye(5)==eye(5)执行后，则下列说法中正确的是（）。 得分/总分 ? A. x是5阶全1矩阵，且采用稀疏存储方式 ? B. x是5阶全1矩阵，且采用完全存储方式 ? C. x是5阶单位矩阵，且采用稀疏存储方式 ? D. x是5阶单位矩阵，且采用完全存储方式
一、Maple V 系统 Maple V是由Waterloo大学开发的数学系统软件，它不但具有精确的数值处理功能，而且具有无以伦比的符号计算功能。Maple V的符号计算能力还是MathCAD和MATLAB等软件的符号处理的核心。Maple提供了2000余种数学函数，涉及范围包括：普通数学、高等数学、线性代数、数论、离散数学、图形学。它还提供了一套内置的编程语言，用户可以开发自己的应用程序，而且Maple自身的2000多种函数，基本上是用此语言开发的。 Maple采用字符行输入方式，输入时需要按照规定的格式输入，虽然与一般常见的数学格式不同，但灵活方便，也很容易理解。输出则可以选择字符方式和图形方式，产生的图形结果可以很方便地剪贴到Windows应用程序内。 二、MATLAB 系统 MATLAB原是矩阵实验室（Matrix Laboratory）在70年代用来提供Linpack和Eispac k软件包的接口程序，采用C语言编写。从80年代出现3.0的DOS版本，逐渐成为科技计算、视图交互系统和程序语言。MATLAB可以运行在十几个操作平台上，比较常见的有基于W indows 9X/NT、OS/2、Macintosh、Sun、Unix、Linux等平台的系统。 MATLAB程序主要由主程序和各种工具包组成，其中主程序包含数百个内部核心函数，工具包则包括复杂系统仿真、信号处理工具包、系统识别工具包、优化工具包、神经网络工具包、控制系统工具包、μ分析和综合工具包、样条工具包、符号数学工具包、图像处理工具包、统计工具包等。而且5.x版本还包含一套几十个的PDF文件，从MATLAB的使用入门到其他专题应用均有详细的介绍。 MATLAB是数值计算的先锋，它以矩阵作为基本数据单位，在应用线性代数、数理统计、自动控制、数字信号处理、动态系统仿真方面已经成为首选工具，同时也是科研工作人员和大学生、研究生进行科学研究的得力工具。MATLAB在输入方面也很方便，可以使用内部的E ditor或者其他任何字符处理器，同时它还可以与Word6.0/7.0结合在一起，在Word的页面里直接调用MATLAB的大部分功能，使Word具有特殊的计算能力。 三、MathCAD 系统 MathCAD是美国Mathsoft公司推出的一个交互式的数学系统软件。从早期的DOS下的1. 0和Windows下的4.0版本，到今日的8.0版本，功能也从简单的数值计算，直至引用Map le强大的符号计算能力，使得它发生了一个质的飞跃。 MathCAD是集文本编辑、数学计算、程序编辑和仿真于一体的软件。MathCAD7.0 Profe ssional（专业版）运行在Win9X/NT下，它的主要特点是输入格式与人们习惯的数学书写格式很近似，采用WYSWYG（所见所得）界面，特别适合一般无须进行复杂编程或要求比较特殊的计算。MathCAD 7.0 Professional 还带有一个程序编辑器，对于一般比较短小，或者要求计算速度比较低时，采用它也是可以的。这个程序编辑器的优点是语法特别简单。 MathCAD可以看作是一个功能强大的计算器，没有很复杂的规则；同时它也可以和Wor d、Lotus、WPS2000等字处理软件很好地配合使用，可以把它当作一个出色的全屏幕数学公式编辑器。 四、Mathematica 系统 Mathematica是由美国物理学家Stephen Wolfram领导的Wolfram Research开发的数学系统软件。它拥有强大的数值计算和符号计算能力，在这一方面与Maple类似，但它的符
一、符号积分 符号积分由函数int来实现。该函数的一般调用格式为： int(s)：没有指定积分变量和积分阶数时，系统按findsym函数指示的默认变量对被积函数或符号表达式s求不定积分； int(s,v)：以v为自变量，对被积函数或符号表达式s求不定积分； int(s,v,a,b)：求定积分运算。a,b分别表示定积分的下限和上限。该函数求被积函数在区间[a,b]上的定积分。a和b可以是两个具体的数，也可以是一个符号表达式，还可以是无穷(inf)。当函数f关于变量x在闭区间[a,b]上可积时，函数返回一个定积分结果。当a,b中有一个是inf时，函数返回一个广义积分。当a,b中有一个符号表达式时，函数返回一个符号函数。 例： 求函数x^2+y^2+z^2的三重积分。内积分上下限都是函数，对z积分下限是sqrt(x*y)，积分上限是x^2*y；对y积分下限是sqrt(x)，积分上限是x^2；对x的积分下限1，上限是2，求解如下： >>syms x y z %定义符号变量 >>F2=int(int(int(x^2+y^2+z^2,z,sqrt(x*y),x^2*y),y,sqrt(x),x^2),x,1,2) %注意定积分的书写格式 F2 = 1610027357/6563700-6072064/348075*2^(1/2)+14912/4641*2^(1/4)+64/225*2 ^(3/4) %给出有理数解 >>VF2=vpa(F2) %给出默认精度的数值解 VF2 = 224.92153573331143159790710032805 二、数值积分 1.数值积分基本原理 求解定积分的数值方法多种多样，如简单的梯形法、辛普生(Simpson)?法、牛顿－柯特斯(Newton-Cotes)法等都是经常采用的方法。它们的基本思想都是将整个积分区间[a,b]分成n个子区间[xi,xi+1]，i=1,2,…,n，其中x1=a，xn+1=b。这样求定积分问题就分解为求和问题。 2.数值积分的实现方法 基于变步长辛普生法，MATLAB给出了quad函数来求定积分。该函数的调用格式为： [I,n]=quad('fname',a,b,tol,trace) 基于变步长、牛顿－柯特斯(Newton-Cotes)法，MATLAB给出了quadl函数来求定积分。该函数的调用格式为： [I,n]=quadl('fname',a,b,tol,trace) 其中fname是被积函数名。a和b分别是定积分的下限和上限。tol用来控制积分精度，缺省时取tol=0.001。trace控制是否展现积分过程，若取非0则展现积分过程，取0则不展现，缺省时取trace=0。返回参数I即定积分值，n为被积函数的调用次数。 例： 求函数'exp(-x*x)的定积分，积分下限为0，积分上限为1。 >>fun=inline('exp(-x.*x)','x'); %用内联函数定义被积函数fname
东北大学秦皇岛分校 数值计算课程设计报告 数值积分算法及MATLAB实现 学院数学与统计学院 专业信息与计算科学 学号5133201 姓名陈悦 指导教师姜玉山张建波 成绩 教师评语： 指导教师签字： 2015年07月14日
1 绪论 数值分析是计算数学的一个主要部分，计算数学是数学科学的一个分支，它研究用计算机求解各种数学问题的数值检索方其理论与软件的实现.而数值分析主要研究数值计算. 现科学技术的发展与进步提出了越来越多的复杂的数值计算问题，这些问题的圆满解决已远人工手算所能胜任，必须依靠电子计算机快速准确的数据处理能力.这种用计算机处理数值问题的方法，成为科学计算.今天，科学计算的应用范围非常广泛，天气预报、工程设计、流体计算、经济规划和预测以及国防尖端的一些科研项目，如核武器的研制、导弹和火箭的发射等，始终是科学计算最为活跃的领域. 1.1 数值积分介绍 数值积分是数值分析的重要环节，实际问题当中常常需要计算积分，有些数值方法，如微分方程和积分方程的求解，也都和积分计算相联系. 求某函数的定积分时，在多数情况下，被积函数的原函数很难用初等函数表达出来，因此能够借助微积分学的牛顿-莱布尼兹公式计算定积分的机会是不多的.另外，许多实际问题中的被积函数往往是列表函数或其他形式的非连续函数，对这类函数的定积分，也不能用不定积分方法求解.由于以上原因，数值积分的理论与方法一直是计算数学研究的基本课题.对微积分学做出杰出贡献的数学大师，如I.牛顿、L.欧拉、C.F.高斯、拉格朗日等人都在数值积分这个领域作出了各自的贡献，并奠定了这个分支的理论基础. 构造数值积分公式最通常的方法是用积分区间上的n 次插值多项式代替被积函数，由此导出的求积公式称为插值型求积公式.特别在节点分布等距的情形称为牛顿-科特斯公式，例如梯形公式(Trapezoidal Approximations)与抛物线公式(Approximations Using Parabolas)就是最基本的近似公式.但它们的精度较差.龙贝格算法是在区间逐次分半过程中，对梯形公式的近似值进行加权平均获得准确程度较高的积分近似值的一种方法，它具有公式简练、计算结果准确、使用方便、稳定性好等优点，因此在等距情形宜采用龙贝格求积公式(Rhomberg Integration).当用不等距节点进行计算时，常用高斯型求积公式计算，它在节点数目相同情况下，准确程度较高，稳定性好，而且还可以计算无穷积分.数值积分还是微分方程数值解法的重要依据.许多重要公式都可以用数值积分方程导出.现探讨数值积分算法以及运用MATLAB软件的具体实现
第四讲绘图功能
作为一个功能强大的工具软件，Matlab 具有很强的图形处理功能，提供了大量的二维、三维图形函数。由于系统采用面向对象的技术和丰富的矩阵运算，所以在图形处理方面即常方便又高效。
4.1 二维图形 一、plot函数 函数格式：plot(x,y)其中x和y为坐标向量函数功能：以向量x、y为轴，绘制曲线。【例1】在区间0≤X≤2 内，绘制正弦曲线Y=SIN（X），其程序为： x=0:pi/100:2*pi; y=sin(x); plot(x,y)
一、plot函数 【例2】同时绘制正、余弦两条曲线Y1=SIN（X）和Y2=COS（X），其程序为： x=0:pi/100:2*pi; y1=sin(x); y2=cos(x); plot(x,y1,x,y2) plot函数还可以为plot(x,y1,x,y2，x,y3，…)形式，其功能是以公共向量x为X轴，分别以y1，y2，y3，…为Y轴，在同一幅图内绘制出多条曲线。
一、plot函数 （一）线型与颜色 格式：plot(x,y1,’cs’,...) 其中c表示颜色，s表示线型。 【例3】用不同线型和颜色重新绘制例4.2图形，其程序为：x=0:pi/100:2*pi; y1=sin(x); y2=cos(x); plot(x,y1,'go',x,y2,'b-.') 其中参数'go'和'b-.'表示图形的颜色和线型。g表示绿色，o表示图形线型为圆圈；b表示蓝色，-.表示图形线型为点划线。
一、plot函数 （二）图形标记 在绘制图形的同时，可以对图形加上一些说明，如图形名称、图形某一部分的含义、坐标说明等，将这些操作称为添加图形标记。 title(‘加图形标题'); xlabel('加X轴标记'); ylabel('加Y轴标记'); text(X,Y,'添加文本');
练习题 1.求函数在指定点的数值导数 x=sym('x'); >> y=[x x.^2 x.^3;1 2*x 3*x.^2;0 2 6*x]; >> x=1; >> eval(diff(y)) ans = 1 2 3 0 2 6 0 0 6 >> x=2; >> eval(diff(y)) ans = 1 4 12 0 2 12 0 0 6 >> x=3; >> eval(diff(y)) ans = 1 6 27 0 2 18 0 0 6 2.求下列函数导数 （1） x=sym('x'); >> y=x^10+10^x+(log(10))/log(x); >> diff(y) ans = 10*x^9+10^x*log(10)-2592480341699211/1125899906842624/log(x)^2/x （2） x=sym('x');
>> y=log(1+x); >> x=1; >> eval(diff(y,2)) %在x=1的条件下对y表达式求两次导数后导函数的值 ans = -0.2500 3.用数值方法求下列积分 首先先讲一下trapz的用法，如下题 t=0:0.001:1; >> y=t; >> trapz(t,y) ans = 0.5000 （1） >> x=1:0.01:5; >> y=(x.^2).*sqrt(2*x.^2+3); >> trapz(x,y) ans = 232.8066 (2) x=pi/4:0.01:pi/3; >> y=x./(sin(x).^2); >> trapz(x,y) ans = 0.3810 第三题拟合曲线题 x=[0:0.1:1]; >> y=[-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2]; >> a=polyfit(x,y,2); >> x=[0.05:0.2:1.05]; >> y=a(3)+a(2)*x+a(1)*x.^2 %注意x要在y前先赋值，不然y不会运行为最新的x对呀的y值 y =
《Matlab与科学计算》作业 第一章MATLAB环境 1、MATLAB通用操作界面窗口包括哪些？命令窗口、历史命令窗口、当前目录窗口、工作空间窗口各有哪些功能？ 答：MATLAB通用操作界面窗口包括：命令窗口、历史命令窗口、当前目录浏览器窗口、工作空间窗口、变量编辑器窗口、M文件编辑/调试器窗口、程序性能剖析窗口、MATLAB帮助。 命令窗口是MATLAB命令操作的最主要窗口，可以把命令窗口当做高级的“草稿纸”。在命令窗口中可以输入各种MATLAB的命令、函数和表达式，并显示除图形外的所有运算结果。 历史命令窗口用来记录并显示已经运行过的命令、函数和表达式，并允许用户对它们进行选择、复制和重运行，用户可以方便地输入和修改命令，选择多行命令以产生M文件。 当前目录窗口用来设置当前目录，可以随时显示当前目录下的M、MKL等文件的信息，扬文件类型、文件名、最后个修改时间和文件的说明信息等，并可以复制、编辑和运行M文件及装载MAT数据文件。 工作空间窗口用来显示所有MATLAB工作空间中的变量名、数据结构、类型、大小和字节数。 2、熟悉课本中表格1.4、1.5、1.6、1.7、1.8的内容。 3、如何生成数据文件？如何把数据文件中的相关内容输入到工作空间中，用实例进行操作。 生成数据文件：
把数据文件中的相关内容输入到工作空间中： 结果： 4、在工作空间中可以通过哪些命令管理变量，写出每种语法的具体操作过程。答：（1）把工作空间中的数据存放到MAT数据文件。 语法：save filename 变量1 变量2 ……参数。 （2）从数据文件中取出变量存放到工作空间。 语法：load filename 变量1 变量2 ……。
4.1数值微积分 4.1.1近似数值极限及导数 Matlab 数值计算中，没有求极限指令，也没有求导指令，而是利用差分指令： 用一个简单矩阵表现diff和gradient指令计算方式。 差分： Dx=diff(X) 对向量： Dx=X(2:n)-X(1:n-1) 对矩阵： DX=X(2:n,:)-X(1:n-1,:) 长度小1. DIFF(X), for a vector X, is [X(2)-X(1) X(3)-X(2) ... X(n)-X(n-1)]. DIFF(X), for a matrix X, is the matrix of row differences, （结果缺少一行） [X(2:n,:) - X(1:n-1,:)]. DIFF(X,N,DIM) is the Nth difference function along dimension DIM. If N >= size(X,DIM), DIFF returns an empty array （N阶差分）
梯度： FX=gradient(F) Fx(1)=Fx(2)-Fx(1); F=[1,2,3;4,5,6;7,8,9] Dx=diff(F) （按行） Dx_2=diff(F,1,2) （按列） [FX,FY]=gradient(F) Fx(1)=Fx(2)-Fx(1), Fx(end)=F(end)-F(end-1) FX与F维数相同。 [FX_2,FY_2]=gradient(F,0.5) %采样间隔0.5 即： Fx(1)=(Fx(2)-Fx(1))/2 F = 1 2 3 4 5 6 7 8 9 Dx = 3 3 3 3 3 3 Dx_2 = 1 1 1 1 1 1 FX = 1 1 1
实验一Matlab基本操作与微积分计算 实验目的 1.进一步理解导数概念及其几何意义. 2.学习matlab的求导命令与求导法. 3.通过本实验加深理解积分理论中分割、近似、求和、取极限的思想方法. 4.学习并掌握用matlab求不定积分、定积分、二重积分、曲线积分的方法. 5.学习matlab命令sum、symsum与int. 实验内容 一、变量 1、变量 MA TLAB中变量的命名规则是: （1）变量名必须是不含空格的单个词; （2）变量名区分大小写; （3）变量名最多不超过19个字符; （4）变量名必须以字母打头,之后可以是任意字母、数字或下划线,变量名中不允许使用标点符号. 1、创建简单的数组 x=[a b c d e f ]创建包含指定元素的行向量 x=first:step: last创建从first起,逐步加step计数,last结束的行向量, step缺省默认值为1 x=linspace(first,last,n）创建从first开始,到last结束,有n个元素的行向量 x=logspace(first,last,n）创建从first开始,到last结束,有n个元素的对数分隔行向量. 注:以空格或逗号分隔的元素指定的是不同列的元素,而以分号分隔的元素指定了不同行的元素. 2、数组元素的访问 （1）访问一个元素: x(i)表示访问数组x的第i个元素. （2）访问一块元素: x(a :b :c)表示访问数组x的从第a个元素开始,以步长为b到第c个元素（但
不超过c）,b可以为负数,b缺损时为1. （3）直接使用元素编址序号: x ([a b c d]) 表示提取数组x的第a、b、c、d个元素构成一个新的数组[x (a) x (b) x(c) x(d)]. 3、数组的运算 （1）标量-数组运算 数组对标量的加、减、乘、除、乘方是数组的每个元素对该标量施加相应的加、减、乘、除、乘方运算. 设:a=[a1,a2,…,an], c=标量, 则: a+c=[a1+c,a2+c,…,an+c] a .*c=[a1*c,a2*c,…,an*c] a ./c= [a1/c,a2/c,…,an/c](右除） a .\c= [c/a1,c/a2,…,c/an] (左除） a .^c= [a1^c,a2^c,…,an^c] c .^a= [c^a1,c^a2,…,c^an] （2）数组-数组运算 当两个数组有相同维数时,加、减、乘、除、幂运算可按元素对元素方式进行的,不同大小或维数的数组是不能进行运算的. 设:a=[a1,a2,…,an], b=[b1,b2,…,bn], 则: a +b= [a1+b1,a2+b2,…,an+bn] a .*b= [a1*b1,a2*b2,…,an*bn] a ./b= [a1/b1,a2/b2,…,an/bn] a .\b=[b1/a1,b2/a2,…,bn/an] a .^b=[a1^b1,a2^b2,…,an^bn] 三、矩阵 1、矩阵的建立 矩阵直接输入:从“[ ” 开始,元素之间用逗号“,”(或空格),行之间用分号“;”(或回车),用“ ]”结束. 特殊矩阵的建立: a=[ ] 产生一个空矩阵,当对一项操作无结果时,返回空矩阵,空矩阵的大小为零. b=zeros (m,n) 产生一个m行、n列的零矩阵 c=ones (m,n) 产生一个m行、n列的元素全为1的矩阵 d=eye (m,n) 产生一个m行、n列的单位矩阵 eye (n) %生成n维的单位向量 eye (size (A)) %生成与A同维的单位阵 2、矩阵中元素的操作 （1）矩阵A的第r行A（r,:） （2）矩阵A的第r列A（:,r） （3）依次提取矩阵A的每一列,将A拉伸为一个列向量A（:） （4）取矩阵A的第i1~i2行、第j1~j2列构成新矩阵:A(i1:i2, j1:j2) （5）以逆序提取矩阵A的第i1~i2行,构成新矩阵:A(i2:-1:i1,:） （6）以逆序提取矩阵A的第j1~j2列,构成新矩阵:A(:, j2:-1:j1 ） （7）删除A的第i1~i2行,构成新矩阵:A(i1:i2,:)=[ ] （8）删除A的第j1~j2列,构成新矩阵:A(:, j1:j2)=[ ] （9）将矩阵A和B拼接成新矩阵:[A B];[A;B] 3、矩阵的运算 （1）标量-矩阵运算同标量-数组运算. （2）矩阵-矩阵运算 a. 元素对元素的运算,同数组-数组运算.(A/B %A右除B; B\A%A左除B) b. 矩阵运算: 矩阵加法:A+B 矩阵乘法:A*B 方阵的行列式:det（A） 方阵的逆:inv（A）
用MATLAB 计算多元函数的积分 三重积分的计算最终是化成累次积分来完成的，因此只要能正确的得出各累次积分的积分限，便可在MA TLAB 中通过多次使用int 命令来求得计算结果。但三重积分的积分域Ω是一个三维空间区域，当其形状较复杂时，要确定各累次积分的积分限会遇到一定困难，此时，可以借助MATLAB 的三维绘图命令，先在屏幕上绘出Ω的三维立体图，然后执行命令 rotate3d on ↙ 便可拖动鼠标使Ω的图形在屏幕上作任意的三维旋转，并且可用下述命令将Ω的图形向三个坐标平面进行投影： view(0,0),向XOZ 平面投影； view(90,0),向YOZ 平面投影； view(0,90),向XOY 平面投影. 综合运用上述方法，一般应能正确得出各累次积分的积分限。 例11.6.1计算zdv Ω ???，其中Ω是由圆锥曲面222z x y =+与平面z=1围成的闭区域 解 首先用MA TLAB 来绘制Ω的三维图形，画圆锥曲面的命令可以是： syms x y z ↙ z=sqrt(x^2+y^2); ↙ ezsurf(z,[-1.5,1.5]) ↙ 画第二个曲面之前，为保持先画的图形不会被清除，需要执行命令 hold on ↙ 然后用下述命令就可以将平面z=1与圆锥面的图形画在一个图形窗口内： [x1,y1]=meshgrid(-1.5:1/4:1.5); ↙ z1=ones(size(x1)); ↙ surf(x1,y1,z1) ↙ 于是得到Ω的三维图形如图:
由该图很容易将原三重积分化成累次积分： 111zdv dy -Ω=???? 于是可用下述命令求解此三重积分： clear all ↙ syms x y z ↙ f=z; ↙ f1=int(f,z.,sqrt(x^2+ y^2),1); ↙ f2=int(f1,x,-sqrt(1- y^2), sqrt(1- y^2)); ↙ int(f2,y,-1,1) ↙ ans= 1/4*pi 计算结果为4 π 对于第一类曲线积分和第一类曲面积分，其计算都归结为求解特定形式的定积分和二重积分，因此可完全类似的使用int 命令进行计算，并可用diff 命令求解中间所需的各偏导数。 例11.6.2用MATLAB 求解教材例11.3.1 解 求解过程如下 syms a b t ↙ x=a*cos(t); ↙ y=a*sin(t); ↙ z=b*t; ↙ f=x^2 +y^2+z^2; ↙ xt=diff(x,t); ↙ yt=diff(y,t); ↙ zt=diff(z,t); ↙ int(f*sqrt(xt^2 +yt^2+zt^2),t,0,2*pi) ↙ ans= 2/3*( a^2 +b^2)^1/2*a^2*pi+8/3*( a^2 +b^2)^1/2*b^2*pi^3 对此结果可用factor 命令进行合并化简： factor （ans ） ans= 2/3*( a^2 +b^2)^1/2*pi*(3* a^2 +4*b^2*pi^2) 例11.6.3用MATLAB 求解教材例11.4.1 解 求解过程如下 syms x y z1 z2↙ f= x^2 +y^2; ↙ z1=sqrt(x^2 +y^2); ↙ z2=1; ↙ z1x=diff(z1,x); ↙ z1y=diff(z1,y); ↙ z2x=diff(z2,x); ↙ z2y=diff(z2,y); ↙
盐城师范学院《MATLAB与科学计算》期末论文 2016－2017学年度第一学期 用MATLAB解决解析几何的图形问题 学生姓名吴梦成 学院数学与统计学院 专业信息与计算科学 班级数15(5)信计 学号 15213542
用MATLAB 解决解析几何的图形问题 摘 要 将 MATLAB 的图形和动画功能都用于解析几何教学，可使教学形象生动。以图形问题为例，详细给出了实例的程序编写和动画实现过程 。在解析几何教学中有一定的应用价值。 【关键词】： MATLAB ； 解析几何 ；图形 ； 动 画；编程 1 引 言 在解析几何的教学中，使用传统的教学方法。许多曲线及曲面的形成过程与变换过程只通过传统的教师讲授静态图示就很难形象生动地表示出来 。在解析几何教学中使用MATLAB 软件辅助教学，不仅可以很容易绘制出复杂的立体图形。把曲线、曲面的形成和变化过程准确地模拟出来 ，而且还能够对它们进行翻转 、旋转 ，甚 至还能够轻而易举地实现图形的动画效果!这对提高教学效率和培养学生的空间想象能力可起到事半功倍的效果。下面结合实例从几个方面说明MATLAB 在解析几何画图方面的应用。 2 利用 MATLAB 绘制三维曲线 在空间解析几何中，各种曲线和曲面方程的建立都离不开图形 ，而空间曲线和曲面图形既难画又费时。借助MATLAB 的绘图功能 ，可以快捷 、 准确地绘出图形，使教学变得形象 、生动 。有利于学生观察三维空间图形的形状 ， 掌握图形的性质 。 一 般地 ，MATLAB 可用plot3，ezplot3，comet3等函数来各种三维曲线 。 例如画螺旋曲线的图形，其参数方程设为 ：t at cos x =，t b sin t y -=，ct =z 。使用 plot3语句画螺旋曲线图形的方法如下( 设a =2 ,b=4,c=3)： );*3),sin(*.*4),cos(*.*2(3;*10:50/:0t t t t t plot pi pi t -= MATLAB 用两条简单的语句就可以画出螺旋 曲线(图1)，但上述方法是静态的 ，为了体
M a t l a b数值积分与数值微分 Matlab数值积分 1.一重数值积分的实现方法 变步长辛普森法、高斯-克朗罗德法、梯形积分法 1.1变步长辛普森法 Matlab提供了quad函数和quadl函数用于实现变步长 辛普森法求数值积分.调用格式为: [I,n]=Quad(@fname,a,b,tol,trace) [I,n]=Quadl(@fname,a,b,tol,trace) Fname是函数文件名，a,b分别为积分下限、积分上限； tol为精度控制，默认为1.0×10-6，trace控制是否展 开积分过程，若为0则不展开，非0则展开，默认不展开. 返回值I为积分数值；n为调用函数的次数. --------------------------------------------------------------------- 例如：求 ∫e0.5x sin(x+π )dx 3π 的值. 先建立函数文件 fesin.m function f=fesin(x) f=exp(-0.5*x).*sin(x+(pi/6));再调用quad函数
[I,n]=quad(@fesin,0,3*pi,1e-10) I= 0.9008 n= 365 --------------------------------------------------------------------- 例如：分别用quad函数和quadl函数求积分 ∫e0.5x sin(x+π 6 )dx 3π 的近似值，比较函数调用的次数. 先建立函数文件 fesin.m function f=fesin(x) f=exp(-0.5*x).*sin(x+(pi/6)); formatlong [I,n]=quadl(@fesin,0,3*pi,1e-10) I= n= 198 [I,n]=quad(@fesin,0,3*pi,1e-10) I= n= 365 --------------------------------------------------------------------- 可以发现quadl函数调用原函数的次数比quad少，并 且比quad函数求得的数值解更精确. 1.2高斯-克朗罗德法
单元测验已完成成绩：分 1 【单选题】 MATLAB一词来自（）的缩写。 A、 Mathematica Laboratory B、 Matrix Laboratory C、 MathWorks Lab D、 Matrices Lab 我的答案：B得分：分 2 【单选题】 下列选项中能反应MATLAB特点的是（）。 A、 算法最优 B、 不需要写程序 C、 程序执行效率高 D、 编程效率高
我的答案：D得分：分 单元测验已完成成绩：分 1 【单选题】 当在命令行窗口执行命令时，如果不想立即在命令行窗口中输出结果，可以在命令后加上（）。 A、 冒号（:） B、 逗号（,） C、 分号（;） D、 百分号（%） 我的答案：C得分：分 2 【单选题】 fix(264/100)+mod(264,10)*10的值是（）。 A、 86 B、 62 C、 423 D、
42 我的答案：D得分：分 3 【单选题】 在命令行窗口输入下列命令后，x的值是（）。 >> clear >> x=i*j A、 不确定 B、 -1 C、 1 D、 i*j 我的答案：B得分：分 4 【单选题】 使用语句x=linspace(0,pi,6)生成的是（）个元素的向量。 A、 8 B、 7 C、 6
D、 5 我的答案：C得分：分 5 【单选题】 ceil的结果为（）。 A、 -2 B、 -3 C、 1 D、 2 我的答案：A得分：分 6 【单选题】 eval('sqrt(4)+2')的值是（）。 A、 sqrt(4)+2 B、 4 C、 2 D、
2+2 我的答案：B得分：分 7 【单选题】 已知a为3×5矩阵，则执行完a(:,[2,4])=[]后（）。 A、 a变成行向量 B、 a变为3行2列 C、 a变为3行3列 D、 a变为2行3列 我的答案：C得分：分 8 【单选题】 在命令行窗口输入以下命令 >> A=[1:3;4:6]; >> D=sub2ind(size(A),[1,1],[2,3]) D的值为（）。 A、 3 6 B、 2 5 C、 4 5
matlab 微积分基本运算 §1 解方程和方程组解 1. 线性方程组求解 对于方程 AX = B ,其中 A 是( m ×n )的矩阵有三种情形： 1)当n=m 且A 非奇异时,此方程为“恰定”方程组。 2)当 n > m 时,此方程为“超定”方程组。 3)当n 0.3188 两种方法所求方程组的解相同。 (2)MATLAB 解超定方程AX=B 的方法 对于方程 AX = B ,其中 A 是( m ×n )的矩阵， n > m ，如果A 列满秩，则此方程是没有精确解的。然而在实际工程应用中，求得其最小二乘解也是有意义的。基本解法有： 1)采用求伪逆运算解方程 x=pinv(A)*B 说明：此解为最小二乘解x=inv(A ’*A)*A*B,这里pinv(A) =inv(A ’*A)*A. 2)采用左除运算解方程 x=A\B 例2 “求伪逆”法和“左除”法求下列方程组的解 ??? ??=+=+=+1 221421 221 2121x x x x x x 命令如下： >> a=[1 2;2 4;2 2]; >> b=[1,1,1]'; >> xc=a\b %用左除运算解方程 运行得结果: xc = 0.4000 0.1000 >> xd=pinv(a)*b %用求伪逆运算解方程 运行得结果: xd = 0.4000 0.1000 >> a*xc-b %xc 是否满足方程ax=b 运行得结果: ans = -0.4000 0.2000 0.0000 可见xc 并不是方程的精确解。 (3) MATLAB 解欠定方程AX=B 的方法 欠定方程从理论上说是有无穷多个解的，如果利用求“伪逆”法和“左除”法来求解，只能得到其中一个解。基本方法： 1)采用求伪逆运算解方程 x=pinv(A)*B 2)采用左除运算解方程 Matlab 与科学计算考试样题（客观题） 1 下面的MATLAB 语句中正确的有： a) 2a ＝pi 。 b) record_1=3+4i c) a=2.0, d) c=1+6j 2. 已知水的黏度随温度的变化公式如下，其中a=0.03368，b=0.000221，计算温度t 为20，30，40度时的粘度分别是： 2 1at bt μμ=++0μ为0℃水的黏度，值为31.78510-?；a 、b 为常数，分别为0.03368、0.000221。 3. 请补充语句以画出如图所示的图形： [x,y]=meshgrid(-2:0.1:2, -2:0.1:2)。 Z=x.*exp(-x.^2-y.^2)。 。 a) Plot3(x,y,Z) b) plot3(x,y,Z) c) mesh(x,y,Z) d) plot3(x,y,z) 2 a) 0.4900 1.2501 0.8560 b) 0.8560 1.2501 0.4900 c) -0.6341 3.8189 -3.7749 d) 3.8189 -3.7749 2.8533 解释说明： >> x=0.5:0.5:3.0。 >> y=[1.75,2.45,3.81,4.80,8.00,8.60]。 >> a=polyfit(x,y,2) a = 0.4900 1.2501 0.8560 >> x1=[0.5:0.25:3.0]。 >> y1=a(1)*x1.^2+a(2)*x1+a(3) >> plot(x,y,'*') >> hold on >> plot(x1,y1,'--r') 5. 求方程在 x=0.5附近的根. 21 x x += a) 0.6180 b) -1.1719e-25 c) －1 d) -1.6180 6. 用Newton-Cotes方法计算如下积分 1 5 x? （a）133.6625 (b)23.8600 (c) 87.9027 (d) -1.6180 7. y=ln(1+x)，求x=1时y" a) -0.25 b) 0.5 c) -0.6137 d) -1.6137 8.某公司用3台轧机来生产规格相同的铝合金薄板。取样测量薄板的 厚度，精确至‰厘M。得结果如下： 轧机1：0.236 0.238 0.248 0.245 0.243 轧机2：0.257 0.253 0.255 0.254 0.261 轧机3：0.258 0.264 0.259 0.267 0.262 计算方差分析结果，并判定各台轧机所生产的薄板的厚度有无显著的差异？ a) p＝1.3431e-005，没有显著差异。Matlab与科学计算样题(加主观题答案)
相关主题
文本预览