数值计算程序

附录一用Mathcad进行数值计算
Mathcad是世界流行的一个优秀的交互式应用数学软件。

其初版由美国Mathsoft公司于1986年推出，1997年问世的7.0版本可在Windows95/98或Windows NT下运行。

Mathcad是集数理计算、图形和文字处理等功能于一体的多媒体综合应用软件，它支持OLE2 技术，可以和其它应用软件协同工作，还具有网络浏览功能。

同时，它又是一个大众化的数学工具软件，用户只须学会简单的操作，无须记住很多的计算机知识和数学技巧，就可利用它进行各种复杂的数学计算、数据处理、作图和文本编辑。

Mathcad因功能强大且使用简便而享誉全世界，它拥有约100多万用户。

无论对于大学师生还是对于企业界的技术人员和管理者,Mathcad都是解决应用数学问题的好帮手。

下面主要针对前几章所涉及的数学模型问题简要介绍如何使用Mathcad7.0进行常用的科学计算。

1.Mathcad 基本用法
Mathcad 7.0分标准版，学生版和专业版，其中专业版(Mathcad7 Professional)的功能最齐全。

1.1 工作页和工具栏
进入Windows95/98后,在”开始/程序”的菜单中单击”Mathcad 7 Professional”项,
即可启动该系统.屏幕显示启动画面后进入用户界面,这时屏幕中央有一个”Tip of the Day”窗口,单击ok后,屏幕显示如图1
145
图1 Mathcad用户界面
在系统标题行下的主菜单中含8个主菜项，即File(文件管理)、Edit(编辑)、View(视图)、Insert(插入)、Format(格式)、Math(数学计算)、Symbolics(符号计算)、Window(窗口管理)、Help(帮助)。

其中大部分菜单功能与WORD97中同名菜单项类似。

单击Help可看到其中有若干项帮助命令，常用的是在线帮助(Mathcad help)以及资源中心(Resource Center)。

当你打开资源中心后可选择进入浏览(Overview)、教程(Tutorial)、快速参考活页(QuickSheets)等多项帮助内容,获取有针对性的指导及示例。

字体、常用工具栏的大部分功能也同WORD类似。

常用工具栏中的f(x)按钮可用于选择系统的内建函数输入到表达式中去。

Mathcad7有200多个内建函数可供我们使用，如果记住了函数名及其变量表的格式，也可以直接用键盘输入所要的函数。

数学工具栏见图2，它含8个按钮，从左到右依次为：
(1)计算器(2)逻辑运算(3)函数作图(4)矩阵和向量运算
(5)微积分运算(6)编程工具(7)希腊字母(8)符号运算
用每个按钮可打开一块工具板，从中可选择所需要的数学符号、运算符或作图命令。

Mathcad工作页的排版是二维的，它用起来
象一张向右向下都无限伸展的白纸，只是用边界
线(虚线)分割不同的页面。

往工作页上键入式子
时，一个式子常可以跨页书写，但若无必要则应图2 数学工具栏
避免式子跨页。

因为打印时，界线两侧分属不同页面，跨页的式子便被分割开了。

光标在工作页的空白处呈红十字(见图1)，可以用箭头键来移动光标。

整个工作页是一个Mathcad文件，当需要保存它时，用主菜单“文件/保存”命令即可，默认的文件名后缀为.mcd。

下次使用该文件时用主菜单“文件/打开”命令便可打开它。

1.2 数学区、文本区、作图区
在工作页中可建立数学区(也称计算区)、文本区、作图区。

在工作页的任何空白处键入任一字符都被默认为开辟一个数学区(也可用Insert/Math Region菜单命令来开辟)，并用矩形框围住，红十字光标变成直角状兰色编辑线将字符半包围，如图3。

当光标移出该区后又变成红十字，数学区域框隐去。

在数学区可以方便地键入表达式并进行计算，也可键入汉字，要注意：(1)每个数学区只能有一个独立的表达式；(2)一个文件中所有数学区的字体格式及字号均是一致的，需要改变时，只能将光标放入任一数学区内通过主菜单Format/Equation进行统一修改；(3)键入空格的效果是：或者移动并扩大编辑线包围的范围，或者将所在数学区变成文本区。

对数学区的计算命令有自动和手动两种模式，用主菜单Math/Automatic Calculation命令转换两种模式。

通常默认为自动计算模式，在信息栏可见到”Auto”信息。

在Auto模式下工作页是”活”的，成为一个广义的“Mathcad程序”，只要修改一处式子，其后与之有计算关系的
146
147 表达式的计算结果以及相关图形也自动随之改变。

在手动模式下，按F9可下达计算命令。

有时需要将某个数学区“关闭”，即禁止计算，方法是光标放入该区，用菜单Format/Properties/Calculation 命令，选择Disable Evaluation 并单击OK 。

被关闭的数学区右上角出现■记号。

要把被关闭的数学区打开的方法类似，取消Disable Evaluation 选择即可。

文本区可在工作页上任意空白处开辟。

方法
有三种：
(1)键入双引号”。

(2)单击主菜单Insert/Text Region 。

(3)键入任一字符后按空格键。

光标在文本区内呈红色竖线，即插入线，这
时可见到文本区的边界，是一个带有三个“把手”的矩形框，如图3所示。

文本区一般是不被计算的，但需要时可在文本区内插入数学区。

方法是让红色插入线位于插入处，点主菜单Insert/Math Region 。

插入文本区的数学区可以被计算，除非把它关闭。

作图区的开辟可用Insert/Graph 菜单命令，再选择坐标系。

简单的办法是利用数学工具栏上的作图工具板进行直观的选择。

也可用快捷键，如键入@便可开辟一个直角坐标作图区。

不论是数学区、文本区还是作图区，都可以用鼠标左键按住其边框拖拽移动到任何地方（这时鼠标呈手掌形）。

如果想同时拖拽若干个区，只须先同时选中它们即可，方法是用左键按住待选区域附近的空白处，拖出一个较大的虚框，松开左键，该虚框所及之处的所有区域分别被小虚框围住，这就表示被选中了。

被选中的区域除能被移动外，还可被复制、删除等，单击右键可弹出文中菜单进行这些处理。

当作了某些删除或移动后，可以用Ctrl+R 键来刷新屏幕。

若需将某些区域横向对齐或竖向对齐，只须选中这些区，再单击常用工具栏上f(x)左侧的两个按钮之一。

1.1 在数学区键入和计算表达式
利用Mathcad 提供的快捷键和大量的数学工具板，可以非常方便地输入数学表达式，并且与通常书写的形式差不多。

例如按从上到下的次序键入下面左列的字符串（ 表示键入空格），每个式子键入完毕后将光标移出所在区域或按回车键结束。

显示结果如右列所示。

键 入 显 示（设系统在默认Auto 模式下）
x^-5 x -5
x| |x|
x:5 x:=5
sin(x)= sin(x)=-0.959
x\= x =2.236
3*x= x 3⋅=15
3+x /5=
5
x 3+=1.6 3+x/5= 5
x 3+=4 y[0:1 y 0:=1 图3 数学区和文本区
148
x!= x!=120
要注意“:=”是赋值号，只须键入“：”即可产生；单独的“=”是表示计算命令的等号，其右端为机器自动给出的计算结果。

上面第三行将x 赋值5，所以后面凡与x 有关的算式均自动得出计算结果。

你可将x:5改为x:6，看结果会如何变化？此外，要注意空格键的作用。

上面的y 0:=1表示已定义y 为向量且第0个分量赋值1，Mathcad 中向量及矩阵的起始下标为0，必要时可用主菜单Math/Options 命令，在Array 栏中选其它整数作为起始下标。

数据的显示长度默认为到3位小数，若需显示更多位，只须光标位于某数学区内时用主菜项Format/Number 命令，在Precision 栏中修改Displayed Precision 的值即可。

一些特殊符号可利用数学工具栏的工具板，特别是微积分符号必须用工具板输入才能在计算时被系统认识。

打开数学工具栏中微积分运算工具板，当光标在数学区或工作页空白处时，单击n 、dx d 、⎰b a 、∑=m n 1
按钮之一，便分别显示如图4的算符： 其中小黑块是占位符，在占位符中分别键入适
当字符便可完成式子的输入，当占位符空缺而光标
移出该区时占位符呈红色，表示该式尚未写完。

在Auto 模式下，系统自动从上到下逐行扫描，
每行又从左到右、不受页边限制地逐个扫描数学区
并依次执行计算。

所以键入计算公式时要注意先后次序。

下面各节举例说明如何用Mathcad 解决常用数学模型的计算问题。

2. 求解一元方程
求根函数root(f(x),x)可用于求解一元方程f(x)=0.调用root 函数前须先定义f(x)，并对x 赋以初始近似值（如果求复根则初始近似值应为复数）。

此外还要给系统常量TOL 指定迭代的精度，否则默认TOL 的值为10-3。

例1 求方程x －cos x=0在0.5附近的根，要求误差不超过10-6。

解 令 f(x):=x-cos(x) x:=0.5 TOL:=10-6
用主菜单Format/Number 命令在Precision 栏修改数据显示为8位小数，再调用求根函数得所求的根为
root(f(x),x)=0.73908513
上例的求解过程中，“:=”右边是人工键入的，而单独的“=”右边是计算机给出的结果。

以下各例也相同，不再一一指出。

例2 利用图形估计方程x －cos x =0的初始近似根。

解 设 f(x):=x-cos(x) 再作以下步骤：
图4 特殊数学符号
149
(1)在上式下方空白处键入@，出现作图
区，内有矩形坐标框。

(2)在坐标框下面占位符处键入x ，左面占
位符处键入f(x)，鼠标左键单击作图区外任一
点，函数图形便自动生成。

(3)双击该图弹出对话框，在菜单的Axes Style 栏中选Crossed ，单击确定
按钮。

函数图形便如图5所示。

(4)观察曲线与x 轴的交点，估计初始近似
根可取为0.5。

(5)若需更精确些，单击
Format/Graph/Trace ，再单击曲线与x 轴的交点，出现过此点的十字交叉虚线，交点的坐标在弹出的对话框内，其中x 的值可取为初始近似根。

3. 线性代数计算
出对话框，填入行数、列数，单击确定，然后在占位符中键入具体字符即可。

用函数lsolve(A,b)可实现列主元高斯消去法求解方程组Ax=b ，结果输出所求的解向量。

例3 求解例3.2的线性方程组（设未知量改为x 0,x 1,x 2）。

解 令
A:=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡----5.04525.015.0201.0 b:= ⎥⎥⎥⎦
⎤⎢⎢⎢⎣⎡-955 x:=lsolve(A,b) 得 x T =[0 -2 2] 即 x 0=0 x 1=-2 x 2=2
（x T 中的T 只能用矩阵运算工具板的M T 键入，不能用^键入，否则系统不能正确计算。

） 例4 求例3.4中矩阵A 的行列式及逆阵。

解 A:= ⎥⎥⎥⎦
⎤⎢⎢⎢⎣⎡643542321 A -1=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡----0125.05.15.1102 |A|=-2 例5 对例3.7的系数矩阵进行选主元的LU 分解。

解 令 A:= ⎥⎥⎥⎦
⎤
⎢⎢⎢⎣⎡---196314112 M:=lu(A) (调用内建的lu 函数) 图5 函数f(x)的图形
150
得 M=⎥⎥⎥⎦
⎤⎢⎢⎢⎣⎡---714.1001143.05.00015.55.100015.1100314001010 （M 为3个3阶方阵并列组成） 令 P:=submatrix(M,0,2,0,2) L:= submatrix(M,0,2,3,5) Q:= submatrix(M,0,2,6,8)
则 P=⎥⎥⎥⎦
⎤⎢⎢⎢⎣⎡001100010 L=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡1143.05.0015.1001 U=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---714.1005.55.100314 满足 A P ⋅=U L ⋅ （这个等号是逻辑等号，用数学工具栏的逻辑运算工具板键入。

） 可分别计算A P ⋅和U L ⋅，验算上述等式。

上例中内建函数submatrix(M,a,b,c,d)是从矩阵M 中取行下标a~b ，列下标c~d 的子矩阵。

例6求例3.9的系数矩阵的LL T 分解。

解 令A:=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡1102021211 L:=cholesky(A) 得 L=⎥⎥⎥⎦
⎤⎢⎢⎢⎣⎡-732.122011001 例7求例4.4的矩阵A 的QR 分解。

解 令 A:= ⎥⎥⎥⎦
⎤
⎢⎢⎢⎣⎡--142222241 M:=qr(A) 得 M=⎥⎥⎥⎦
⎤
⎢⎢⎢⎣⎡------300333.0667.0667.0060667.0333.0667.0003667.0667.0333.0 (M 由Q 和R 并列构成) 令 Q:=submatrix(M,0,2,0,2) R:=submatrix(M,0,2,3,5)
则 A =R Q ⋅ （可验算之）
例8求例4.3的矩阵的特征值、特征向量。

解 令 A:= ⎥⎥⎥⎦
⎤
⎢⎢⎢⎣⎡----110111011 得 eigenvals(A)= ⎥⎥⎥⎦
⎤
⎢⎢⎢⎣⎡-414.21414.0 eigenvecs(A)= ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---5.0707.05.0707.00707.05.0707.05.0
151
即 三个特征值为 -0.414, 1 ,2.414，矩阵eigenvecs(A)的各列依次为相应的特征向量。

4. 插值和曲线拟合
4.1 代数多项式插值
函数linterp(X,Y,x)可求出插值点x 处的线性插值，其中X 、Y 是向量，分别存放插值基点及其对应的函数值。

例9 已知函数表5.4 ,用线性插值求f(10.5)的近似值。

解 令X 0:=10 X 1:=11 Y 0:=2.303 Y 1:=2.398 x:=10.5
所求线性插值为 linterp(X,Y ,x)=2.351
一般，已知n+1个插值基点及其对应的函数值，求插值点处的n 次插值，可以利用条件函数if (a,b,c)自行定义一个拉格郎日插值函数。

if 函数的值是：当表达式a 为真时if 函数值为b ，否则为c 。

看下例：
例10 已知函数表5.5，试用拉格郎日插值法计算f(1.5)的近似值。

解 自定义拉格郎日插值函数为
Lagr(X,Y,n,x):= 1),x x x x i,if(j Y n 0j j i j n 0i i ∏∑==--≠⋅
令 X:=[0 1 2]T Y:=[2 -1 4]T n:=2
P2(x):=Lagr(X,Y ,n,x)
则f(1.5)的近似值为 P2(1.5)=0.5
[说明] [0 1 2]T 的输入法：先建立1行3列的矩阵[0 1 2]，再按空格键使编辑线将
它从右侧半包围，呈[0 1 2]| 形，单击数学工具栏上矩阵运算按钮，从工具板选
即得。

4.2 三次样条插值
内建函数cspline(X,Y),pspline(X,Y),lspline(X,Y)用于计算不同边界条件下的三次样条插值函数在各基点的二阶导数。

其中X 、Y 向量分别存放n+1个插值基点和相应的函数值。

这三个函数的输出值均为n+4个元素的向量，其中第四个到第n+4个分量依次给出所求的二阶导数值（头三个分量表示插值类型，供内部函数识别用）。

如需计算插值点x 处的三次样条插值，则将上述函数与interp 函数配合使用。

下面仅以计算自然三次样条插值的lspline 函数为例.
例11 已知y=f(x)的函数表5.11，求自然三次样条插值函数在各点的二阶导数，并求f(3)。

解 令 X:=[1 2 4 5]T Y:=[1 3 4 2]T n:=3
152 S:=lspline(X,Y) m:=submatrix(S,3,n+3,0,0)
则各基点的二阶导数为 m 0=0 m 1=-0.75 m 2=-2.25 m 3=0
令 f(x):=interp(S,X,Y ,x)
则 f(3)=4.25
4.3 曲线拟合
例12 设已知数据表5.13，求经验直线拟合这组数据。

解 令 X:=[2 4 6 8]T Y:=[2 11 28 40]T
c 0:=intercept(X,Y) c 1:=slope(X,Y)
得 c 0=-12.5 c 1=6.55
故经验直线为 F(x):=c 0+c 1·x
此例中函数intercept 和slope 的作用是显而易见的。

若用m 次代数多项式作曲线拟合，可调用regress 函数，将它与interp 函数配合可得到所求的多项式函数。

regress 函数的调用形式是
regress(X,Y,m)
其中X 、Y 为存放n+1个数据点的两个向量，m 为拟合多项式的次数，函数的输出值是 n+4个元素的向量，其第4到n+4个元素为所求多项式的系数（按升幂次序排列）。

例13 已知数据表5.14,求二次多项式拟合这组数据。

解 令 X:=[0 1 2 3 4 5]T Y:=[5 2 1 1 2 3]T
n:=5 m:=2 S:=regress(X,Y ,m)
则所求多项式为 F(x):=interp(S,X,Y ,x)
此式可具体用于计算。

若要看到多项式的系数，可令
c:=submatrix(S,3,n+3,0,0)
则有 c T =[4.714 -2.786 0.5]
这样就可以写 F(x):=4.714-2.786x+0.5x 2
例14设在数据表 5.14中，y 的对应值改为9,7,3,5,5,3，试用最小二乘法求形如221
0x c x 11c c g(x)+++=的拟合函数，并作图。

解 令X 同上例, Y:=[9 7 3 5 5 3]T
F(x):=[1 x
11+ 2x ]T C:=linfit(X,Y ,F) （linfit 函数的值是一个由拟合函数的系数构成的列向量） 则令 g(x):=C ·F(x)（系数向量C 与函数向量F(x)作内积）
即为所求。

再为作图做准备，令
I:=0..5 (“..”必须用分号“;”输入)
x:=0,0.1..5（表示x 取0,0.1,0.2,…,5,步长为5.）
153 图形如图6所示。

图中记号×表示已知的数据点。

作图步骤为：（ 表示键入空格）
（1）键入@，出现作图区。

在坐标框
下面占位符处键入“X[i ,x ”，左面占位符
键入“Y[i ,x ”。

（2）双击该图，在弹出的对话框内单
击Traces 子菜单，可见到各曲线的当前设
置栏，其项目有
Legend,Symbol,Line,Color,Type,Weight 共
6项。

当前设置栏的下面有各单项选框，
单击下拉按钮可进行后5项的选择。

(3）单击trace1的当前设置，依次在Symbol 等5项下面的选框中选择x’s,solid,red,points,1。

(4）单击trace2的置，类似（5）的做法，在后5项的选框中依次选none,solid,blue,lines,1。

单击确定。

5. 定积分数值计算
定积分的计算非常方便,只须利用数学工具栏上的微积分工具板,单击定积分按钮,便可产生定积分符号模板，在占位符中填入上下限、被积函数、积分变量，再键入“=”号，便可得到计算结果。

例15 求下列定积分的值：
（1）dx x ⎰-202sin 4π (2) dx e x ⎰-102 (3) dx x ⎰+1
0214 (4)dx x x ⎰++1021)1ln(
解 先修改数据显示格式，使能显示到5位小数。

（1）dx sin(x)42π02⎰
-=2.93492 (2) dx e 10x 2⎰-=0.74682 (3) dx x 14102⎰+=3.14159 (4) dx x 1x)ln(1102⎰++=0.2722
说明：数学表达式x 2sin 在键入时是：sin(x) ^2，显示为2sin(x)。

2
x 1x)ln(1++是键入ln(1+x)后按空格键，直至编辑线从右、下部将整个式子半包围，再键入/1+x^2。

第(4)个积分是0.27220，但末尾的0不显示，如想显示，可修改数据显示格式设置。

若被积函数没有具体表达式，只有函数表，那也很容易构造复合求积公式进行计算。

例16根据数据表6.3，分别用复合梯形公式和复合抛物线公式计算定积分I=dx y ..⎰1
1502。

图6 曲线拟合
154 解 被积函数是x 的函数，但没有给出具体表达式，直接用定积分计算工具不能计算。

令 y:=[0.4804 0.5669 0.6490 0.7262 0.7985 0.8658 0.9281]T
a:=0.5 b:=1.1 n:=6 k:=0..n z k :=(y k )2
用复合梯形公式为
I:= )z 2z (z n 2a b 1n 1
i i n 0∑-=⋅++⋅⋅- 结果为 I=0.32032318 （设数据显示格式为8位小数。

）
用复合抛物线公式为
m:= 2n I:= )z 2z 4z (z n 3a b 1m 1i 2i 1m 0i 1i 2n 0∑
∑-=-=+⋅⋅+⋅++⋅⋅- 结果得 I=0.32010557
(k:=0..n 的输入法是: k:0;n )
6. 求解一阶常微分方程初值问题
求解一阶常微分方程初值问题
b x x a y x y y x f y ≤≤===000,)(,),('
最简便的办法是利用内建的龙格-库塔函数
rkfixed(y,a,b,n,F)
其中y 是存放所求数值解的向量，其第一个分量y 0必须预先赋以初值。

a,b 为求解区间的端点，n 为区间等分数，F 为方程的右端函数f(x,y)，但必须先定义为
F(x,y):=f(x,y 0)
函数rkfixed(y,a,b,n,F)的函数值是一个两列矩阵，其中列出i i y ,x 的对应数值。

例17 求解初值问题
2.0,10,1)0(,2'=≤≤=-=h x y y x y y
解 令 a:=0 b:=1 n:=5
f(x,y):= y
2x y - y 0:=1 F(x,y):=f(x,y 0) 取数据格式为显示到5位小数。

结果是
155 rkfixed(y,a,b,n,F)= ⎥
⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡73214.1161251.18.048328.16
.034167.14.018323.12.010 也可令 R:=rkfixed(y,a,b,n,F) X:=><0R Y:=><1R
（><0R 的上标不能用Shift+6，必须是键入R 后用Ctrl+6，或者在数学工具栏的矩阵工具板单击><M ，再填入0。

它表示取矩阵R 的第1列。

><1R 表示取第2列，余类推。

）
令 k:=0..length(X)-1 (length(X)的函数值为向量X 的元素个数)
结果为(在同一行分别键入X[k= 和Y[k= 便可显示如下表格，适当移动还可使两表靠拢。

)。