19.8.4 交互程序系统和计算机代数系统的应用

19.8.4.1 Matlab

商用程序系统 Matlab (矩阵实验室) 是一个数学求解公式问题的交互环境, 同时也是科学技术计算的高级原本语言. 优先设置线性代数的问题和算法. Matlab 很方便地将实现数值程序及数据和结果的先进的图形表示统一起来. 根据 IEEE 标准 (表 19.4) 计算大多以双精度浮点数运行. 能与 Matlab 兼容的进一步的选择是可免费下载的 Scilab 和 Octave 系统.

1. 函数概述

下面是 Matlab 中程序和函数的简要概述.

Allgemeine mathematische Funktionen

1. Trigonometrie 5. Koordinatentransformationen 2. Exponentialfunktion, Logarith- 6. Rundung und gebrochener Teil men 3. Spezielle Funktionen 7. Diskrete Mathematik 4. Komplexe Arithmetik 8. Mathematische Konstanten

1. Manipulation von Feldern und

Matrizen

2. Spezielle Matrizen

3. Matrix-Analysis (Normen, Kondi-

tion)

4. Lineare Gleichungssysteme

---

5. Eigenwerte und Singulärwerte

6. Matrix-Faktorisierungen

7. Matrixfunktionen

8. Verfahren für dünn besetzte Ma-trizen

---

Numerische Verfahren

1. Berechnung statistischer Daten 7. Ermittlung der konvexen Hülle 2. Korrelation und Regression 8. Numerische Integration 3. Diskrete Fouriertransformation 9. Gewöhnliche Differentialgleichun- 4. Polynome und Splines gen 5. Ein- und mehrdimensionale Inter- 10. Partielle Differentialgleichungen polation 11. Nichtlineare Gleichungen 6. Triangulationen und Zerlegungen 12. Minimierung von Funktionen

此外, Matlab 有许多用来求解不同类型的特殊数学问题的程序包, 即所谓工具箱. 作为某些例子, 这里提到曲线拟合、滤波器、商业数学、时间序列分析、信号和模式处理、神经网络、优化、偏微分方程、样条、统计和小波.

在下面的段落中, 以简单的例子来说明 Matlab 的应用. 在后面 Mathematica 和 Maple 的数值应用中, 将同样讨论这些问题.

2. 数值线性代数

在开启 Matlab 命令提示符 “>>” 后, 命令窗口显示等待接受命令. 若命令不以分号 “; ” 结束,则在命令窗口出现结果. 求解线性方程组 $\mathbf{A}\underset{―}{\mathbf{x}}=\underset{―}{\mathbf{b}}$ 的基本命令 (参见第 1242 页 19.2.1) 是反斜杠算符 “\”.

$◼$ 给定矩阵 $\mathbf{A}=\left(\begin{array}{lll}1& 0& 3\\ 2& 1& 1\\ 1& 2& 3\end{array}\right)$ 和向量 $\underset{―}{b}={(-2,3,2)}^{\mathrm{T}}$ . 输入

$A=\left[\begin{array}{lll}1& 0& 3\end{array};\begin{array}{lll}2& 1& 1\end{array};\begin{array}{lll}1& 2& 3\end{array}\right],b=\left[-2;3;2\right],x=A\setminus b,\mathrm{norm}\left(A\ast x-b\right)$

输出为

$$\begin{array}{l}A=\begin{array}{rrrrrrrrr}1& 0& 3& & & -2& & & 1.0000\\ 2& 1& 1& b=& 3& x=& 2.0000& \text{ ans }=8.8818& \\ 1& 2& 3& & & 2& & -1.0000& \end{array}\text{ ans }=8.8818\mathrm{e}-016\end{array}$$

残量的欧几里得范数表明, 得到的解 (并未显示所有的数字) 以允许的计算机浮点数表示的精度满足方程组.

若矩阵 $\mathbf{A}$ 是非奇异的方阵,则线性方程组有唯一解. 由反斜杠算符 “\” 用通常的列主元高斯消元法,得三角分解 $\mathbf{PA}=\mathbf{LA}$ (参见第 1242 页 19.2.1.1). - 矩阵 $\mathbf{A}$ 的三角分解也可由输入

$$>>\left[\mathrm{L},\mathrm{R},\mathrm{P}\right]=\mathrm{lu}\left(\mathrm{A}\right)$$

实现. 得到输出

$$\mathrm{L}=\left(\begin{array}{rrr}1.0000& 0& 0\\ 0.5000& 1.0000& 0\\ 0.5000& -0.3333& 1.0000\end{array}\right),\mathrm{R}=\left(\begin{array}{rrr}2.0000& 1.0000& 1.0000\\ 0& 1.5000& 2.5000\\ 0& 0& 3.3333\end{array}\right),$$$$P=\left(\begin{array}{lll}0& 1& 0\\ 0& 0& 1\\ 1& 0& 0\end{array}\right)$$

(这里为了避免混乱用括号表示矩阵).

反斜杠算符首先测试系数矩阵 $\mathbf{A}$ 的性质,若 $\mathbf{A}$ 是一个三角矩阵排列,则求解相应的梯队形式. 若是对称矩阵 $\mathbf{A}$ ,则应用楚列斯基方法 (参见第 1245 页 19.2.1.2).

若系数矩阵 $\mathbf{A}$ 的条件数太大,则求解时可能产生数值问题. 因此在 Matlab 运行期间, 估计条件数的倒数值, 若此数过小, 给出警告提示.

$n=13$ 阶的希尔伯特矩阵可作为例子,其中 $h_{ik}=1/\left(i+k-1\right)$ .

$\gg \mathrm{x}=\mathrm{hilb}\left(13\right)\setminus$ ones(13,1);

Warning: Matrix is close to singular or badly scaled.

Results may be inaccurate. RCOND $=2.409320\mathrm{e}-017$ .

在超定线性方程组的情况下, 相应的线性拟合问题由正交化过程处理, 即由正交变换得到 $QR$ 分解

$$\mathbf{A}=\mathbf{QR}$$

(参见第 1246 页 19.2.1.3).

$A=\left[\begin{array}{lll}1& 0& 3\end{array};\begin{array}{lll}2& 1& 1\end{array};\begin{array}{lll}1& 2& 3\end{array};\begin{array}{lll}1& 1& -1\end{array}\right];b=\left[-2;3;2;1\right];x=A\setminus b,$

$$\mathrm{norm}\left(\mathrm{A}\ast \mathrm{x}-\mathrm{b}\right)$$

0.4673

$\mathrm{x}=1.4393$ ans $=2.0508$

-0.4953

$$>>\left[\mathrm{Q},\mathrm{R}\right]=\mathrm{qr}\left(\mathrm{A}\right)$$$$Q=\left(\begin{array}{rrrr}-0.3780& -0.4583& 0.6466& 0.4785\\ -0.7559& -0.2750& -0.2321& -0.5469\\ -0.3780& 0.8250& 0.4145& -0.0684\\ -0.3780& 0.1833& -0.5969& 0.6836\end{array}\right),$$$$R=\left(\begin{array}{rrr}-2.6458& -1.8898& -2.6458\\ 0& 1.5584& 0.6417\\ 0& 0& 3.5480\\ 0& 0& 3.5480\end{array}\right)$$

反斜杠算符也对欠定和缺秩的线性方程组给出了有意义的结果. 如何处理大型稀疏矩阵的方法的详细介绍可见相应的 Matlab 文件和说明 [19.15],[19.30].

3. 方程的数值求解

在 Matlab 中用行向量 $\left(a_n,a_{n-1},\cdots ,a_1,a_0\right)$ 表示多项式

$$p\left(x\right)=a_n{x}^n+a_{n-1}{x}^{n-1}+\cdots +a_{1}x+a_0$$

的系数. 有多个函数可用来处理多项式.

$◼$ 作为例子,确定多项式 $p\left(x\right)=x^6+3x^2-5$ 在 1 处的值、导数 (即导数多项式的系数向量) 及根.

$>>\mathrm{p}=\left[\begin{array}{lllllll}1& 0& 0& 0& 3& 0& -5\end{array}\right];$

$>>$ polyval(p,1)ans $=-1$

$\gg$ polyder(p)ans $=600060$

$\gg$ roots(p) ans $=0.8673+1.1529\mathrm{i},0.8673-1.1529\mathrm{i},1.0743$ ,

$$-0.8673+1.1529i,-0.8673-1.1529i,-1.0743$$

根作为相应伴随矩阵的特征值被确定. 命令 fzero 用来求非线性标量方程的近似解.

$◼$ 计算方程 ${\mathrm{e}}^{-x^3}-4x^2=0$ 的三个解.

$>>\mathrm{fzero}\left(\left(a\right)\exp \left(-x^{\wedge }\right)-4\ast x^{\wedge }2,1\right)$ ans $=0.4741$

$\gg \mathrm{fzero}\left(@\left(x\right)\exp \left(-x^{\wedge }\right)-4\ast x^{\wedge }2,0\right)$ ans $=-0.5413$

$\gg \mathrm{fzero}\left(@\left(x\right)\exp \left(-x^{\wedge }\right)-4\ast x^{\wedge }2,-1\right)$ ans $=-1.2085$

在命令 fzero 下输入方程作为未命名的函数. 显然它依赖于第二个变量的给定初值, 其解是近似的. 迭代过程用二分法和试位法的组合完成 (参见第 1236 页 19.1.1.3).

4. 插值

基于给定数据的函数拟合可以通过插值 (参见第 1276 页 19.6.1) 或最佳函数逼近 (参见第 1280 页 19.6.2.2) 得到. 在 Matlab 中, 命令 plot 是以图形表示数据集合的最简单的方式. 在自打开图表窗口的菜单中包含编辑图形 (线型、符号、标题和说明) 的工具, 用于输出和打印基本拟合.

基本拟合是工具箱的子程序, 提供了多种插值方法和不同阶的最佳逼近多项式. 由函数 “interp1” 和 “polyfit” 来实现.

$◼$ 根据输入

$\gg$ plot $([1.70045,1.2523,0.638803,0.423479,0.249091,0.160321,$

$0.0883432,0.0570776,0.0302744,0.0212794]);$

数据值被置于数据位置 $1,2,\cdots ,10$ 且用图形表示. 图 $19.19\left(\mathrm{a}\right)$ 展示了相应的三次插值样条和四阶最佳逼近多项式的数据集.

函数 interp2 提供在二维矩形网格上给定数据作插值的适当方法 (参见第 1295 页 19.7.2.1). 对于非规则分布的数据作插值调用 griddata.

$◼$ 命令序列

$>>\left[X,Y\right]=$ meshgrid $\left(-2:1:2\right);F=4-\mathrm{sqrt}\left(16-X\cdot 2-Y\cdot 2\right)$ ;

$\gg \left[\mathrm{Xe},\mathrm{Ye}\right]=$ meshgrid $\left(-2:0.1:2\right);\mathrm{S}=\mathrm{interp}2\left(\mathrm{X},\mathrm{Y},\mathrm{F},\mathrm{Xe},\mathrm{Ye},{}^{\prime }{\mathrm{spline}}^{\prime }\right)$ ;

$>>\mathrm{surf}\left(\mathrm{Xe},\mathrm{Ye},\mathrm{S}\right)$

$\gg$ hold on; stem3(X, Y, F,'fill')

实现了函数 $f\left(x,y\right)=\sqrt{16-x^2-y^2}$ 在给定网格上的二元三次样条插值. 插值样条在较细的矩形网格上求值. 图 19.19(b) 描绘了插值函数而且也显示了数据点.

5. 数值积分

在 Matlab 中程序 quad 和 quadl 可用于数值积分. 这两个程序基于自适应选取步长递归应用插值求积. Quad 基于辛普森公式, quadl 应用高阶洛巴托公式 (参见第 1253 页 19.3.2). 在被积函数充分光滑的情况及较高的精度要求下, quadl 比 quad 更有效.

作为第一个例子,考虑近似计算定积分 $I={\int }_0^1\frac{\sin x}{x}\mathrm{d}x$ (正弦积分参见第 681 页 8.2.5,1.).

$\gg$ format long; [I, fwerte]=quad(@(x)(sin(x)./x),0,1)

Warning: Divide by zero. $>\mathrm{In}@\left(x\right)\left(\sin \left(x\right)./x\right)\mathrm{In}$ quad at 62

$I=0.94608307007653$ fwerte $=14$

$\gg$ format long; [I, fwerte]=quadl((@(x)(sin(x)./x),0,1) Warning: Divide by zero. $>\mathrm{In}@\left(x\right)\left(\sin \left(x\right)./x\right)\mathrm{In}$ quadl at 64 $I=0.94608307036718$ fwerte $=19$

两个程序显然都认知被积函数在区间左端的不连续性, 但得到积分的近似值并无困难. 基于同一例子在 1257 页 19.3.4.2 中的结果, 函数求值的数目看来很大, 但可自适应递归确定.

$\gg$ format long; $\left[I,\text{fwerte}\right]=\mathrm{quad}\left(Q\left(x\right)\left(\sin \left(x\right)./x\right),0,1,1e-14\right)$

$I=0.94608307036718$ fwerte $=258$

$\gg$ format long; [I, fwerte]=quadl(( $@\left(x\right)\left(\sin \left(x\right)./x\right),0,1,1\mathrm{e}-14$ )

$I=0.94608307036718$ fwerte $=19$

(这里没有重复警告信息). 作为进一步的论证给出确定的精度为 ${10}^{-14}$ (容忍误差为 ${10}^{-6}$ ),在这一情况下显然可见 quad1 的优越性.

$◼$ 确定定积分 $I={\int }_{-1000}^{1000}{\mathrm{e}}^{-x^2}\mathrm{d}x$ . >> format long; [I, fwerte]=quad(@(x)(exp(-x.^2)),-1000,1000,1e-10) $I=1.77245385094233$ fwerte $=585$

$\gg$ format long; [I, fwerte]=quadl((@(x)(exp(-x.^2)),-1000,1000,1e-10) $I=1.77245385090571$ fwerte $=768$

被积函数在积分区间的非常宽广的部分形状平坦,而在 $x=0$ 处有相当陡峭的峰值, 使得 quad 在这一情况下较好.

6. 微分方程的数值解

Matlab 提供了许多数值求解一阶常微分方程初值问题的程序. 标准的程序为 ode45, 其中应用带自适应步长选取的 4 阶和 5 阶龙格-库塔法 (参见第 1260 页 19.4.1.2). 为达到高精度, 使用预估-校正型线性多步法的程序 ode113 更有效 (参见第 1262 页 19.4.1.4). 此外, 还有些程序对刚性微分方程组特别有效 (参见第 1264 页 19.4.1.5,4.).

$◼$ 用龙格-库塔法 (参见第 1260 页 19.4.1.2) 求解区间 $\left[0,1\right]$ 上的问题 $y^{\prime }=\frac{1}{4}(x^2+$ $y^2),y\left(0\right)=0$ ,输入:

$\gg \left[x,y\right]=\mathrm{ode}45\left(\left(x,y\right)\left(\left(x{.}^2+y{.}^2\right)./4\right),\left[\begin{array}{ll}0& 1\end{array}\right],0\right);\mathrm{plot}\left(x,y\right)$

得到的解如图 19.20(a) 所示.

在区间 $0\le t\le 50$ 上求解特殊的洛伦茨方程组 (参见第 1153 页 17.2.4.3)

$$x_1^{\prime }=10\left(x_2-x_1\right),x_2^{\prime }=28x_1-x_2-x_1{x}_3,x_3^{\prime }=x_1{x}_2-\frac{8}{3}{x}_3$$

其初值条件为 $x\left(0\right)={(0,1,0)}^{\mathrm{T}}$ . 使用如下命令:

$>>\left[t,x\right]=\mathrm{ode}45(@\left(t,x\right)([10\ast \left(x\left(2\right)-x\left(1\right)\right);$

$$28\ast x\left(1\right)-x\left(2\right)-x\left(1\right)\ast x\left(3\right);x\left(1\right)\ast x\left(2\right)-8\ast x\left(3\right)/3]),\left[050\right],\left[0;1;0\right]);$$

$>>\mathrm{plot}\left(x\left(:,1\right),x\left(:,3\right)\right)$

最后的命令给出了在 $x_1,x_3$ 平面的相图 (图 19.20(b)).

19.8.4.2 Mathematica

1. 求解数学问题的工具

计算机代数系统 Mathematica 是能用来求解大量数值数学问题的非常有效的工具. Mathematica 的数学程序完全不同于符号计算. 根据某些预先给定的原则, 类似于图像表示的情况, Mathematica 确定所考虑函数的数值表, 并利用这些数值确定要求的解. 因为点数必须是有限的, 这可能是一个 “坏” 性态函数的问题. 尽管 Mathematica 在问题区域内试图选取更多的节点, 我们必须假定在所考虑区域上的某种连续性. 这可能在最终的结果中产生误差. 建议尽可能使用所考虑问题的许多信息, 且若有可能, 在符号形式下进行计算, 即使这仅对子问题有可能.

在表 19.5 中, 我们介绍了用于数值计算的算子.

NIntegrate	计算定积分
NSum	计算和 $\sum _{i=1}^{n}f\left(i\right)$
NProduct	计算积
NSolve	数值计算求解代数方程
NDSolve	数值计算求解微分方程

在启动 Mathematica 的 "Prompt" 之后显示 $In\left[l\right]\text{:=}$ ; 这表示系统准备接受输入. Mathematica 用 out [l] 表示相应的输出. 一般在行中 $\mathrm{In}\left[n\right]\text{:=}$ 表示输入, $out\left[n\right]\text{:=}$ 表示 Mathemetica 返回的答案. 表达式中的箭头 $\to$ 表示 $x$ 用值 $a$ 代替.

2. 曲线拟合与插值

(1) 曲线拟合 Mathematica 可用最小二乘法选取函数拟合数据集 (参见第 609 页 6.2.5), 且平均逼近离散问题 (参见第 1280 页 19.6.2.2). 一般的指令为

$$\begin{array}{}\text{(19.291)}& \mathrm{Fit}\left[\left\{y_1,y_2,\cdots \right\},\text{ funkt },x\right]\end{array}$$

这里数值 $y_i$ 组成数据表, funkt 为用于数据拟合的可选函数表, $x$ 表示独立变量的相应区域. 若 funkt 选定,例如为 Table $\left[x^{\wedge }i,\{i,0,n\}\right]$ ,则是用 $n$ 次多项式进行拟合.

• 设给定如下数据表:

$$\mathrm{In}\left[1\right]\text{:=}l=\{1.70045,1.2523,0.638803,0.423479,0.249091,0.160321,0.0883432\text{,}$$$$0.0570776,0.0302744,0.0212794\}$$

输入

$$\mathrm{In}\left[2\right]\text{:=}f1=\mathrm{Fit}\left[l,\left\{1,x,x^{\wedge }2,x^{\wedge }3,x^{\wedge }4\right\},x\right]$$

若 $l$ 的元素对应 $x$ 赋值 $1,2,\cdots ,10$ . 其结果为如下四次逼近多项式

$$\text{ Out }\left[2\right]=2.48918-0.853487x+0.0998996x^2-0.00371393x^3-0.0000219224x^4$$

用命令

$$\mathrm{In}\left[3\right]\text{:=}\mathrm{Plot}\left[\mathrm{ListPlot}\left[l,\{x,10\}\right],f1,\{x,1,10\},\text{ AxesOrigin }\to \{0,0\}\right]$$

可得数据表示及图 19.21(a) 给出的近似曲线.

这完全满足给定的数据. ${\mathrm{e}}^{1-0.5x}$ 的级数展开取了前四项.

(2) 插值 Mathematica 提供特殊的算法确定插值函数. 它们表示为类似于纯函数构成的所谓插值函数类. 表 19.6 给出了其使用指南. 取代单独的函数值 $y_i$ ,可在给定点给出一列函数值与指定的导数值.

$◼$ 用 $\mathrm{In}\left[3\right]\text{:=}\mathrm{Plot}\left[\mathrm{Interpolation}\left[l\right]\left[x\right],\{x,1,10\}\right]$ 得图 19.21(b). 显然 Mathemat-ica 给出了精确对应于数据表的结果.

表 19.6 插值命令

Interpolation $\left[\left\{y_1,y_2,\cdots \right\}\right]$	对 $i$ 为整数, $x_i$ 取 $y_i$ 值,给出近似函数
Interpolation $\left[\left\{\left\{x_1,y_1\right\},\left\{x_2,y_2\right\},\cdots \right\}\right]$	对点列 $\left(x_i,y_i\right)$ 给出近似函数

3. 多项式方程的数值解

如同在第 1348 页 20.3.2.1 中所示, Mathematica 能数值确定多项式的根. 命令是

$$\text{NSolve}\left[p\left(x\right)==0,x,n\right]$$

其中 $n$ 为计算需达到的精度. 若忽略 $n$ ,则计算达到机器精度. 若输入 $m$ 次多项式, 则能得到全部 $m$ 个根.

$$In\left[1\right]\text{:=}\text{NSolve}\left[x^{\wedge }6+3x^{\wedge }2-5==0\right]$$$$\text{Out}\left[1\right]=\{x->-1.07432\},\{x->-0.867262-1.15292\mathrm{I}\}\text{,}$$$$\{x\to -0.867262+1.15292\mathrm{I}\},\{x\to 0.867262-1.15292\mathrm{I}\},$$$$\{x\to 0.867262+1.15292\mathrm{I}\},\{x\to 1.07432\}$$

4. 数值积分

Mathematica 为数值积分提供了程序 NIntegrate. 不同于符号方法, 它用被积函数的数值表工作. 考虑以两个非正常积分为例 (参见第 673 页 8.2.3).

$A:\mathrm{In}\left[1\right]\text{:=}$ NIntegrate $\left\lbrack {\operatorname{Exp}\left\lbrack {-{x}^{ \land }2}\right\rbrack ,{ x, - \operatorname{Infinity},\operatorname{Infinity}} }\right\rbrack \operatorname{Out}\left\lbrack 1\right\rbrack =$1.77245.

$\mathbf{B}:\mathrm{In}\left[2\right]\text{:=}$ NIntegrate $\left[1/x^{\wedge }2,\{x,-1,1\}\right]$

Power::infy: Infinite expression $\frac{1}{0}$ encountered.

NIntegrate::inum: Integrand ComplexInfinity is not numerical at $\{x\}=\{0\}$ .

在例 B 中, Mathematica 认知在 $x=0$ 处不连续并给出警告. Mathematica 使用问题区域内大量点的数值表, 并认知极点. 然而, 答案可能仍然是错的.

Mathematica 使用某些预定的选项进行数值积分, 而在某些特殊情况下是不充分的. 可以用参数 MinRecursion 和 MaxRecursion 确定 Mathematica 在问题区域工作的最小和最大递归步数. 默认选项为 0 和 6 . 若此值增加, 则即使 Mathematica 运行变慢, 但结果更好.

$In\left[3\right]\text{:=}$ NIntegrate $\left[\mathrm{Exp}\left[-x\wedge 2\right],\{x,-1000,1000\}\right]$ ,因为积分区域太大, Math-ematica 找不到 $x=0$ 的峰值,答案是

NIntegrate::ploss:

数值积分因失去精度停止. 无法达到要求的精度或准确度的目标; 怀疑有如下情况: 高振荡的被积函数或积分的真值为 0 .

$$\mathrm{Out}\left[3\right]=1.34946\cdot {10}^{-26}$$

若要求为

$$\mathrm{In}\left[4\right]\text{:=}\text{NIntegrate}[\mathrm{Exp}\left[-x^{\wedge }2\right],\{x,-1000,1000\}\text{,}$$$$\text{MinRecursion}\to 3\text{, MaxRecursion}\to 10]$$

则结果是

$$\text{Out}\left[4\right]=1.77245$$

类似地, 可用命令

$$\begin{array}{}\text{(19.292)}& \mathtt{\text{NIntegrate}}\left[\text{fun,}\left\{x,x_a,x_1,x_2,\cdots ,x_e\right\}\right]\end{array}$$

得到接近于积分真实值的结果. 给出积分上下限间的奇点 $x_i$ ,可使 Mathematica 求值更准确.

5. 微分方程的数值解

Mathematica 介绍了用 InterpolatingFunction 数值求解常微分方程和方程组的结果, 可得到解在给定区间任意点的数值及解函数的简略图形表示. 最常用的命令见表 19.7.

表 19.7 微分方程数值解命令

$\mathtt{\text{NDSolve}}\left[dgl,y,\left\{x,x_a,x_e\right\}\right]$	在 $x_a$ 和 $x_e$ 之间的区域计算微分方程数值解
$\mathtt{\text{InterpolatingFunction}}[$ liste $][x]$	给出点 $x$ 处的解
Plot $\left[\text{Evaluate}\left[y\left[x\right]/\text{. lois}\right]],\left\{x,x_a,x_e\right\}\right]$	解的图形表示

求解描述重物在有摩擦的介质中运动的微分方程. 二维运动方程为

$$\ddot{x}=-\gamma \sqrt{{\dot{x}}^2+{\dot{y}}^2}\cdot \dot{x},\ddot{y}=-g-\gamma \sqrt{{\dot{x}}^2+{\dot{y}}^2}\cdot \dot{y}.$$

设摩擦正比于速度. $g=10,\gamma =0.1$ ,初值为 $x\left(0\right)=0,y\left(0\right)=10,\dot{x}\left(0\right)=100,\dot{y}\left(0\right)=$ 200,则下面的命令可用来求解运动方程:

$$\mathrm{In}\left[1\right]\text{:=}dg=\mathrm{NDSolve}[\{x^{\prime \prime }\left[t\right]==-0.1\mathrm{Sqrt}\left[x^{\prime }{\left[t\right]}^{\wedge }2+y^{\prime }{\left[t\right]}^{\wedge }2\right]x^{\prime }\left[t\right],y^{\prime \prime }\left[t\right]==-10$$$$-0.1\mathrm{Sqrt}\left[x^{\prime }{\left[t\right]}^{\wedge }2+y^{\prime }{\left[t\right]}^{\wedge }2\right]y^{\prime }\left[t\right],x\left[0\right]==y\left[0\right]==0,x^{\prime }\left[0\right]==100,y^{\prime }\left[0\right]==200\},$$$$\{x,y\},\{t,15\}]$$

Mathematica 通过插值函数给出答案:

$$\text{Out [1]}=\{\{x\to \text{InterpolatingFunction}\left[\{0.,15.\},<>\right]\text{,}$$$$y\to \text{InterpolatingFunction}\left[\{0.,15.\},<>\right]\}\}$$

则解

$In\left[2\right]\text{:=}$ ParametricPlot $\left[\{x\left[t\right],y\left[t\right]\}/.dg,\{t,0,2\},\text{PlotRange}\to \text{All}\right]$

表示为参数曲线 (图 19.22(a)).

NDSolve 接受影响结果精度的多种选择.

计算的精度可由命令 AccuracyGoal 给出. 类似运行命令 PrecisionGoal. 计算过程中 Mathematica 根据所谓的 WorkingPrecision 运行, 在要求更高精度的计算中应增加五个单位.

Mathematica 在所考虑区域内的工作步数规定为 500 . Mathematica 一般在问题区域的邻域增加节点数, 在奇点邻域可能耗尽步数极限. 在这种情况下, 可用 MaxSteps 增加步数, 也可对 MaxSteps 用命令 Infinity.

• 傅科摆 (Foucault) 方程为

$$\ddot{x}\left(t\right)+{\omega }^{2}x\left(t\right)=2\mathrm{\Omega }\dot{y}\left(t\right),\ddot{y}\left(t\right)+{\omega }^{2}y\left(t\right)=-2\mathrm{\Omega }\dot{x}\left(t\right),$$

其中 $\omega =1,\mathrm{\Omega }=0.025$ ,初值条件为 $x\left(0\right)=0,y\left(0\right)=10,\dot{x}\left(0\right)=\dot{y}\left(0\right)=0$ ,其解为

$$\mathrm{In}\left[3\right]\text{:=}\mathrm{dg}3=\mathrm{NDSolve}[\{x^{\prime \prime }\left[t\right]==-x\left[t\right]+0.05y^{\prime }\left[t\right],y^{\prime \prime }\left[t\right]==-y\left[t\right]-0.05x^{\prime }\left[t\right]\text{,}$$$$x\left[0\right]==0,y\left[0\right]==10,x^{\prime }\left[0\right]==y^{\prime }\left[0\right]==0\},\{x,y\},\{t,0,40\}]$$$$\text{Out}\left[3\right]=\{\{x\to \text{InterpolatingFunction}\left[\{0.,40.\},<>\right]\text{,}$$

$y\to$ InterpolatingFunction $\left[\{0.,40.\},<>\right]\}\}$

加上

$$In\left[4\right]\text{:=}\text{ParametricPlot}\left[\{x\left[t\right],y\left[t\right]\}/\cdot dg3,\{t,0,40\},\text{AspectRatio}\to 1\right]$$

得图 19.22(b).

19.8.4.3 Maple

计算机代数系统 Maple 能用内置逼近法求解许多数值数学问题. 计算要求的节点数,可由规定整体变量 Digits 的任意值 $n$ 确定. 但应记住选取比规定值大的 $n$ 会导致计算速度降低.

1. 表达式和函数的数值计算

在启动 Maple 后, 显示符号 “Prompt”, 这表示准备输入. 输入和输出通常连在一行中表示,由箭头算子 “ $\to$ ” 分开.

(1) 算子 evalf 能作为实数求值的包含内置和用户定义函数的表达式的数值, 可用命令

$$\begin{array}{}\text{(19.293)}& \mathrm{evalf}\left(\mathrm{expr},n\right)\end{array}$$

计算. $\mathrm{expr}$ 为其值待定的表达式,参数 $n$ 可选,估计达到 $n$ 位精度. 默认的精度由全局变量 Digits 确定.

$◼$ 求函数 $y=f\left(x\right)=\sqrt{x}+\ln x$ 的数值表. 首先由箭头算子定义函数:

$$>f\text{:=}z\to \mathrm{sqrt}\left(z\right)+\ln \left(z\right);\to f\text{:=}z\to \sqrt{x}+\ln x$$

则用命令 $\mathrm{evalf}\left(f\left(x\right)\right)$ 可得到要求的函数值; 其中 $x$ 应以数值代替.

在 1 和 4 之间步长为 0.2 的函数值表可用命令

$$>\text{for}x\text{from 1 by 0.2 to 4 do print}\left(f\left[x\right]=\mathrm{evalf}\left(f\left(x\right),12\right)\right)\text{od}$$

得到. 这里要求用 12 位字节计算.

Maple 以形如 $f_{\left[3.2\right]}=2095200519181$ 的元素的一栏表格形式给出结果.

(2) 算子 $\mathrm{evalhf}\left(\exp r\right)$ 除了算子 $\mathrm{evalf}$ 外还有算子 $\mathrm{evalhf}$ . 其用法类似于 evalf. 其变量也是实值表达式. 使用计算机上可行的硬件浮点双精度计算, 求符号表达式的数值. 回到 Maple 浮点值. 在大多数情况下使用 evalhf 加速计算, 但同时使用 evalf 和 Digits 会损失定义的精度. 例如, 在 1303 页 19.8.2 的问题中, 可能产生相当大的误差.

2. 方程的数值解

大多数情况下使用 Maple 可以数值求解方程或方程组. 其命令是 fsolve. 其句法为

$$\begin{array}{}\text{(19.294)}& \text{fsolve(eqn, var, option)}\end{array}$$

该命令确定实数解. 若 eqn 是多项式形式, 则结果都是实根. 若 eqn 不是多项式形式, 则很可能 fsolve 仅返回一个解. 可用的选项在表 19.8 中给出.

complex	确定一个复根 (或多项式的所有根)
maxsols $=n$	确定至少 $n$ 个根 (仅对多项式方程)
fulldigits	保证不减少计算中用到的位数
intervall	在给定区间内求根

$◼\mathbf{A}$: 确定多项式方程 $x^6+3x^3-5=0$ 的所有解. 由

$$>eq\text{:=}{x}^{\wedge }6+3\ast x^{\wedge }2-5=0$$

得

$$>\text{fsolve}\left(eq,x\right);\to -1.074323739,1.074323739$$

Maple 仅得到两个实根. 选择 complex, 也得到复根:

$>$ fsolve(eq, x, complex);

$-1.074323739,-0.8672620244-1.152922012\mathrm{I},-0.8672620244+1.152922012\mathrm{I}$ ,

$0.8672620244-1.152922012\mathrm{I},0.8672620244+1.152922012\mathrm{I},1.074323739$

$◼\mathbf{B}$: 确定超越方程 ${\mathrm{e}}^{-x^3}-4x^2=0$ 的两个解. 定义方程

$$>eq\text{:=}\exp \left(-x^{\wedge }3\right)-4\ast x^{\wedge }2=0$$

正根的结果为

$$>\text{fsolve}\left(eq,x\right);\to 0.4740623572$$

由

$$>\text{fsolve}\left(eq,x,x=-2..0\right);\to -0.5412548544$$

Maple 也确定第二个 (负) 根.

3. 数值积分

通常只可能用数值方法计算定积分的值. 当被积函数太复杂, 或原函数不能表示成基本函数时, 正是这种情况. 在 Maple 中确定定积分的命令为 evalf:

$$\begin{array}{}\text{(19.295)}& \mathrm{evalf}\left(\mathrm{int}\left(f\left(x\right),x=a..b\right),n\right)\end{array}$$

Maple 通过使用逼近公式计算定积分.

$◼$ 计算定积分 ${\int }_{-2}^2{\mathrm{e}}^{-x^3}\mathrm{d}x$ . 因为原函数未知,对于积分命令,得到如下答案:

$$>\mathrm{int}\left(\exp \left(-x^{\wedge }3\right),x=-2..2\right);\to {\int }_{-2}^2{\mathrm{e}}^{-x^3}\mathrm{d}x$$

若

$$>\mathrm{evalf}\left(\mathrm{int}\left(\exp \left(-x^{\wedge }3\right),x=-2..2\right),15\right)$$

则答案为: 277.745841695583.

Maple 使用内置逼近法计算有 15 位数字的数值积分.

在某些情况下, 尤其若积分区间太大, 该方法失效. 此时可试用图书馆的另一个逼近程序

$$\text{readlib(’evalf/int’):}$$

该程序使用自适应牛顿法.

输入

$$>\mathrm{evalf}\left(\mathrm{int}\left(\exp \left(-x^{\wedge }2\right),x=-1000..1000\right)\right)$$

得到误差信息. 由

$$>\text{readlib}\left({}^{\prime }\text{evalf/int}{}^{\prime }\right)\text{:}$$$$\text{> {}^{prime }text{evalf/int}{}^{prime }left( {exp left( {-x^{ land }2}right) , x = - {1000}.{.1000},{10},text{\_NCrule}}right) text{;}}$$$$\text{1.772453851}$$

得到正确的结果. 第三个参量规定精度, 最后一个规定逼近法的内部记号.

4. 微分方程的数值解

常微分方程可用 Maple 运算 dsolve 求解. 然而在大多数情况下, 不可能确定封闭形式的解. 在这些情况下必须给定相应的初值条件, 可以数值求解.

为了做到这一点, 使用如下形式的命令 dsolve:

$$\begin{array}{}\text{(19.296)}& \text{dsolve(deqn, var, numeric),}\end{array}$$

其中选项 numeric 作为第三参量. 这里参量 deqn 包含实际的微分方程和初值条件. 运算结果是一个程序,若记为 $\mathrm{f}$ ,则使用命令 $f\left(t\right)$ ,得到独立变量在 $t$ 处的解函数值.

Maple 应用龙格-库塔法得到该结果 (参见第 1260 页 19.4.1.2). 相对误差和绝对误差的默认精度为 ${10}^{-\text{Digits}+3}$ . 用户可通过整体符号_RELERR 和_ABSERR 修改这些默认容许误差. 若在计算中有某些问题, 则 Maple 给出不同的错误信息.

在解 19.4.1.2 龙格-库塔法一节中给出的问题时, Maple 给出:

$$>r\text{:=}\mathrm{dsolve}\left(\{\mathrm{diff}\left(y\left(x\right),x\right)=\left(1/4\right)\ast \left(x^{\wedge }2+y{\left(x\right)}^{\wedge }2\right),y\left(0\right)=0\},y\left(x\right),\text{ numeric }\right);$$$$r\text{:=}{\text{ proc }}^{\prime }\text{dsolve/numeric/result2 }{}^{\prime }\left(x,1592392,\left[1\right]\right)\text{ end }$$

其中

$$>\mathrm{r}\left(0.5\right);\to \{x\left(.5\right)=0.5000000000,y\left(x\right)\left(.5\right)=0.01041860472\}$$

我们能确定解,例如在 $x=0.5$ 处的值.