数值分析上机试验指导书讲解

1、数值计算方法”上机实验指导书 实验一 误差分析 实验 1.1(病态问题) 实验目的 ：算法有“优”与“劣”之分，问题也有“好”与“坏”之别。对数值方法 的研究而言，所谓坏问题就是问题本身对扰动敏感者，反之属于好问题。通过本实验可获 得一个初步体会。 数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及 方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要 付出一些代价(如耗用更多的机器时间、占用更多的存储空间等) 。 问题提出 ：考虑一个高次的代数多项式 20 p(x) ( x 1)(x 2) ( x 20) (x k) (1.1) k1 显然该

2、多项式的全部根为 1,2,20共计 20个，且每个根都是单重的。现考虑该多项式的 一个扰动 p( x) x19 0 (1.2) 其中 是一个非常小的数。这相当于是对( 1.1)中 x19 的系数作一个小的扰动。我们希望 比较( 1.1)和( 1.2)根的差别，从而分析方程( 1.1)的解对扰动的敏感性。 实验内容 ：为了实现方便，我们先介绍两个 MATLAB 函数：“roots”和“ poly”。 u roots(a) 其中若变量 a 存储 n+1 维的向量，则该函数的输出 u 为一个 n 维的向量。设 a 的元素依次 为 a1,a2, ,an 1，则输出 u 的各分量是多项式方程 a1xn

3、a2xn 1an x an 1 0 的全部根；而函数 b poly (v 的输出 b 是一个 n+1 维向量，它是以 n 维向量 v 的各分量为根的多项式的系数。 可见“ roots” 和“ poly ”是两个互逆的运算函数。 ess 0.000000001; ve zeros(1,21); ve( 2) ess; roots( poly (1 : 20) ve) 上述简单的 MATLAB 程序便得到 （1.2）的全部根，程序中的“ess”即是（1.2）中的 实验要求 ： （1）选择充分小的 ess，反复进行上述实验，记录结果的变化并分析它们。如果扰动 项的系数 很小，我们自然感觉（ 1.1）

4、和（ 1.2）的解应当相差很小。计算中你 有什么出乎意料的发现？表明有些解关于如此的扰动敏感性如何？ （2）将方程（ 1.2）中的扰动项改成 x18 或其它形式，实验中又有怎样的现象出现？ （3）（选作部分）请从理论上分析产生这一问题的根源。注意我们可以将方程（1.2） 写成展开的形式， p（ x, ） x20 x190 （1.3） 同时将方程的解 x 看成是系数 的函数，考察方程的某个解关于 的扰动是否敏感， 与研究它关于 的导数的大小有何关系？为什么？你发现了什么现象，哪些根关于 的变 化更敏感？ 思考题一 ：（上述实验的改进） 在上述实验中我们会发现用 roots 函数求解多项式方程的精

5、度不高， 为此你可以考虑用 符号函数 solve 来提高解的精确度，这需要用到将多项式转换为符号多项式的函数 poly2sym, 函数的具体使用方法可参考 MATLAB 的帮助。 思考题二 ：（二进制产生的误差） 1000 用 MATLAB 计算 0.1 100 。结果居然有误差！因为从十进制数角度分析，这一计算 i1 应该是准确的。实验反映了 计算机内部的二进制本质 。 思考题三 ：（一个简单公式中产生巨大舍入误差的例子） 可以用下列式子计算自然对数的底数 这个极限表明随着 n的增加，计算 e值的精度是不确定的。现编程计算 f（n） （1 1）n与 n exp（1）值的差。 n大到什么程度的

6、时候误差最大？你能解释其中的原因吗？ 相关 MATLAB 函数提示： poly(a)求给定的根向量 a 生成其对应的多项式系数(降序)向量 roots(p)求解以向量 p 为系数的多项式(降序)的所有根 poly2sym(p) 将多项式向量 p 表示成为符号多项式(降序) sym(arg) 将数字、字符串或表达式 arg 转换为符号对象 syms arg1 arg2 argk 将字符 arg1,arg2,argk定义为基本符号对象 solve(eq1) 求符号多项式方程 eq1 的符号解 实验二 插值法 实验 2.1（多项式插值的振荡现象） 问题提出 ：考虑一个固定的区间上用插值逼近一个函数。

7、显然拉格朗日插值中使用的 节点越多，插值多项式的次数就越高。我们自然关心插值多项式的次数增加时，Ln（x） 是 否也更加靠近被逼近的函数。龙格（ Runge）给出一个例子是极著名并富有启发性的。设 区间-1,1 上函数 f ( x) 1 2 1 25x2 实验内容：考虑区间 -1,1的一个等距划分，分点为 xi1 2i , i 0,1,2, ,n n 则拉格朗日插值多项式为 n Ln(x) i0 1 2 1 2 5xi2 li x( 其中的 li（x）,i 0,1,2, ,n是 n 次拉格朗日插值基函数。 实验要求 ： （1） 选择不断增大的分点数目 n=2,3.，画出原函数 f（x）及插值多

8、项式函数 Ln（x） 在 -1，1上的图像，比较并分析实验结果。 （ 2）选择其他的函数，例如定义在区间 -5 ，5上的函数 h(x) 1 x4 , g ( x ) arctan x 重复上述的实验看其结果如何 3）区间 a,b上切比雪夫点的定义为 xk ,k 1,2, ,n 1 b a b acos(2k 1) 2 2 2(n 1) 以x1,x2, xn 1为插值节点构造上述各函数的拉格朗日插值多项式，比较其结果,试分析 原因。 实验 2.2（样条插值的收敛性） 问题提出 ：多项式插值是不收敛的，即插值的节点多，效果不一定就好。对样条函数 插值又如何呢？理论上证明样条插值的收敛性是比较困难的

9、，但通过本实验可以验证这一 理论结果 实验内容 ：请按一定的规则分别选择等距或者非等距的插值节点，并不断增加插值节点的个数。考虑实验 2.1 中的函数或选择其他你有兴趣的函数。 实验要求 ： （1）随节点个数增加，比较被逼近函数和样条插值函数误差的变化情况。分析所得结 果并与拉格朗日多项式插值比较 （可以用 MATLAB 的函数“ spline”作此函数的三次样条插 值，取 n=10、20，分别画出插值函数及原函数的图形 ）。 （2）样条插值的思想是早产生于工业部门。作为工业应用的例子考虑如下问题：某汽 车制造商用三次样条插值设计车门的曲线，其中一段的数据如下： xk 0 1 2 3 4 5

10、6 7 8 9 10 yk 0.0 0.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29 yk 0.8 0.2 要求：i.自己编程计算（用三弯矩、三转角方程均可） ii. 主函数 myspline（x,y,边界类型，边界值， xi ） 其中： x 节点 y 节点上的函数值 xi 未知节点 返回： S（xi） iii. 三对角方程组用追赶法求解 （书 P160） 实验 2.3：（ 一维插值的应用画图 ） 画你自己的手的形状，在 MATLAB中输入 figure（position,get（0,screensize） axes（position,0 0

11、1 1） x,y=ginput; 将你的手掌张开放在计算机屏幕上，然后使用计算机鼠标选取一系列点勾勒出手的轮 廓，按回车键结束 ginput过程，这样就获得了一系列你的手掌外形数据点 （x,y） 。也可以这 样获得数据点 （ x, y） ，先把手放在一张白纸上， 并用笔画出它的轮廓线， 然后将纸贴在计算 机屏幕上，透过纸能看到平面上的鼠标，并通过 ginput记录下轮廓上的点。 将x和 y坐标值看作是两个独立变量的函数，独立变量的取值为从1到记录的点的数 目。利用MATLAB 的插值函数进行插值，并画出你的手掌外形轮廓。 思考题 ：(二维插值) 1. 在一丘陵地带测量高程， x 和 y 方向每

12、隔 100 米测一个点，得高程数据 如下。试用 MATLAB 的二维插值函数“ interp2”进行插值，并由此找出最高点 和该点的高程。 y x 100 200 300 400 100 636 697 624 478 200 698 712 630 478 300 680 674 598 412 400 662 626 552 334 2. 画出山区地貌图。 利用 MATLAB 的 peaks 函数生成某山区的一些地点及其高度三维数据(单 位： m)。命令格式： x,y,z=peaks(n)，生成的 n 阶矩阵 x,y,z 为测量的山区地点 三维数据( n=30)。根据 peaks 函数生成

13、的数据，利用 Matlab 二维插值画出该 山区的地貌图和等值线图 (提示函数： interp2、 meshgrid、plot3 等)。 plot(x,y) subplot(m,n,k) yi=interp1(x,y,xi) pp=spline(x,y) 相关 MATLAB 函数提示： 作出以数据 (x(i),y(i) 为节点的折线图，其中 x,y 为同长度的 向量 将图形窗口分为 m*n 个子图，并指向第 k 幅图 根据数据 (x,y) 给出在 xi 的分段线性插值结果 yi 返回样条插值的分段多项式 (pp)形式结构 pp=csape(x,y, 边界类型 ,边界值 ) 生成各种边界条件的三

14、次样条插值 yi=ppval(pp,xi)pp 样条在 xi 的函数值 ZI=interp2(x,y,z,xi,yi) x,xi 为行向量， y,yi 为列向量， z 为矩阵的双线性二维插 值 ZI=interp2( ,spline) ZI=griddata(x,y,z,xi,yi) 使用二元三次样条插值 x,y,z 均为向量(不必单调)表示数据， xi,yi 为网格向量 的三角形线性插值(不规则数据的二维插值) 实验三 函数逼近与曲线拟合 实验 3.1(曲线逼近方法的比较) 问题提出 ：曲线的拟合和插值，是逼近函数的基本方法，每种方法具有各 自的特点和特定的适用范围，实际工作中合理选择方法是

15、重要的。 实验内容 ：考虑实验 2.1 中的著名问题。下面的 MATLAB 程序给出了该函 数的二次和三次拟合多项式。 x=-1:0.2:1; y=1/(1+25*x.*x); xx=-1:0.02:1; p2=polyfit(x,y,2); yy=polyval(p2,xx); plot(x,y, o ,xx,yy); xlabel( x ); ylabel( y ); hold on; p3=polyfit(x,y,3); yy=polyval(p3,xx); plot(x,y, o ,xx,yy); hold off; 适当修改上述 MATLAB 程序，也可以拟合其他你有兴趣的函数。 实

16、验要求 ： (1) 将拟合的结果与拉格朗日插值及样条插值的结果比较。 ( 2)归纳总结数值实验结果， 试定性地说明函数逼近各种方法的适用范围， 及实际应用中选择方法应注意的问题。 实验 3.2：(最小二乘拟合的经验公式和模型) 1.(已知经验公式) ：某类疾病发病率为 y 和年龄段 x (每五年为一段，例 bx 如 05 岁为第一段， 610 岁为第二段 )之间有形如 y ae 的经验关系，观 测得到的数据表如下 x 1 2 3 4 5 6 7 8 9 y 0.898 2.38 3.07 1.84 2.02 1.94 2.22 2.77 4.02 x 10 11 12 13 14 15 16

17、17 18 19 y 4.76 5.46 6.53 10.9 16.5 22.5 35.7 50.6 61.6 81.8 实验要求 ： bx ( 1)用最小二乘法确定模型 y ae 中的参数 a和b (提示函数： lsqcurvefit ， lsqnonlin )。 bx ( 2)利用 MATLAB 画出离散数据及拟合函数 y ae 图形。 3)利用 MATLAB 画出离散点处的误差图，并计算相应的均方误差。 2(最小二乘拟合模型未知) 某年美国轿车价格的调查资料如表，其中 xi 表示轿车的使用年数， yi 表示 相应的平均价格， 实验要求 ：试分析用什么形式的曲线来拟合表中的数据， 并预 测

18、使用 4.5 年后轿车的平均价格大致为多少？ xi 1 2 3 4 5 6 7 8 9 10 yi 2615 1943 1494 1087 765 538 484 290 226 204 实验 3.3(研究最佳平方逼近多项式的收敛性质) 实验内容和要求： 取函数 f(x) ex ，在-1 ，1 上以勒让德多项式为基函数， 对于 n 0,1, ,10构造最佳平方逼近多项式 pn(x)，令 n(x) f(x) pn(x) ，将 n(x) x的曲线画在一个图上。 令En max f(x) pn(x) ，画出 En n的曲线。做出 En n之间的最小二乘 1x1 曲线，能否提出关于收敛性的猜测。 思考

19、题一 ：(病态)考虑将0,130 等分节点，用多项式 y 1 xx5生成 数据，再用 polyfit 求其 3次、5次、10 次、 15次拟合多项式，并分析误差产生 的原因。 思考题二： 调研 Matlab 拟合工具箱的使用。 相关 MATLAB 函数提示： p=polyfit(x,y,k) 用 k 次多项式拟合向量数据 (x,y) ，返回多项式的降幂系数，当 k=n-1 时，实现多项式插值，这里 n 是向量维数 实验四 数值积分与数值微分 实验 4.1 (高斯数值积分方法用于积分方程求解) 问题提出 ：线性的积分方程的数值求解，可以被转化为线性代数方程组的 求解问题。 而线性代数方程组所含未

20、知数的个数， 与用来离散积分的数值方法的 节点个数相同。在节点数相同的前提下，高斯数值积分方法有较高的代数精度， 用它通常会得到较好的结果。 实验内容 ：求解第二类 Fredholm 积分方程 b y(t) k(t,s)y(s)ds f (t),a t b a 首先将积分区间 a,b等分成 n 份，在每个子区间上离散方程中的积分就得到 线性代数方程组。 实验要求 ：分别使用如下方法，离散积分方程中的积分 1.复化梯形方法；2.复化辛甫森方法； 3.复化高斯方法。求解如下的积分方程。 2 1 t t t 1) y(t) 2 0et y(s)ds et ，方程的准确解为 et ； e 1 0 11

21、 2) y(t) 0est y(s)ds 1 1e2t et ，方程的准确解为 et ； 比较各算法的计算量和误差以分析它们的优劣。 实验 4.2(高维积分数值计算的蒙特卡罗方法) 问题提出 ：高维空间中的积分，如果维数不很高且积分区域是规则的或者 能等价地写成多重积分的形式， 可以用一元函数积分的数值方法来计算高维空间 的积分。蒙特卡罗方法对计算复杂区域甚至不连通的区域上的积分并没有特殊的 困难。 实验内容 ：对于一般的区域 ，计算其测度(只要理解为平面上的面积或 空间中的体积)的一般方法是：先找一个规则的区域 A 包含 ，且 A 的测度是 已知的。生成区域 A 中 m 个均匀分布的随机点

22、pi ,i 1,2, ,m ，如果其中有 n 个 落在区域 中，则区域 的测度 m( )为 n/m。函数 f(x) 在区域 上的积分可以 近似为：区域 的测度与函数 f(x)在 中 n个随机点上平均值的乘积，即 1 f(x)dx m( ) 1f (pk ) n pk 实验要求 ：假设冰琪淋的下部为一锥体而上面为一半球，考虑冰琪淋体积问 题：计算锥面 z2 x2 y2 上方和球面 x2 y2 (z 1)2 1内部区域的体积。如果 使用球面坐标，该区域可以表示为如下的积分： 2cos 2 2 sin d d d 用蒙特卡罗方法可以计算该积分。 另一方面，显然这样的冰琪淋可以装在如下立方体的盒子里

23、1 x 1, 1 y 1,0 z 2 而该立方体的体积为 8。只要生成这个盒子里均匀分布的随机点，落入冰琪淋锥 点的个数与总点数之比再乘以 8 就是冰琪淋锥的体积。 比较两种方法所得到的结 果。 类似的办法可以计算复杂区域的测度(面积或体积) 。试求由下列关系所界 定区域的测度： 0 x 1,1 y 2, 1 z 3 (1) ex y sin( z) y 0 1 x 3, 1 y 4 33 (2) x3 y3 29 y ex 2 0 x 1,0 y 1,0 z 1 (3) x2 sin y z x z ex 1 相关 MATLAB 函数提示： diff(x)如果 x 是向量，返回向量 x 的差

24、分；如果 x 是矩阵，则按各列作 差分 diff(x,k)k 阶差分 q=polyder(p) 求得由向量 p 表示的多项式导函数的向量表示 q Fx=gradient(F,x) 返回向量 F 表示的一元函数沿 x 方向的导函数 F(x) ，其中 x 是与 F 同维数的向量 z=trapz(x,y) x 表示积分区间的离散化向量； y 是与 x 同维数的向量， 表示被积 函数； z 返回积分的近似值 z=guad(fun,a,b,tol) 自适应步长 Simpson 积分法求得 Fun 在区间 a,b 上的定积 分， Fun为M 文件函数句柄， tol 为积分精度 z=dblquad(fun,

25、a,b,c,d,tol,method) 求得二元函数 Fun(x,y)的重积分 z=triplequad(fun,a,b,c,d,e,f,tol,method) 求得三元函数 Fun(x,y,z)的重积分 10 实验五 解线性方程组的直接方法 实验 5.1 （主元的选取与算法的稳定性） 问题提出 ：Gauss消去法是我们在线性代数中已经熟悉的。 但由于计算机的 数值运算是在一个有限的浮点数集合上进行的， 如何才能确保 Gauss消去法作为 数值算法的稳定性呢？ Gauss消去法从理论算法到数值算法，其关键是主元的选 择。主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。 实验内

26、容 ：考虑线性方程组 Ax b, A Rn n,b Rn 编制一个能自动选取主元，又能手动选取主元的求解线性方程组的 Gauss 消去过程。 实验要求 ： 61 861 （ 1）取矩阵 A 861 86 7 15 b ，则方程有解 x* (1,1, ,1)T 15 14 取 n=10 计算矩阵的条件数。让程序自动选取主元 （顺序消元 ） ，结果如何？ （2）现选择程序中手动选取主元的功能。每步消去过程总选取按模最小或 按模尽可能小的元素作为主元， 观察并记录计算结果。 若每步消去过程总选取按 模最大的元素作为主元，结果又如何？分析实验的结果。 （ 3）取矩阵阶数 n=20 或者更大，重复上述实

27、验过程，观察记录并分析不 同的问题及消去过程中选择不同的主元时计算结果的差异， 说明主元素的选取在 消去过程中的作用。 （4）选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。重复上 述实验，观察记录并分析实验结果。 实验 5.2（线性代数方程组的性态与条件数的估计） 问题提出 ：理论上，线性代数方程组 Ax b 的摄动满足 x co n（dA） A b x 1 A 1 A A b 矩阵的条件数确实是对矩阵病态性的刻画，但在实际应用中直接计算它显 然不现实，因为计算 A 1 通常要比求解方程 Ax b 还困难。 实验内容 ：MATLAB 中提供有函数 “condest”可以用来估计矩阵的条

28、件数， 它给出的是按 1-范数的条件数。首先构造非奇异矩阵 A 和右端，使得方程是可 以精确求解的。再人为地引进系数矩阵和右端的摄动 A和 b，使得 A 和 b 充 11 分小 实验要求： 1）假设方程 Ax=b 的解为 x，求解方程 （A A）x? b b ，以 1-范数， 给出 x x? x 的计算结果。 给出 x x 的计算结果。 （2）选择一系列维数递增的矩阵 （可以是随机生成的） ，比较函数“condest” 所需机器时间的差别 .考虑若干逆是已知的矩阵，借助函数“eig”很容易给出 cond2（A） 的数值。将它与函数“ cond（A,2）”所得到的结果进行比较。 （3）利用“co

29、ndest”给出矩阵 A 条件数的估计， 针对（ 1）中的结果给出 x x 的理论估计，并将它与（ 1）给出的计算结果进行比较，分析所得结果。注意， 如果给出了 cond（A） 和 A 的估计，马上就可以给出 A 1 的估计。 （ 4）估计著名的 Hilbert 矩阵的条件数。 1 H (hi,j)n n, hi,j i j 1, i,j 1,2, ,n 思考题一：（ Vadermonde矩阵） 1 x0 2 x0 x0n 1 2 n x1 x1 x1 1 2 n x2 x2 x2 2 n 1 xn xn xn A ，b n xi x0 i0 n x1i i0 n x2i i0 n xni i

30、0 其中， xk 1 0.1k,k 0,1, ,n ， （1）对 n=2,5,8，计算 A 的条件数；随 n 增大，矩阵性态如何变化？ （2）对 n=5，解方程组 Ax=b ；设 A 的最后一个元素有扰动 10-4，再求解 Ax=b （ 3）计算（ 2）扰动相对误差与解的相对偏差，分析它们与条件数的关系。 （4）你能由此解释为什么不用插值函数存在定理直接求插值函数而要用拉格朗 日或牛顿插值法的原因吗？ 12 相关 MATLAB 函数提示： zeros(m,n) 生成 m 行， n 列的零矩阵 ones(m,n) 生成 m 行，n 列的元素全为 1 的矩阵 eye(n) 生成 n 阶单位矩阵 r

31、and(m,n) 生成 m 行,n 列(0,1)上均匀分布的随机矩阵 diag(x) 返回由向量 x 的元素构成的对角矩阵 tril(A) 提取矩阵 A 的下三角部分生成下三角矩阵 triu(A) 提取矩阵 A 的上三角部分生成上三角矩阵 rank(A) 返回矩阵 A 的秩 det(A) 返回方阵 A 的行列式 inv(A) 返回可逆方阵 A 的逆矩阵 V ，D=eig(A) 返回方阵 A 的特征值和特征向量 norm(A,p) 矩阵或向量的 p 范数 cond(A,p) 矩阵的条件数 L ， U，P=lu(A) 选列主元 LU 分解 R=chol(X) 平方根分解 Hi=hilb(n) 生成

32、 n 阶 Hilbert 矩阵 13 实验六 解线性方程组的迭代法 实验 6.1（病态的线性方程组的求解） 问题提出 ：理论的分析表明，求解病态的线性方程组是困难的。实际情况 是否如此，会出现怎样的现象呢？ 实验内容 ：考虑方程组 Hx=b 的求解，其中系数矩阵 H 为 Hilbert 矩阵， 1 H （hi,j ）n n, hi ,j, i, j 1,2, ,n i j 1 这是一个著名的病态问题。 通过首先给定解 （例如取为各个分量均为 1）再 计算出右端 b 的办法给出确定的问题。 实验要求 ： （1）选择问题的维数为 6，分别用 Gauss消去法、列主元 Gauss消去法、 J 迭代法

33、、GS迭代法和 SOR 迭代法求解方程组， 其各自的结果如何？将计算结果 与问题的解比较，结论如何？ （ 2）逐步增大问题的维数（至少到 100），仍然用上述的方法来解它们，计 算的结果如何？计算的结果说明了什么？ （3）讨论病态问题求解的算法 实验 6.2 书上 P211 计算实习题 2，其中 N=100，第二小问改为用 Jacobi迭 代、G-S 迭代、红黑排序的 G-S 迭代求解，并比较他们之间的收敛速度； 进一步， 用 BSOR 迭代求解，试找出最优松弛因子。 14 实验七 非线性方程求根 实验 7.1（迭代法、初始值与收敛性） 实验目的 ：初步认识非线性问题的迭代法与线性问题迭代法的

34、差别，探讨 迭代法及初始值与迭代收敛性的关系。 问题提出 ：迭代法是求解非线性方程的基本思想方法，与线性方程的情况 一样，其构造方法可以有多种多样， 但关键是怎样才能使迭代收敛且有较快的收 敛速度。 实验内容 ：考虑一个简单的代数方程 x2 x 1 0 针对上述方程，可以构造多种迭代法，如 xn 1 xn2 1 （7.1） 1 xn 1 1 （7.2） xn xn 1 xn 1 （7.3） 在实轴上取初始值 x0，请分别用迭代（ 7.1）-（7.3）作实验，记录各算法 的迭代过程。 实验要求 ： （ 1）取定某个初始值，分别计算（ 7.1） -（7.3）迭代结果，它们的收敛性 如何？重复选取不

35、同的初始值， 反复实验。 请自选设计一种比较形象的记录方式 （如利用 MATLAB 的图形功能），分析三种迭代法的收敛性与初值选取的关系。 （2）对三个迭代法中的某个，取不同的初始值进行迭代，结果如何？试分 析迭代法对不同的初值是否有差异？ （3）线性方程组迭代法的收敛性是不依赖初始值选取的。比较线性与非线 性问题迭代的差异，有何结论和问题。 相关 MATLAB 函数提示： x=fzero（fun,x0） 返回一元函数 fun 的一个零点，其中 fun 为函数句柄， x0 为标 量时，返回在 x0 附近的零点；x0 为向量 a,b时，返回函数在 a,b 中的零点 x,f,h=fsolve（fu

36、n,x0） 返回一元或多元函数 x0 附近 fun 的一个零点， 其中 fun 为 函数句柄， x0 为迭代初值； f 返回 fun 在 x 的函数值，应该接近 0； h 返回值如果大于 0，说明计算结果可靠，否则不可靠 15 实验八 常微分方程初值问题数值解法 实验 8.1（Lorenz 问题与混沌） 问题提出 ：考虑著名的 Lorenz 方程 (8.1) ddxt s(y x) dy rx y xz ddtz dz xy bz dt 其中 s,r,b 为变化区域有一定限制的实参数。该方程形式简单，表面上看并 无惊人之处，但由该方程提示出的许多现象，促使“混沌”成为数学研究的崭新 领域，在实

37、际应用中也产生了巨大的影响。 实验内容：先取定初值 y0=（0,0,0）,参数 s=10,r=28,b=8/3,用MATLAB 的数值 求常微分方程函数 ods45编程对（ 8.1）进行求解 实验要求 ： （1）对目前取定的参数值 s,r 和 b，选取不同的初值 y0 进行运算，观察计 算的结果有什么特点？解的曲线是否有界？解的曲线是不是周期的或趋于某个 固定点？ （2）在问题允许的范围内适当改变其中的参数值 s,r,b，再选取不同的初始 值 y0 进行运算，观察并记录计算的结果有什么特点？是否发现什么不同的现 象？ 思考题一 ：考虑如下常微分方程 x rx bsi na（x） k0 k1c

38、o sm（t） 先固定其中参数： r=0.9,b=3,a=1,k0=0,m=1,初始值为原点， t 的范围为 0-500 分别取不同的 k1,如 1，1.5，2，2.5，3，3.5 等，绘图观察相位图（ x,x）的 变化情况，能找到其中的规律吗？ 试试改变其它参数呢？ 16 相关 MATLAB 函数提示： t,y=ode45(odefun,tspan,y0) odefun 表示 f(t,y) 的函数句柄， t 是标量， y 是标量 或向量； tspan 是二维向量 t0,tf ，表示自变量初值 t0 和终值 tf ；y0 表示初值向量，若无输出参数，则作出图形 17 附录 MATLAB 简介

39、这里介绍 MATLAB 一些入门知识，包括 MATLAB 桌面和窗口， MATLAB 命令格式、数据格式、数据文件和变量管理， MATLAB 的数组和矩阵运算， MATLAB 的字符串、元胞和结构等数据类型， MATLAB 的程序设计方法， MATLAB 作图方法在线帮助的使用和程序文件和目录的管理等。 表一 MATLAB 的基本命令 主题词 含义 主题词 含义 format 设置数据显示格式 feval 函数求值 who 显示变量名 input 提示输入 whos 显示变量信息 disp 输出 clear 清除内存变量 tic 启动秒表 save 保存工作变量到文件 toc 时间读数(秒)

40、load 从文件装载变量 help 帮助 linspace 区间等分 lookfor 查找 length 获取数组长度 type 列程序清单 size 矩阵大小 which 查找文件目录 max 最大值 double 双精度 min 最小值 str2num 字符串转化为数值 sum 求和 num2str 数值转化为字符串 find 条件检索 MATLAB 桌面 启动 MATLAB 后，就进入 MATLAB 的桌面，图 1 为 MATLAB6.1 的默认 ( Default)桌面。第一行为菜单栏， 第二行为工具栏， 下面是三个最常用的窗口。 右边最大的是命令窗口( Command Window

41、)，左上方前台为发行说明书窗口 ( Launch pad)，后台为工作空间( Workspace)，左下方为命令历史( Command History )后台为当前目录( Current Directory)。 1. 窗口 (1)命令窗口 该窗口是进行 MATLAB 操作最主要的窗口。窗口中“ ”为命令输入提示 符，其后输入运算命令，按回车键就可执行运算，并显示运算结果 .。 18 图1 （2）发行说明书窗口 发行说明书窗口是 MATLAB 所特有的，用来说明用户所拥有的 Mathworks 公司产品的工具包、演示以及帮助信息。 （3）工作空间 在默认桌面，位于左上方窗口前台，列出内存中 MA

42、TLAB 工作空间的所有 变量的变量名、 尺寸、字节数。 用鼠标选中变量， 击右键可以打开、 保存、删除、 绘图等操作。 （4）当前目录 在默认桌面， 位于左下方窗口后台， 用鼠标点击可以切换到前台。 该窗口列 出当前目录的程序文件（ .m）和数据文件（ .mat）等。用鼠标选中文件，击右键 可以进行打开、运行、删除等操作。 （5）命令历史（ Command History） 该窗口列出在命令窗口执行过的 MATLAB 命令行的历史记录。用鼠标选中 命令行，击右键可以进行复制、执行（ Evaluate Selection）、删除等操作。 除上述窗口外， MATLAB 常用窗口还有编程器窗口、图

43、形窗口等。 、数据和变量 1.表达式 在命令窗口作一些简单的计算， 就如同使用一个功能强大的计算器， 使用变 量无须预先定义类型。 19 例如，设球半径为 r=2，求球的体积 V 4 r3 。 3 r=2%表达式将 2 赋予变量 r r=%系统返回 r 的值 2 v=4/3*pi*r3%pi 为内置常量 ，乘方用 表示 v= 33.5103 几个表达式可以写在一行，用分号（； ）或逗号（，）分割，用分号（；）使 该表达式运算结果不显示，而逗号（， ）则显示结果。也可以将一个长表达式分 在几行上写，用三点（ ）续行。 若需要修改已执行过的命令行， 可以在命令历史中找到该命令行复制， 再粘 贴至命

44、令窗口修改。也可以直接使用键盘 调出已执行过的命令行修改。 2. 数据显示格式 MATLAB 默认的数据显示格式为短格式（ short）：当结果为整数，就作为整 数显示；当结果是实数， 以小数点后四位的长度显示。 若结果的有效数字超出一 定范围，以科学计数法显示（如 3.2000e-006表示3.2 10 6 ）。数据显示格式可使 用命令 Format 改变。例如： format long;v%长格式， 16 位 v = 33.51032163829112 format short;v%短格式 v = 33.5103 format rational;v %有理格式，近似分数 v = 6501/

45、194 3. 复数 MATLAB 中复数可以如同实数一样，直接输入和计算。例如： a=1+2i;b=5-4*i;c=a/b c = -0.0732 + 0.3415i 4. 预定义变量 MATLAB 有一些预定义变量（表 1），启动时就已赋值，可以直接使用，如 前我们使用的圆周率 pi 和虚数单位 i. 表1 常用预定义变量 变量名 说明 i或j 虚数单位 1 pi 圆周率 3.14159 eps 浮点数识别精度 2(-52)= 2.2204 10 16 realmin 最小正实数 2.2251 10 308 realmax 最大正实数 1.7977 10308 inf 无穷大 20 NaN

46、没有意义的数 预定义变量在工作空间观察不到。 如果预定义变量被用户重新赋值， 则原来 的功能暂不能使用。当这些用户变量被清除( clear)或 MATLAB 重新启动后， 这些功能得以恢复。 5. 用户变量 MATLAB 变量名总以字母开头，以字母、数字或下划线组成，区分大小写， 有效字符长度为 63 个。如 A,a,a1,a_b都是合法的，且 a 与 A 表示不同变量。在 Command Window 中使用的变量一旦被赋值，就会携带这个值存在于工作空间， 直到被清除或被赋予新的值。 ans是系统一个特别的变量名。若一个表达式运算结果没有赋予任何变量， 系统自动用 ans存放答案。例如： A

47、=5+4i;b=5-4*i;B=1;A*b%没有定义 A*b 的输出变量 ans = 41%ans 来接受计算结果，注意这是大写 A 与 小写 b 的乘积，尽管我们可以使用工作空间来查询和清除变量， 但使用下列命令 方式更快捷： whos%查询 Workspace中的变量列表 Name Size Bytes Class A 1x1 16 double array (complex) B 1x1 8 double array a 1x1 16 double array (complex) ans 1x1 8 double array b 1x1 16 double array (complex)

48、 c 1x1 16 double array (complex) Grand total is 6 elements using 80 bytes A%查询变量 A 的值 A = 5.0000 + 4.0000i clear A%清除变量 A A%再查询 A 的值，已经不存在了 ? Undefined function or variable A. clear%清除 Workspace中所有变量 whos%Workspace中已没有任何变量了 三、数组和矩阵运算 MATLAB 基本数据单元是无需指定维数的数组。 数组运算是 MATLAB 最鲜 明的特点，一方面可以使得计算程序简明易读，另一方面

49、可以提高计算速度。 1. 数组的输入 最常用的数组是双精度数值数组( double array)。一维数组相当于向量，二 维数组相当于矩阵， 一维数组可以视为二维数组的特例。 二维数组的第一维称为 “行”，第二维称为“列” 。MATLAB 数组无需预先定义维数。直接输入数组的 21元素，用中括号( )表示一个数组，同行元素间用空格或逗号分隔，不同行间 用分号或回车分隔，例如： linspace生 clear;a=1,2,3;4,5,6;7,8,9 a = 1 2 3 4 5 6 7 8 9 或 a=1 2 3 %这种方式特别适用于大型矩阵 4 5 6 7 8 9 a = 1 2 3 4 5 6

50、 7 8 9 对于等差数列构造的一维数组，可用冒号运算生成，也可用函数 成。 b=0:3:10 %初值： 增量：终值 b = 0 3 6 9 b=0:10 %增量为 1 可省略 b = 0 1 2 3 4 5 6 7 10 b=10:-3:0 b = 10 7 4 1 b=linspace(0,10,4) b = 0 3.3333 6.6667 length(b) ans = 4 b(3) ans = 6.6667 b(1,end) %递减 %将区间0， 10等分为 4-1=3 份 10.0000 %查询 b 的长度 %查询 b 的第三个元素 %查询 b 的首和尾元素 ans = 0 10 二

51、维数组元素双下标编址按通常方式，单下标编址按列排序。 size(a)%查询数组 a 的尺寸 ans = 3 3 a(3,2),a(6) 22 ans = c=a(1 3,2 3) 块矩阵） c = 2 3 8 9 d=a(2,:) d = 4 5 6 a(:) ans = 1 4 ans = 8 %提取 a 的第一、第三行和第二、第三列（分 %提取 a 的第二行 %将 a 所有元素按单下标顺序排为列向量 7 2 5 8 3 6 9 一些特殊的二维数组可以用函数产生，例如： a=zeros(2,4) %生成 2行 4列零矩阵 a = 0 0 0 0 0 0 0 0 b=ones(1,4) %生成

52、 1行 4列1矩阵 b = 1 1 1 1 c=a;b %拼接 c = 0 0 0 0 0 0 0 0 1 1 1 1 c(2,1)=100 %修改部分元素 c = 0 0 0 0 100 0 0 0 1 1 1 1 reshape(c,2,6) %按 2 行 6 列重排矩阵元素 ans = 0 1 0 0 1 0 23 100 0 1 0 0 1 注意：数组下标对应矩阵的行和列，编址一律从 1 开始，不能用 0. 矩阵输入也可用“ load”命令从外部数据文件导入 2. 数组运算 数组运算是指数组对应元素之间的运算， 也称点运算。 矩阵的乘法、 乘方和 除法有特殊的数学含义， 并不是数组对应

53、元素的运算， 所数组乘法、 乘方和除法 的运算符前特别加了一个点。 特别要区分数组运算在乘法、 乘方和除法上的意义 和表示上与矩阵运算的不同。 表2 数组运算符 运算 符号 说明 数组加与减 A+B 与 A-B 对应元素之间加减 数乘数组 k*A 或 A*k k乘 A的每个元素 数与数组加减 k+A 或 k-A k 加（减） A 的每个元素 数组乘数组 A.*B 点运算只有点乘、点乘方、点除三个，表 示对应元素之间的运算； （.* ）是一个整 体，点（.）不能漏掉，（.）和（ *）之间也 不能有空格 数组乘方 A.k,k.A 数除以数组 k./A 数组除法 左除 A.B,右除 B./A cle

54、ar;A=1 -1;0 2;B=0 1;1 -1; A.*B%注意不是 A*B ans = 0 -1 0 -2 A.B,A./B Warning: Divide by zero. ans = 0 -1.0000 Inf -0.5000 Warning: Divide by zero. ans = Inf -1 0 -2 A.2 ans = 1 1 0 4 1./A Warning: Divide by zero. ans = 1.0000 -1.0000 Inf 0.5000 3. 矩阵运算 矩阵是一个二维数组， 所以矩阵的加、 减、数乘等运算与数组运算是一致的 但是有两点需要注意： 24 （

55、1）对于乘法、乘方和除法等三种运算，矩阵运算与数组运算的运算符及 含义不同：矩阵运算按线性变换定义， 使用通常符号； 数组运算按对应元素运算 定义，使用点运算符； （2）数与矩阵加减、矩阵除法在数学上是没有意义的，在MATLAB 中为 简便起见，定义了这两类运算，其含义见表 3. 表3 矩阵运算符 运算 符号 说明 转置 A 加与减 A+B 与 A-B 同数组运算 数乘矩阵 k*A 或 A*k 同数组运算 矩阵乘法 A*B 矩阵乘方 Ak 数与矩阵加减 k+A 与 k-A k+A 等价于 k*ones(size(A)+A 矩阵除法 左除 AB ，右除 B/A 它们分别为矩阵方程 AX=B 和

56、XA=B 的解 A=1 2;3 4;B=4 3;2 1; 100+A ans = 101 102 103 104 A*B,A.*B%注意矩阵运算和数组运算的区别 ans = 8 5 20 13 ans = 4 6 6 4 AB,B/A,A.B,B./A%注意矩阵运算和数组运算的区别 ans = -6.0000 -5.0000 5.0000 4.0000 ans = -3.5000 2.5000 -2.5000 1.5000 ans = 4.0000 1.5000 0.6667 0.2500 ans = 4.0000 1.5000 0.6667 0.2500 4.数学函数 数组的数学函数也是按每

57、个元素的运算， 使用通常的函数符号， 常用数学函 数见表 4 25 表 4 数学函数 函数 意义 函数 意义 sin 正弦 fix 向 0 取整 cos 余弦 mod 模余 tan 正切 rem 除法余数 cot 余切 abs 绝对值(模) asin 反正弦 exp 指数函数 acos 反余弦 log 自然对数 sqrt 开方 log10 以 10 为底的对数 A=4 -1;3 2; B=exp(A) B = 54.5982 0.3679 20.0855 7.3891 C=fix(B) C = 54 0 20 7 D=sin(C) D = -0.5588 0 0.9129 0.6570 E=l

58、og(D) Warning: Log of zero. E = -0.5820 + 3.1416i -Inf -0.0911 -0.4201 5.关系与逻辑运算 MATLAB 的关系运算和逻辑运算符都是对于元素的操作，其结果是特殊的 逻辑数组( logical array)表 5，“真”用 1 表示，“假”用 0 表示，而逻辑运算中， 所有非零元素作为 1(真)处理。 表5 关系运算和逻辑运算 运算符 含义 运算符 含义 小于 b=sqrt(a); a=,num2str(a),a 的开方 =,num2str(b) ans = a=12,a的开方 =3.4641 MATLAB 命令可以定义成一个

59、字符串，使用 eval 可以使该字符串所表达的 MATLAB 命令得到执行。 fun=x.2.*sin(x); x=1;eval(fun) ans = 0.8415 x=1:3;eval(fun) ans = 0.8415 3.6372 1.2701 2. 元胞和结构 不管是数值数组还是字符数组， 其数据结构必须是整齐的。 首先数值和字符 不能混合，其次小数组拼接成大数组时，其尺寸 (size)必须相符 (agree)。例如： A=first;second%错误 ? Error using = vertcat All rows in the bracketed expression must

60、have the same number of columns. 将不同类型、不同尺寸的数组，加大括号( )，可构成一个元胞。几个元 胞可以构成元胞数组。 Ac1=first;1:3;Ac2=second;1 2;3 4; Ac=Ac1,Ac2 Ac = first second 1x3 double2x2 double size(Ac) ans = 2 2 Ac(2,1)%小括号，查询 Ac 的第二行第一列元素 ans = 1x3 double 28 Ac2,1 内容 %大括号，查询 Ac 的第二行第一列元素的具体 ans = 1 2 3 一个结构通过“域”来定义，比元胞更丰富、更灵活。几个