Transcript 概率论预备知识

第七章
计算机代数
7.1 引言
• 数值计算系统
Fortran语言，…
例： …X=1 Y=2 Z=X+Y…
• 非数值计算系统
Mathematica语言，…
例：勒让德多项式定义
Pn ( x )
递推公式：
P0 ( x )  1, P1 ( x )  x , Pn ( x )  [(2 n  1) xPn 1 ( x )  ( n  1) Pn  2 ( x )] / n
1
计算机代数系统的发展历史
二十世纪六十年代最早的计算机代数系统几乎完全是基于
LISP表处理语言.它是用来处理表链的。它对于早期符号计算
程序的重要性，就好比同一时期处理数值计算的程序FORTRAN
系统。在这个阶段，REDUCE 程序对高能物理已经表现出一些
特殊的用途。
 SCHOONSCHIP 是M. Veltman 用汇编语言写的，专门应用于
粒子物理领域。汇编代码的应用导致了难以置信的高速计算程
序（相对于最初的解释代码），从而使计算更复杂的高能物理
散射过程成为可能。由于人们逐渐认识到这个程序的重要性，
因而，1998 年M.Veltman 因此获得了诺贝尔物理奖。

2
同时值得一提的是基于Franz LISP 的MACSYMA 系统，
它引发了算法的重要发展。
 从1980 年以来,新的计算机代数系统开始采用C 语言
编写。这个时期还出现了最早的商业计算机代数系统，
其中Mathematica 和Maple 最为著名。
 近几年，有关大型程序可维护性的问题变得越来越重
要。全部的设计范例都由过程设计变到了面向对象设计。
反映在编程语言上从C 变到C++。这样GiNaC 库随之发展
起来，它支持C++环境下的符号计算。

3
7.2 Mathematica 软件运算和编程
Mathematica: 该系统是美国Wolfram 公司开发的一个
功能强大的计算机通用数学系统。其基本系统主要是
用C 语言开发的。它是当前运用十分广泛的符号代数
处理系统。
Mathematica系统可以在交互式的状态下运行。作为
一种计算机语言，它对于各种数、变量、函数、代数表达
式和语句都有比较严格的要求和规定。本节将仅介绍利用
Mathematica语言进行运算和编程的基本知识。
4
Mathematica概述
5
6
7
Mathematica的基本量
8
I
I
I
9
I
10
11
12
13
例如：按Mathematica的格式输出下列表达式
(1) ( 56 + 38) ´ (39 - 23) ¸ 29
2
(2) 67 ´
(3) (
97
39 ¸ 9
2
+
1372
(4)
3
2
a
57 ) ¸ [
37
19
+ b ( x - 9)
11
6
´ (
12
3
) ]
7
2
2 + 3x
14
15
16
计算下列各式：
(1)
(2)
3
2009  2006  2010
2
4  7 i
(  4 i )(2  5 i )
(3) a  2; b  11; ln( a )  ( L og 3 b )
3
(4) x 
2
6
; y 
5
2
 ; sin ( x )  co s ( y )
2
2
6
(5) x = 3; a = 1; b = 2;
a
+ b ( x - 9)
2
2 + 3x
17
^
18
19
20
定义以下函数并计算函数值：
2
(1) f ( x ) 
3
sin ( 2 x )
并计算： f ( x 
x 1
)
6
(2) f ( x )  e cos(2 x )
x
(3) f ( x , y )  ln( x

5 y
)
并计算： f ( x  3.22)
并计算：f ( x  2, y= 3)
ìï
0,
x< 0
ï
ï 2
ï x / 2,
0£ x< 1
ï
并计算： f ( x  0.5)
(4) f ( x ) = í
2
ï 2 x - x / 2 - 1, 1 £ x < 2
ï
ï
ï
1,
x³ 2
ïî
21
stop
22
如： In [1] := T ab le[i, {i, 1, 5}]
O u t[1] := {1, 2, 3, 4, 5}
In [2] := T ab le[i, {i, 0, 0.5, 0.1}]
O u t[2] := {0, 0.1, 0.2, 0.3, 0.4, 0.5}
In [3] := v = T ab le[2i + 1, {i, 1, 5}]
O u t[3] := {3, 5, 7, 9, 11}
In [4] := v[[3]]
O u t[4] := 7
23
1.用Table命令生成如下的表：
(3, 6, 9,12,15)
并计算其中第二个元素的2倍后再加上第三个元素的结果。
2.用Table命令生成如下的表：
(1, 3), (2, 5), (3, 7), (4, 9), (5,11)
3.矩阵运算
24
(1)求矩阵乘积
(2)求矩阵逆
(3)求矩阵行列式
25
注意：两个逻辑运算符： & &
表示and,两者都； ||
表示or,或者。
26
Mathematica的基本运算
27
^
28
（1）将下式分解质因式：
2  x  2 y  xy  y  y
2
3
（2）将下式展开：
( x  2)(6 x  3 x  5)
2
( x  y  z)
3
（3）将下式在x=1附近作泰勒级数展开：
1
1 x
29
30
31
32
33
求解下列方程：
(1) 3 x  9 x  x  7  0
3
(2) e
x
2
 2x  x  0
3
1  x 2  y 2  0
(3) 
2
2
2  x  y  0
34
Mathematica的函数作图
35
36
37
38
39
绘制下列函数图：
(1) f ( x )  x  6 x  3 x  26 x  24, x  [  3, 5]
4
3
(2) f ( x )  e cos(2 x ),
x
2
x  [0, 4]
(3) f 1 ( x )  sin x , f 2 ( x )  sin x
1.5
,
实线
虚线
x  [0, 4]
画在一幅图中，f1(x)用红色，f2(x)用蓝色，
在横轴上标识“time”，在纵轴上标识“position”。
40
41
42
43
44
1.将以下各点连线作图：
(1, 3), (2, 5), (3, 7), (4, 9), (5,11)
2.作以下参数方程的图线：
 x  3 sin t

 y  cos t
2.作以下三维图形：
f ( x , y )  sin( x  y )
2
2
45
微积分的基本操作
46
47
48
49
1.计算极限： lim sin x
x
2x 1
2.计算一阶和二阶微分：
y x
4/5
 sin
x
3.计算以下积分：

1
1
2 1  x dx
2
50
微分方程的求解
未知待求的
51
52
求解下列方程：
(1) y   y ,
2
y (0)  1
(2) y    y  x  1,
y (0)  1, y (0)  0
53
Mathematica程序设计
1．控制选择 （i） 顺序控制 分号“;”
（ii） 条件控制
If语言结构。
If[逻辑表达式，表达式]
If[逻辑表达式，表达式1，表达式2]
If[逻辑表达式，表达式1，表达式2，表达式3]
例如：
f[x_] := If [(x > 0) || (x=0) , N[Sqrt[x]], Print[“x is negative。”],
Print[“x is not numerical.”]]
54
Which语句结构。
Which[条件1，表达式1，条件2，表达式2，…, 条件n，表达式n]
例： Which[ 2= =3,x,3= =3,y ]
其结果为y。这是由于条件2==3的结果为False,
而条件3==3的结果为True.
55
Switch语句结构。
Switch[判别表达式，模式1，表达式1，模式2，表达
式2，…, 模式n，表达式n]
例： i=1
Switch[i^2,0,x,1,y,2,z]
最后结果为y。
56
（iii）循环控制
Do语句结构。
Do[表达式，{循环描述}]
例： Do[ Print[2^i],{i,1,5} ]
该指令的结果是循环打印出，
（i=1,2,3,4,5）的值2，4，8，16，32。
57
For语句结构。
For[初始表达式，条件，步进表达式，循环表达式]
例如： For[ i=0,i<=10,i++,Print[i] ]
值为0，1，2，3，4，5，6，7，8，9，10。
While语句结构。
While[条件，循环表达式]
例如： i=0
While[ i<=10,Print[i];i++ ]
58
2. 子程序
一般形式为：
Module[{<局部变量名表>},表达式1；表达式2；…表达式n]
例如：
unit[x_, y_, z_] := Module[{len}, len = Sqrt[x^2 + y^2 + z^2];
N[{x/len, y/len, z/len }] ]
59
3.程序包结构
一般的程序包具有如下的基本框架：
BeginPackage[“ 程序包名称 ”]
名字::usage=“字符串，程序包中定义在包外可以使用
的函数、变量等的名字和使用说明。”
……
Begin[“ Private ”]
程序包主体.(包括外部可用的函数及变量,一
些内部函数变量的定义.)
……
End[ ]
EndPackage[ ]
60
编程练习
y
1.用投点法求圆周率：
2r
方法：随机、均匀地在正方形
范围内投点，若共投了N个点，
而落在圆内的点数为M，当N足
够大时，有：
Sr
Ss
»
M
其中： S r   r 2 ,
N
所以： p » 4
r
o
r
S s  4r
2r
x
2
M
N
61
m = 0,
n = 50000
i= 0
x = 2 r h1 ,
y = 2rh2
N
2
2
(x - r) + ( y - r) < r
Y
m = m+1
i < n, i + +
Y
N
pi = 4 m / n
62
2.用直接抽样法对满足以下分布的随机变量抽样：
ìï 2 x , 0 £ x £ 1,
f ( x ) = ïí
ïïî 0, 其 它
解：用直接抽样法，首先求出分布函数
F ( x) = x
2
则满足分布的随机变量抽样为：
x=
x
其中 x 是[0,1]区间均匀分布的随机数。
63
1
L
f ( 1 )
3.用第一类舍选法对满足以下分布的随机
变量抽样：
ì
ïï 2 x , 0 £ x £ 1,
f ( x) = í
ïïî 0, 其 它
解：
64
4.用Metropolis方法对正态分布密度函数 f ( x ) = e
进行抽样，其中 - 10 £ x £ 10.
解：引入过渡概率： ( x  x )  m in {1,
x[ 0 ,4 ]
Metropolis计算步骤：
f ( x )
- x
2
}
f ( x)
1）选择初始位置： 0  x 0  0, f max ( x  0)  1
2）取 x1 Î [0,1], L = 10, 定义步长：di = - L + L x1 , di Î [- 5, 5]
3）取  2  [0,1], 判断： 2   ( xi  xi 1 ),
2
4）如果是，则 xi + 1 = xi + di , h k = xi + 1 然后返回步骤1
5）如果否，则 h k = x i , 返回步骤2
65
x1 , x 2
d n = - L / 2 + L x1 ,
x try = x n + dn
x 2 £ w ( x n ® x try )
N
Y
x n + 1 = x try
66
编程练习题
1.用投点法求圆周率：
2.用直接抽样法对满足以下分布的随机变量抽样：
ìï 2 x , 0 £ x £ 1,
f ( x ) = ïí
ïïî 0, 其 它
3.用第一类舍选法对满足以下分布的随机变量抽样：
ìï 2 x , 0 £ x £ 1,
f ( x ) = ïí
ïïî 0, 其 它
4.用Metropolis方法对正态分布密度函数 f ( x ) = e
进行抽样，其中 - 10 £ x £ 10.
- x
2
67
练习： 用改进Euler方法解方程
ìï y ' = y 2
ï
í
ï y (0 ) = 1
ïî
( x Î [0, 0 .4 ])
(0)
y i+ 1 = y i + h y i
(1)
y i+ 1 = y i +
h
2
2
2
(0)
2
( y i + ( y i+ 1 ) )
注意：取步长h=0.1，计算 y1 , y 2 , y 3 , y 4 。
68
1.《数学运算大师 Mathematica 4》, 洪维恩，人民邮
电出版社。
2.《Visual Fortran 常用数值算法集》, 何光渝，高永
利，科学出版社。
69