数值计算程序.docx

1、数值计算程序 附录一 用Mathcad进行数值计算Mathcad是世界流行的一个优秀的交互式应用数学软件。其初版由美国Mathsoft公司于1986年推出，1997年问世的7.0版本可在Windows95/98或Windows NT下运行。Mathcad是集数理计算、图形和文字处理等功能于一体的多媒体综合应用软件，它支持OLE2 技术，可以和其它应用软件协同工作，还具有网络浏览功能。同时，它又是一个大众化的数学工具软件，用户只须学会简单的操作，无须记住很多的计算机知识和数学技巧，就可利用它进行各种复杂的数学计算、数据处理、作图和文本编辑。Mathcad因功能强大且使用简便而享誉全世界，它拥有约

2、100多万用户。无论对于大学师生还是对于企业界的技术人员和管理者,Mathcad都是解决应用数学问题的好帮手。下面主要针对前几章所涉及的数学模型问题简要介绍如何使用Mathcad7.0进行常用的科学计算。1. Mathcad 基本用法Mathcad 7.0分标准版，学生版和专业版，其中专业版(Mathcad7 Professional)的功能最齐全。 1.1 工作页和工具栏进入Windows95/98后,在”开始/程序”的菜单中单击”Mathcad 7 Professional”项,即可启动该系统.屏幕显示启动画面后进入用户界面,这时屏幕中央有一个”Tip of the Day”窗口,单击ok

3、后,屏幕显示如图1图1 Mathcad用户界面 在系统标题行下的主菜单中含8个主菜项，即File(文件管理)、Edit(编辑)、View(视图)、Insert(插入)、Format(格式)、Math(数学计算)、Symbolics(符号计算)、Window(窗口管理)、Help(帮助)。其中大部分菜单功能与WORD97中同名菜单项类似。单击Help可看到其中有若干项帮助命令，常用的是在线帮助(Mathcad help)以及资源中心(Resource Center)。当你打开资源中心后可选择进入浏览(Overview)、教程(Tutorial)、快速参考活页(QuickSheets)等多项帮助内

4、容,获取有针对性的指导及示例。字体、常用工具栏的大部分功能也同WORD类似。常用工具栏中的f(x)按钮可用于选择系统的内建函数输入到表达式中去。Mathcad7有200多个内建函数可供我们使用，如果记住了函数名及其变量表的格式，也可以直接用键盘输入所要的函数。数学工具栏见图2，它含8个按钮，从左到右依次为：(1) 计算器 (2)逻辑运算 (3)函数作图 (4)矩阵和向量运算图2 数学工具栏(5)微积分运算 (6)编程工具 (7)希腊字母 (8)符号运算用每个按钮可打开一块工具板，从中可选择所需要的数学符号、运算符或作图命令。 Mathcad工作页的排版是二维的，它用起来象一张向右向下都无限伸展

5、的白纸，只是用边界线(虚线)分割不同的页面。往工作页上键入式子时，一个式子常可以跨页书写，但若无必要则应避免式子跨页。因为打印时，界线两侧分属不同页面，跨页的式子便被分割开了。 光标在工作页的空白处呈红十字(见图1)，可以用箭头键来移动光标。 整个工作页是一个Mathcad文件，当需要保存它时，用主菜单“文件/保存”命令即可，默认的文件名后缀为.mcd。下次使用该文件时用主菜单“文件/打开”命令便可打开它。 1.2 数学区、文本区、作图区在工作页中可建立数学区(也称计算区)、文本区、作图区。在工作页的任何空白处键入任一字符都被默认为开辟一个数学区(也可用Insert/Math Region菜单

6、命令来开辟)，并用矩形框围住，红十字光标变成直角状兰色编辑线将字符半包围，如图3。当光标移出该区后又变成红十字，数学区域框隐去。 在数学区可以方便地键入表达式并进行计算，也可键入汉字，要注意：(1)每个数学区只能有一个独立的表达式；(2)一个文件中所有数学区的字体格式及字号均是一致的，需要改变时，只能将光标放入任一数学区内通过主菜单Format/Equation进行统一修改；(3)键入空格的效果是：或者移动并扩大编辑线包围的范围，或者将所在数学区变成文本区。 对数学区的计算命令有自动和手动两种模式，用主菜单Math/Automatic Calculation命令转换两种模式。通常默认为自动计算

7、模式，在信息栏可见到”Auto”信息。在Auto模式下工作页是”活”的，成为一个广义的“Mathcad程序”，只要修改一处式子，其后与之有计算关系的表达式的计算结果以及相关图形也自动随之改变。在手动模式下，按F9可下达计算命令。 有时需要将某个数学区“关闭”，即禁止计算，方法是光标放入该区，用菜单Format/Properties/Calculation命令，选择Disable Evaluation并单击OK。被关闭的数学区右上角出现记号。要把被关闭的数学区打开的方法类似，取消Disable Evaluation选择即可。 图3 数学区和文本区 文本区可在工作页上任意空白处开辟。方法有三种：(

8、1)键入双引号”。(2)单击主菜单Insert/Text Region。(3)键入任一字符后按空格键。 光标在文本区内呈红色竖线，即插入线，这时可见到文本区的边界，是一个带有三个“把手”的矩形框，如图3所示。 文本区一般是不被计算的，但需要时可在文本区内插入数学区。方法是让红色插入线位于插入处，点主菜单Insert/Math Region。插入文本区的数学区可以被计算，除非把它关闭。 作图区的开辟可用Insert/Graph菜单命令，再选择坐标系。简单的办法是利用数学工具栏上的作图工具板进行直观的选择。也可用快捷键，如键入便可开辟一个直角坐标作图区。 不论是数学区、文本区还是作图区，都可以用鼠

9、标左键按住其边框拖拽移动到任何地方（这时鼠标呈手掌形）。如果想同时拖拽若干个区，只须先同时选中它们即可，方法是用左键按住待选区域附近的空白处，拖出一个较大的虚框，松开左键，该虚框所及之处的所有区域分别被小虚框围住，这就表示被选中了。被选中的区域除能被移动外，还可被复制、删除等，单击右键可弹出文中菜单进行这些处理。当作了某些删除或移动后，可以用Ctrl+R键来刷新屏幕。若需将某些区域横向对齐或竖向对齐，只须选中这些区，再单击常用工具栏上f(x)左侧的两个按钮之一。1.1 在数学区键入和计算表达式利用Mathcad提供的快捷键和大量的数学工具板，可以非常方便地输入数学表达式，并且与通常书写的形式差

10、不多。例如按从上到下的次序键入下面左列的字符串（ 表示键入空格），每个式子键入完毕后将光标移出所在区域或按回车键结束。显示结果如右列所示。键 入 显 示（设系统在默认Auto模式下） x-5 x-5 x| |x| x:5 x:=5 sin(x)= sin(x)=-0.959 x= =2.236 3*x= =15 3+x /5= =1.6 3+x/5= =4 y0:1 y0:=1 x!= x!=120要注意“:=”是赋值号，只须键入“：”即可产生；单独的“=”是表示计算命令的等号，其右端为机器自动给出的计算结果。上面第三行将x赋值5，所以后面凡与x有关的算式均自动得出计算结果。你可将x:5改为x

11、:6，看结果会如何变化？此外，要注意空格键的作用。 上面的y0:=1表示已定义y为向量且第0个分量赋值1，Mathcad中向量及矩阵的起始下标为0，必要时可用主菜单Math/Options命令，在Array栏中选其它整数作为起始下标。 数据的显示长度默认为到3位小数，若需显示更多位，只须光标位于某数学区内时用主菜项Format/Number命令，在Precision栏中修改Displayed Precision的值即可。 一些特殊符号可利用数学工具栏的工具板，特别是微积分符号必须用工具板输入才能在计算时被系统认识。打开数学工具栏中微积分运算工具板，当光标在数学区或工作页空白处时，单击、按钮之一

12、，便分别显示如图4的算符： 图4 特殊数学符号其中小黑块是占位符，在占位符中分别键入适当字符便可完成式子的输入，当占位符空缺而光标移出该区时占位符呈红色，表示该式尚未写完。 在Auto模式下，系统自动从上到下逐行扫描，每行又从左到右、不受页边限制地逐个扫描数学区并依次执行计算。所以键入计算公式时要注意先后次序。 下面各节举例说明如何用Mathcad解决常用数学模型的计算问题。2. 求解一元方程 求根函数root(f(x),x)可用于求解一元方程f(x)=0.调用root函数前须先定义f(x)，并对x赋以初始近似值（如果求复根则初始近似值应为复数）。此外还要给系统常量TOL指定迭代的精度，否则默

13、认TOL的值为10-3。 例1 求方程xcos 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 利用图形估计方程xcos x =0的初始近似根。 解 设 f(x):=x-cos(x) 再作以下步骤： 图5 函数f(x)的图形(1)在上式下方空白处

14、键入，出现作图区，内有矩形坐标框。(2)在坐标框下面占位符处键入x，左面占位符处键入f(x)，鼠标左键单击作图区外任一点，函数图形便自动生成。(3)双击该图弹出对话框，在X-Y Axes子菜单的Axes Style栏中选Crossed，单击确定按钮。函数图形便如图5所示。 (4)观察曲线与x轴的交点，估计初始近似根可取为0.5。(5)若需更精确些，单击Format/Graph/Trace，再单击曲线与x轴的交点，出现过此点的十字交叉虚线，交点的坐标在弹出的对话框内，其中x的值可取为初始近似根。3. 线性代数计算 输入矩阵的方法是：单击数学工具栏的矩阵运算按钮，在工具板中选矩阵形式按钮，弹出对话

15、框，填入行数、列数，单击确定，然后在占位符中键入具体字符即可。用函数lsolve(A,b)可实现列主元高斯消去法求解方程组Ax=b，结果输出所求的解向量。例3 求解例3.2的线性方程组（设未知量改为x0,x1,x2）。解 令 A:= b:= x:=lsolve(A,b)得 xT=0 -2 2 即 x0=0 x1=-2 x2=2（xT中的T只能用矩阵运算工具板的MT键入，不能用键入，否则系统不能正确计算。） 例4 求例3.4中矩阵A的行列式及逆阵。 解 A:= A-1= |A|=-2例5 对例3.7的系数矩阵进行选主元的LU分解。解 令 A:= M:=lu(A) (调用内建的lu函数)得 M=

16、（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中取行下标ab，列下标cd的子矩阵。例6求例3.9的系数矩阵的LLT分解。解 令A:= L:=cholesky(A) 得 L=例7求例4.4的矩阵A的QR分解。解 令 A:= M:=qr(A) 得 M= (M由Q和R并列构成)令 Q:=su

17、bmatrix(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:=10 X1:=11 Y0:=2.30

18、3 Y1:=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):= 令 X:=0 1 2T Y:=2 -1 4T n:=2 P2(x):=Lagr(X,Y,n,x)则f(1.5)的近似值为 P2(1.5)=0.5说明 0 1 2T的输入法：先建

19、立1行3列的矩阵0 1 2，再按空格键使编辑线将它从右侧半包围，呈0 1 2| 形，单击数学工具栏上矩阵运算按钮，从工具板选MT 即得。 4.2 三次样条插值内建函数cspline(X,Y),pspline(X,Y),lspline(X,Y)用于计算不同边界条件下的三次样条插值函数在各基点的二阶导数。其中X、Y向量分别存放n+1个插值基点和相应的函数值。这三个函数的输出值均为n+4个元素的向量，其中第四个到第n+4个分量依次给出所求的二阶导数值（头三个分量表示插值类型，供内部函数识别用）。如需计算插值点x处的三次样条插值，则将上述函数与interp函数配合使用。下面仅以计算自然三次样条插值的l

20、spline函数为例.例11 已知y=f(x)的函数表5.11，求自然三次样条插值函数在各点的二阶导数，并求f(3)。解 令 X:=1 2 4 5T Y:=1 3 4 2T n:=3 S:=lspline(X,Y) m:=submatrix(S,3,n+3,0,0)则各基点的二阶导数为 m0=0 m1=-0.75 m2=-2.25 m3=0令 f(x):=interp(S,X,Y,x)则 f(3)=4.254.3 曲线拟合例12 设已知数据表5.13，求经验直线拟合这组数据。解 令 X:=2 4 6 8T Y:=2 11 28 40T c0:=intercept(X,Y) c1:=slope(

21、X,Y)得 c0=-12.5 c1=6.55故经验直线为 F(x):=c0+c1x此例中函数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 5T Y:=5 2 1 1 2 3T n:=5

22、m:=2 S:=regress(X,Y,m)则所求多项式为 F(x):=interp(S,X,Y,x)此式可具体用于计算。若要看到多项式的系数，可令 c:=submatrix(S,3,n+3,0,0)则有 cT=4.714 -2.786 0.5这样就可以写 F(x):=4.714-2.786x+0.5x2 例14设在数据表5.14中，y的对应值改为9,7,3,5,5,3，试用最小二乘法求形如的拟合函数，并作图。 解 令X 同上例, Y:=9 7 3 5 5 3T F(x):=1 T C:=linfit(X,Y,F) （linfit函数的值是一个由拟合函数的系数构成的列向量）则令 g(x):=C

23、F(x)（系数向量C与函数向量F(x)作内积）即为所求。再为作图做准备，令 I:=0.5 (“.”必须用分号“;”输入) x:=0,0.1.5（表示x取0,0.1,0.2,5,步长为5.） 图6 曲线拟合图形如图6所示。图中记号表示已知的数据点。作图步骤为：（ 表示键入空格） （1）键入，出现作图区。在坐标框下面占位符处键入“Xi ,x”，左面占位符键入“Yi ,x”。 （2）双击该图，在弹出的对话框内单击Traces子菜单，可见到各曲线的当前设置栏，其项目有Legend,Symbol,Line,Color,Type,Weight共6项。当前设置栏的下面有各单项选框，单击下拉按钮可进行后5项的

24、选择。(3）单击trace1的当前设置，依次在Symbol等5项下面的选框中选择xs,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.7468

25、2 (3) =3.14159 (4) =0.2722说明：数学表达式在键入时是：sin(x) 2，显示为。是键入ln(1+x)后按空格键，直至编辑线从右、下部将整个式子半包围，再键入/1+x2。第(4)个积分是0.27220，但末尾的0不显示，如想显示，可修改数据显示格式设置。 若被积函数没有具体表达式，只有函数表，那也很容易构造复合求积公式进行计算。例16根据数据表6.3，分别用复合梯形公式和复合抛物线公式计算定积分I=。解 被积函数是x的函数，但没有给出具体表达式，直接用定积分计算工具不能计算。令 y:=0.4804 0.5669 0.6490 0.7262 0.7985 0.8658 0

26、.9281T a:=0.5 b:=1.1 n:=6 k:=0.n zk:=(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,y

27、0)函数rkfixed(y,a,b,n,F)的函数值是一个两列矩阵，其中列出的对应数值。例17 求解初值问题 解 令 a:=0 b:=1 n:=5 f(x,y):= y0:=1 F(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的元素个数)结果为(在同一行分别键入Xk= 和Yk= 便可显示如下表格，适当移动还可使两表靠拢。)Xk Yk 010.21.183230.31.341670.41.483280.51.612510.61.73214