数值计算程序.docx
《数值计算程序.docx》由会员分享，可在线阅读，更多相关《数值计算程序.docx（19页珍藏版）》请在冰豆网上搜索。
数值计算程序
附录一用Mathcad进行数值计算
Mathcad是世界流行的一个优秀的交互式应用数学软件。
其初版由美国Mathsoft公司于1986年推出，1997年问世的7.0版本可在Windows95/98或WindowsNT下运行。
Mathcad是集数理计算、图形和文字处理等功能于一体的多媒体综合应用软件，它支持OLE2技术，可以和其它应用软件协同工作，还具有网络浏览功能。
同时，它又是一个大众化的数学工具软件，用户只须学会简单的操作，无须记住很多的计算机知识和数学技巧，就可利用它进行各种复杂的数学计算、数据处理、作图和文本编辑。
Mathcad因功能强大且使用简便而享誉全世界，它拥有约100多万用户。
无论对于大学师生还是对于企业界的技术人员和管理者,Mathcad都是解决应用数学问题的好帮手。
下面主要针对前几章所涉及的数学模型问题简要介绍如何使用Mathcad7.0进行常用的科学计算。
1.Mathcad基本用法
Mathcad7.0分标准版，学生版和专业版，其中专业版（Mathcad7Professional）的功能最齐全。
1.1工作页和工具栏
进入Windows95/98后,在”开始/程序”的菜单中单击”Mathcad7Professional”项,即可启动该系统.屏幕显示启动画面后进入用户界面,这时屏幕中央有一个”TipoftheDay”窗口,单击ok后,屏幕显示如图1
图1Mathcad用户界面
在系统标题行下的主菜单中含8个主菜项，即File（文件管理）、Edit（编辑）、View（视图）、Insert（插入）、Format（格式）、Math（数学计算）、Symbolics（符号计算）、Window（窗口管理）、Help（帮助）。
其中大部分菜单功能与WORD97中同名菜单项类似。
单击Help可看到其中有若干项帮助命令，常用的是在线帮助（Mathcadhelp）以及资源中心（ResourceCenter）。
当你打开资源中心后可选择进入浏览（Overview）、教程（Tutorial）、快速参考活页（QuickSheets）等多项帮助内容,获取有针对性的指导及示例。
字体、常用工具栏的大部分功能也同WORD类似。
常用工具栏中的f（x）按钮可用于选择系统的内建函数输入到表达式中去。
Mathcad7有200多个内建函数可供我们使用，如果记住了函数名及其变量表的格式，也可以直接用键盘输入所要的函数。
数学工具栏见图2，它含8个按钮，从左到右依次为：
（1）计算器
（2）逻辑运算（3）函数作图（4）矩阵和向量运算
图2数学工具栏
（5）微积分运算（6）编程工具（7）希腊字母（8）符号运算
用每个按钮可打开一块工具板，从中可选择所需要的数学符号、运算符或作图命令。
Mathcad工作页的排版是二维的，它用起来象一张向右向下都无限伸展的白纸，只是用边界线（虚线）分割不同的页面。
往工作页上键入式子时，一个式子常可以跨页书写，但若无必要则应避免式子跨页。
因为打印时，界线两侧分属不同页面，跨页的式子便被分割开了。
光标在工作页的空白处呈红十字（见图1），可以用箭头键来移动光标。
整个工作页是一个Mathcad文件，当需要保存它时，用主菜单“文件/保存”命令即可，默认的文件名后缀为.mcd。
下次使用该文件时用主菜单“文件/打开”命令便可打开它。
1.2数学区、文本区、作图区
在工作页中可建立数学区（也称计算区）、文本区、作图区。
在工作页的任何空白处键入任一字符都被默认为开辟一个数学区（也可用Insert/MathRegion菜单命令来开辟），并用矩形框围住，红十字光标变成直角状兰色编辑线将字符半包围，如图3。
当光标移出该区后又变成红十字，数学区域框隐去。
在数学区可以方便地键入表达式并进行计算，也可键入汉字，要注意：
（1）每个数学区只能有一个独立的表达式；
（2）一个文件中所有数学区的字体格式及字号均是一致的，需要改变时，只能将光标放入任一数学区内通过主菜单Format/Equation进行统一修改；（3）键入空格的效果是：
或者移动并扩大编辑线包围的范围，或者将所在数学区变成文本区。
对数学区的计算命令有自动和手动两种模式，用主菜单Math/AutomaticCalculation命令转换两种模式。
通常默认为自动计算模式，在信息栏可见到”Auto”信息。
在Auto模式下工作页是”活”的，成为一个广义的“Mathcad程序”，只要修改一处式子，其后与之有计算关系的表达式的计算结果以及相关图形也自动随之改变。
在手动模式下，按F9可下达计算命令。
有时需要将某个数学区“关闭”，即禁止计算，方法是光标放入该区，用菜单Format/Properties/Calculation命令，选择DisableEvaluation并单击OK。
被关闭的数学区右上角出现■记号。
要把被关闭的数学区打开的方法类似，取消DisableEvaluation选择即可。
图3数学区和文本区
文本区可在工作页上任意空白处开辟。
方法有三种：
（1）键入双引号”。
（2）单击主菜单Insert/TextRegion。
（3）键入任一字符后按空格键。
光标在文本区内呈红色竖线，即插入线，这时可见到文本区的边界，是一个带有三个“把手”的矩形框，如图3所示。
文本区一般是不被计算的，但需要时可在文本区内插入数学区。
方法是让红色插入线位于插入处，点主菜单Insert/MathRegion。
插入文本区的数学区可以被计算，除非把它关闭。
作图区的开辟可用Insert/Graph菜单命令，再选择坐标系。
简单的办法是利用数学工具栏上的作图工具板进行直观的选择。
也可用快捷键，如键入@便可开辟一个直角坐标作图区。
不论是数学区、文本区还是作图区，都可以用鼠标左键按住其边框拖拽移动到任何地方（这时鼠标呈手掌形）。
如果想同时拖拽若干个区，只须先同时选中它们即可，方法是用左键按住待选区域附近的空白处，拖出一个较大的虚框，松开左键，该虚框所及之处的所有区域分别被小虚框围住，这就表示被选中了。
被选中的区域除能被移动外，还可被复制、删除等，单击右键可弹出文中菜单进行这些处理。
当作了某些删除或移动后，可以用Ctrl+R键来刷新屏幕。
若需将某些区域横向对齐或竖向对齐，只须选中这些区，再单击常用工具栏上f（x）左侧的两个按钮之一。
1.1在数学区键入和计算表达式
利用Mathcad提供的快捷键和大量的数学工具板，可以非常方便地输入数学表达式，并且与通常书写的形式差不多。
例如按从上到下的次序键入下面左列的字符串（表示键入空格），每个式子键入完毕后将光标移出所在区域或按回车键结束。
显示结果如右列所示。
键入显示（设系统在默认Auto模式下）
x^-5x-5
x||x|
x:
5x:
=5
sin（x）=sin（x）=-0.959
x\=
=2.236
3*x=
=15
3+x/5=
=1.6
3+x/5=
=4
y[0:
1y0:
=1
x!
=x!
=120
要注意“:
=”是赋值号，只须键入“：
”即可产生；单独的“=”是表示计算命令的等号，其右端为机器自动给出的计算结果。
上面第三行将x赋值5，所以后面凡与x有关的算式均自动得出计算结果。
你可将x:
5改为x:
6，看结果会如何变化？
此外，要注意空格键的作用。
上面的y0:
=1表示已定义y为向量且第0个分量赋值1，Mathcad中向量及矩阵的起始下标为0，必要时可用主菜单Math/Options命令，在Array栏中选其它整数作为起始下标。
数据的显示长度默认为到3位小数，若需显示更多位，只须光标位于某数学区内时用主菜项Format/Number命令，在Precision栏中修改DisplayedPrecision的值即可。
一些特殊符号可利用数学工具栏的工具板，特别是微积分符号必须用工具板输入才能在计算时被系统认识。
打开数学工具栏中微积分运算工具板，当光标在数学区或工作页空白处时，单击
、
、
、
按钮之一，便分别显示如图4的算符：
图4特殊数学符号
其中小黑块是占位符，在占位符中分别键入适当字符便可完成式子的输入，当占位符空缺而光标移出该区时占位符呈红色，表示该式尚未写完。
在Auto模式下，系统自动从上到下逐行扫描，每行又从左到右、不受页边限制地逐个扫描数学区并依次执行计算。
所以键入计算公式时要注意先后次序。
下面各节举例说明如何用Mathcad解决常用数学模型的计算问题。
2.求解一元方程
求根函数root（f（x）,x）可用于求解一元方程f（x）=0.调用root函数前须先定义f（x），并对x赋以初始近似值（如果求复根则初始近似值应为复数）。
此外还要给系统常量TOL指定迭代的精度，否则默认TOL的值为10-3。
例1求方程x－cosx=0在0.5附近的根，要求误差不超过10-6。
解令f（x）:
=x-cos（x）x:
=0.5TOL:
=10-6
用主菜单Format/Number命令在Precision栏修改数据显示为8位小数，再调用求根函数得所求的根为
root（f（x）,x）=0.73908513
上例的求解过程中，“:
=”右边是人工键入的，而单独的“=”右边是计算机给出的结果。
以下各例也相同，不再一一指出。
例2利用图形估计方程x－cosx=0的初始近似根。
解设f（x）:
=x-cos（x）再作以下步骤：
图5函数f（x）的图形
（1）在上式下方空白处键入@，出现作图区，内有矩形坐标框。
（2）在坐标框下面占位符处键入x，左面占位符处键入f（x），鼠标左键单击作图区外任一点，函数图形便自动生成。
（3）双击该图弹出对话框，在X-YAxes子菜单的AxesStyle栏中选Crossed，单击确定按钮。
函数图形便如图5所示。
（4）观察曲线与x轴的交点，估计初始近似根可取为0.5。
（5）若需更精确些，单击Format/Graph/Trace，再单击曲线与x轴的交点，出现过此点的十字交叉虚线，交点的坐标在弹出的对话框内，其中x的值可取为初始近似根。
3.线性代数计算
输入矩阵的方法是：
单击数学工具栏的矩阵运算按钮，在工具板中选矩阵形式按钮，弹出对话框，填入行数、列数，单击确定，然后在占位符中键入具体字符即可。
用函数lsolve（A,b）可实现列主元高斯消去法求解方程组Ax=b，结果输出所求的解向量。
例3求解例3.2的线性方程组（设未知量改为x0,x1,x2）。
解令
A:
=
b:
=
x:
=lsolve（A,b）
得xT=[0-22]即x0=0x1=-2x2=2
（xT中的T只能用矩阵运算工具板的MT键入，不能用^键入，否则系统不能正确计算。
）
例4求例3.4中矩阵A的行列式及逆阵。
解A:
=
A-1=
|A|=-2
例5对例3.7的系数矩阵进行选主元的LU分解。
解令A:
=
M:
=lu（A）（调用内建的lu函数）
得M=
（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=
L=
U=
满足
=
（这个等号是逻辑等号，用数学工具栏的逻辑运算工具板键入。
）
可分别计算
和
，验算上述等式。
上例中内建函数submatrix（M,a,b,c,d）是从矩阵M中取行下标a~b，列下标c~d的子矩阵。
例6求例3.9的系数矩阵的LLT分解。
解令A:
=
L:
=cholesky（A）得L=
例7求例4.4的矩阵A的QR分解。
解令A:
=
M:
=qr（A）
得M=
（M由Q和R并列构成）
令Q:
=submatrix（M,0,2,0,2）R:
=submatrix（M,0,2,3,5）
则A=
（可验算之）
例8求例4.3的矩阵的特征值、特征向量。
解令A:
=
得eigenvals（A）=
eigenvecs（A）=
即三个特征值为-0.414,1,2.414，矩阵eigenvecs（A）的各列依次为相应的特征向量。
4.插值和曲线拟合
4.1代数多项式插值
函数linterp（X,Y,x）可求出插值点x处的线性插值，其中X、Y是向量，分别存放插值基点及其对应的函数值。
例9已知函数表5.4,用线性插值求f（10.5）的近似值。
解令X0:
=10X1:
=11Y0:
=2.303Y1:
=2.398x:
=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）:
=
令X:
=[012]TY:
=[2-14]Tn:
=2
P2（x）:
=Lagr（X,Y,n,x）
则f（1.5）的近似值为P2（1.5）=0.5
[说明][012]T的输入法：
先建立1行3列的矩阵[012]，再按空格键使编辑线将它从右侧半包围，呈[012]|形，单击数学工具栏上矩阵运算按钮，从工具板选MT即得。
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:
=[1245]TY:
=[1342]Tn:
=3
S:
=lspline（X,Y）m:
=submatrix（S,3,n+3,0,0）
则各基点的二阶导数为m0=0m1=-0.75m2=-2.25m3=0
令f（x）:
=interp（S,X,Y,x）
则f（3）=4.25
4.3曲线拟合
例12设已知数据表5.13，求经验直线拟合这组数据。
解令X:
=[2468]TY:
=[2112840]T
c0:
=intercept（X,Y）c1:
=slope（X,Y）
得c0=-12.5c1=6.55
故经验直线为F（x）:
=c0+c1·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:
=[012345]TY:
=[521123]T
n:
=5m:
=2S:
=regress（X,Y,m）
则所求多项式为F（x）:
=interp（S,X,Y,x）
此式可具体用于计算。
若要看到多项式的系数，可令
c:
=submatrix（S,3,n+3,0,0）
则有cT=[4.714-2.7860.5]
这样就可以写F（x）:
=4.714-2.786x+0.5x2
例14设在数据表5.14中，y的对应值改为9,7,3,5,5,3，试用最小二乘法求形如
的拟合函数，并作图。
解令X同上例,Y:
=[973553]T
F（x）:
=[1
]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.）
图6曲线拟合
图形如图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）
（2）
（3）
（4）
解先修改数据显示格式，使能显示到5位小数。
（1）
=2.93492
（2）
=0.74682
（3）
=3.14159（4）
=0.2722
说明：
数学表达式
在键入时是：
sin（x）^2，显示为
。
是键入ln（1+x）后按空格键，直至编辑线从右、下部将整个式子半包围，再键入/1+x^2。
第（4）个积分是0.27220，但末尾的0不显示，如想显示，可修改数据显示格式设置。
若被积函数没有具体表达式，只有函数表，那也很容易构造复合求积公式进行计算。
例16根据数据表6.3，分别用复合梯形公式和复合抛物线公式计算定积分I=
。
解被积函数是x的函数，但没有给出具体表达式，直接用定积分计算工具不能计算。
令y:
=[0.48040.56690.64900.72620.79850.86580.9281]T
a:
=0.5b:
=1.1n:
=6k:
=0..nzk:
=（yk）2
用复合梯形公式为
I:
=
结果为I=0.32032318（设数据显示格式为8位小数。
）
用复合抛物线公式为
m:
=
I:
=
结果得I=0.32010557
（k:
=0..n的输入法是:
k:
0;n）
6.求解一阶常微分方程初值问题
求解一阶常微分方程初值问题
最简便的办法是利用内建的龙格-库塔函数
rkfixed（y,a,b,n,F）
其中y是存放所求数值解的向量，其第一个分量y0必须预先赋以初值。
a,b为求解区间的端点，n为区间等分数，F为方程的右端函数f（x,y），但必须先定义为
F（x,y）:
=f（x,y0）
函数rkfixed（y,a,b,n,F）的函数值是一个两列矩阵，其中列出
的对应数值。
例17求解初值问题
解令a:
=0b:
=1n:
=5
f（x,y）:
=
y0:
=1F（x,y）:
=f（x,y0）
取数据格式为显示到5位小数。
结果是
rkfixed（y,a,b,n,F）=
也可令R:
=rkfixed（y,a,b,n,F）X:
=
Y:
=
（
的上标不能用Shift+6，必须是键入R后用Ctrl+6，或者在数学工具栏的矩阵工具板单击
，再填入0。
它表示取矩阵R的第1列。
表示取第2列，余类推。
）
令k:
=0..length（X）-1（length（X）的函数值为向量X的元素个数）
结果为（在同一行分别键入X[k=和Y[k=便可显示如下表格，适当移动还可使两表靠拢。
）
XkYk
0
1
0.2
1.18323
0.3
1.34167
0.4
1.48328
0.5
1.61251
0.6
1.73214
升级会员 