关 键 词：

卫国 MATLAB 程序设计 应用

资源描述：

刘卫国-MATLAB程序设计与应用（第2版）,卫国,MATLAB,程序设计,应用

内容简介：

第6章MATLAB数值计算6 1数据处理与多项式计算6 2数值微积分6 3离散傅立叶变换6 4线性方程组求解6 5非线性方程与最优化问题求解6 6常微分方程的数值求解6 7稀疏矩阵 6 1数据处理与多项式计算6 1 1数据统计与分析1 求矩阵最大元素和最小元素MATLAB提供的求数据序列的最大值和最小值的函数分别为max和min 两个函数的调用格式和操作过程类似 1 求向量的最大值和最小值y max X 返回向量X的最大值存入y 如果X中包含复数元素 则按模取最大值 y I max X 返回向量X的最大值存入y 最大值的序号存入I 如果X中包含复数元素 则按模取最大值 求向量X的最小值的函数是min X 用法和max X 完全相同 例求向量x的最大值 命令如下 x 43 72 9 16 23 47 y max x 求向量x中的最大值 y l max x 求向量x中的最大值及其该元素的位置 2 求矩阵的最大值和最小值求矩阵A的最大值的函数有3种调用格式 分别是 max A 返回一个行向量 向量的第i个元素是矩阵A的第i列上的最大值 Y U max A 返回行向量Y和U Y向量记录A的每列的最大值 U向量记录每列最大值的行号 max A dim dim取1或2 dim取1时 该函数和max A 完全相同 dim取2时 该函数返回一个列向量 其第i个元素是A矩阵的第i行上的最大值 求最小值的函数是min 其用法和max完全相同 例6 1分别矩阵A中各列和各行元素中的最大值 并求整个矩阵的最大值和最小值 3 两个向量或矩阵对应元素的比较函数max和min还能对两个同型的向量或矩阵进行比较 调用格式为 U max A B A B是两个同型的向量或矩阵 结果U是与A B同型的向量或矩阵 U的每个元素等于A B对应元素的较大者 U max A n n是一个标量 结果U是与A同型的向量或矩阵 U的每个元素等于A对应元素和n中的较大者 min函数的用法和max完全相同 例求两个2 3矩阵x y所有同一位置上的较大元素构成的新矩阵p 2 求矩阵的平均值和中值求数据序列平均值的函数是mean 求数据序列中值的函数是median 两个函数的调用格式为 mean X 返回向量X的算术平均值 median X 返回向量X的中值 mean A 返回一个行向量 其第i个元素是A的第i列的算术平均值 median A 返回一个行向量 其第i个元素是A的第i列的中值 mean A dim 当dim为1时 该函数等同于mean A 当dim为2时 返回一个列向量 其第i个元素是A的第i行的算术平均值 median A dim 当dim为1时 该函数等同于median A 当dim为2时 返回一个列向量 其第i个元素是A的第i行的中值 3 矩阵元素求和与求积数据序列求和与求积的函数是sum和prod 其使用方法类似 设X是一个向量 A是一个矩阵 函数的调用格式为 sum X 返回向量X各元素的和 prod X 返回向量X各元素的乘积 sum A 返回一个行向量 其第i个元素是A的第i列的元素和 prod A 返回一个行向量 其第i个元素是A的第i列的元素乘积 sum A dim 当dim为1时 该函数等同于sum A 当dim为2时 返回一个列向量 其第i个元素是A的第i行的各元素之和 prod A dim 当dim为1时 该函数等同于prod A 当dim为2时 返回一个列向量 其第i个元素是A的第i行的各元素乘积 例6 2求矩阵A的每行元素的乘积和全部元素的乘积 4 矩阵元素累加和与累乘积在MATLAB中 使用cumsum和cumprod函数能方便地求得向量和矩阵元素的累加和与累乘积向量 函数的调用格式为 cumsum X 返回向量X累加和向量 cumprod X 返回向量X累乘积向量 cumsum A 返回一个矩阵 其第i列是A的第i列的累加和向量 cumprod A 返回一个矩阵 其第i列是A的第i列的累乘积向量 cumsum A dim 当dim为1时 该函数等同于cumsum A 当dim为2时 返回一个矩阵 其第i行是A的第i行的累加和向量 cumprod A dim 当dim为1时 该函数等同于cumprod A 当dim为2时 返回一个向量 其第i行是A的第i行的累乘积向量 5 求标准方差在MATLAB中 提供了计算数据序列的标准方差的函数std 对于向量X std X 返回一个标准方差 对于矩阵A std A 返回一个行向量 它的各个元素便是矩阵A各列或各行的标准方差 std函数的一般调用格式为 Y std A flag dim 其中dim取1或2 当dim 1时 求各列元素的标准方差 当dim 2时 则求各行元素的标准方差 flag取0或1 当flag 0时 按 1所列公式计算标准方差 当flag 1时 按 2所列公式计算标准方差 缺省flag 0 dim 1 例6 4对二维矩阵x 从不同维方向求出其标准方差 6 相关系数MATLAB提供了corrcoef函数 可以求出数据的相关系数矩阵 corrcoef函数的调用格式为 corrcoef X 返回从矩阵X形成的一个相关系数矩阵 此相关系数矩阵的大小与矩阵X一样 它把矩阵X的每列作为一个变量 然后求它们的相关系数 corrcoef X Y 在这里 X Y是向量 它们与corrcoef X Y 的作用一样 例6 5生成满足正态分布的10000 5随机矩阵 然后求各列元素的均值和标准方差 再求这5列随机数据的相关系数矩阵 命令如下 X randn 10000 5 M mean X D std X R corrcoef X 7 排序MATLAB中对向量X是排序函数是sort X 函数返回一个对X中的元素按升序排列的新向量 sort函数也可以对矩阵A的各列或各行重新排序 其调用格式为 Y I sort A dim 其中dim指明对A的列还是行进行排序 若dim 1 则按列排 若dim 2 则按行排 Y是排序后的矩阵 而I记录Y中的元素在A中位置 6 1 2数据插值1 一维数据插值在MATLAB中 实现这些插值的函数是interp1 其调用格式为 Y1 interp1 X Y X1 method 函数根据X Y的值 计算函数在X1处的值 X Y是两个等长的已知向量 分别描述采样点和样本值 X1是一个向量或标量 描述欲插值的点 Y1是一个与X1等长的插值结果 method是插值方法 允许的取值有 linear nearest cubic spline 注意 X1的取值范围不能超出X的给定范围 否则 会给出 NaN 错误 例6 7给出概率积分的数据表如表6 1所示 用不同的插值方法计算f 0 472 例6 8某检测参数f随时间t的采样结果如表5 1 用数据插值法计算t 2 7 12 17 22 17 32 37 42 47 52 57时的f值 2 二维数据插值在MATLAB中 提供了解决二维插值问题的函数interp2 其调用格式为 Z1 interp2 X Y Z X1 Y1 method 其中X Y是两个向量 分别描述两个参数的采样点 Z是与参数采样点对应的函数值 X1 Y1是两个向量或标量 描述欲插值的点 Z1是根据相应的插值方法得到的插值结果 method的取值与一维插值函数相同 X Y Z也可以是矩阵形式 同样 X1 Y1的取值范围不能超出X Y的给定范围 否则 会给出 NaN 错误 例6 9设z x2 y2 对z函数在 0 1 0 2 区域内进行插值 例6 10某实验对一根长10米的钢轨进行热源的温度传播测试 用x表示测量点0 2 5 10 米 用h表示测量时间0 30 60 秒 用T表示测试所得各点的温度 试用线性插值求出在一分钟内每隔10秒 钢轨每隔0 5米处的温度 6 1 3曲线拟合在MATLAB中 用polyfit函数来求得最小二乘拟合多项式的系数 再用polyval函数按所得的多项式计算所给出的点上的函数近似值 polyfit函数的调用格式为 P S polyfit X Y m 函数根据采样点X和采样点函数值Y 产生一个m次多项式P及其在采样点的误差向量S 其中X Y是两个等长的向量 P是一个长度为m 1的向量 P的元素为多项式系数 polyval函数的功能是按多项式的系数计算x点多项式的值 例6 11用一个3次多项式在区间 0 2 内逼近函数 命令如下 X linspace 0 2 pi 50 Y sin X P polyfit X Y 3 得到3次多项式的系数和误差 6 1 4多项式计算1 多项式的四则运算 1 多项式的加减运算 2 多项式乘法运算函数conv P1 P2 用于求多项式P1和P2的乘积 这里 P1 P2是两个多项式系数向量 3 多项式除法函数 Q r deconv P1 P2 用于对多项式P1和P2作除法运算 其中Q返回多项式P1除以P2的商式 r返回P1除以P2的余式 这里 Q和r仍是多项式系数向量 deconv是conv的逆函数 即有P1 conv P2 Q r 2 多项式的导函数对多项式求导数的函数是 p polyder P 求多项式P的导函数p polyder P Q 求P Q的导函数 p q polyder P Q 求P Q的导函数 导函数的分子存入p 分母存入q 上述函数中 参数P Q是多项式的向量表示 结果p q也是多项式的向量表示 3 多项式求值MATLAB提供了两种求多项式值的函数 polyval与polyvalm 它们的输入参数均为多项式系数向量P和自变量x 两者的区别在于前者是代数多项式求值 而后者是矩阵多项式求值 1 代数多项式求值polyval函数用来求代数多项式的值 其调用格式为 Y polyval P x 若x为一数值 则求多项式在该点的值 若x为向量或矩阵 则对向量或矩阵中的每个元素求其多项式的值 例6 14已知多项式x4 8x3 10 分别取x 1 2和一个2 3矩阵为自变量计算该多项式的值 2 矩阵多项式求值polyvalm函数用来求矩阵多项式的值 其调用格式与polyval相同 但含义不同 polyvalm函数要求x为方阵 它以方阵为自变量求多项式的值 设A为方阵 P代表多项式x3 5x2 8 那么polyvalm P A 的含义是 A A A 5 A A 8 eye size A 而polyval P A 的含义是 A A A 5 A A 8 ones size A 例6 15仍以多项式x4 8x3 10为例 取一个2 2矩阵为自变量分别用polyval和polyvalm计算该多项式的值 4 多项式求根n次多项式具有n个根 当然这些根可能是实根 也可能含有若干对共轭复根 MATLAB提供的roots函数用于求多项式的全部根 其调用格式为 x roots P 其中P为多项式的系数向量 求得的根赋给向量x 即x 1 x 2 x n 分别代表多项式的n个根 例6 16求多项式x4 8x3 10的根 命令如下 A 1 8 0 0 10 x roots A 若已知多项式的全部根 则可以用poly函数建立起该多项式 其调用格式为 P poly x 若x为具有n个元素的向量 则poly x 建立以x为其根的多项式 且将该多项式的系数赋给向量P 例6 17已知f x 1 计算f x 0的全部根 2 由方程f x 0的根构造一个多项式g x 并与f x 进行对比 命令如下 P 3 0 4 5 7 2 5 X roots P 求方程f x 0的根G poly X 求多项式g x 6 2数值微积分 6 2 1数值微分1 数值差分与差商2 数值微分的实现在MATLAB中 没有直接提供求数值导数的函数 只有计算向前差分的函数diff 其调用格式为 DX diff X 计算向量X的向前差分 DX i X i 1 X i i 1 2 n 1 DX diff X n 计算X的n阶向前差分 例如 diff X 2 diff diff X DX diff A n dim 计算矩阵A的n阶差分 dim 1时 缺省状态 按列计算差分 dim 2 按行计算差分 例6 18设x由 0 2 间均匀分布的10个点组成 求sinx的1 3阶差分 命令如下 X linspace 0 2 pi 10 Y sin X DY diff Y 计算Y的一阶差分D2Y diff Y 2 计算Y的二阶差分 也可用命令diff DY 计算D3Y diff Y 3 计算Y的三阶差分 也可用diff D2Y 或diff DY 2 例6 19用不同的方法求函数f x 的数值导数 并在同一个坐标系中做出f x 的图像 程序如下 f inline sqrt x 3 2 x 2 x 12 x 5 1 6 5 x 2 g inline 3 x 2 4 x 1 sqrt x 3 2 x 2 x 12 2 1 6 x 5 5 6 5 x 3 0 01 3 p polyfit x f x 5 用5次多项式p拟合f x dp polyder p 对拟合多项式p求导数dpdpx polyval dp x 求dp在假设点的函数值dx diff f x 3 01 0 01 直接对f x 求数值导数gx g x 求函数f的导函数g在假设点的导数plot x dpx x dx x gx 作图 6 2 2数值积分1 数值积分基本原理求解定积分的数值方法多种多样 如简单的梯形法 辛普生 Simpson 法 牛顿 柯特斯 Newton Cotes 法等都是经常采用的方法 它们的基本思想都是将整个积分区间 a b 分成n个子区间 xi xi 1 i 1 2 n 其中x1 a xn 1 b 这样求定积分问题就分解为求和问题 2 数值积分的实现 1 被积函数是一个解析式MATLAB提供了quad函数和quadl函数来求定积分 它们的调用格式为 quad filename a b tol trace quadl filename a b tol trace 例6 20用两种不同的方法求定积分 先建立一个函数文件ex m functionex ex x ex exp x 2 然后在MATLAB命令窗口 输入命令 formatlongI quad ex 0 1 注意函数名应加字符引号I 0 74682418072642I quadl ex 0 1 I 0 74682413398845也可不建立关于被积函数的函数文件 而使用语句函数 内联函数 求解 命令如下 g inline exp x 2 定义一个语句函数g x exp x 2 I quadl g 0 1 注意函数名不加 号I 0 74682413398845formatshort 2 被积函数由一个表格定义在科学实验和工程应用中 函数关系往往是不知道的 只有实验测定的一组样本点和样本值 这时 就无法使用quad函数计算其定积分 在MATLAB中 对由表格形式定义的函数关系的求定积分问题用trapz X Y 函数 其中向量X Y定义函数关系Y f X X Y是两个等长的向量 X x1 x2 xn Y y1 y2 yn 并且x1 x2 xn 积分区间是 x1 xn 例6 21用trapz函数计算定积分 在MATLAB命令窗口 输入命令 X 0 0 01 1 Y exp X 2 trapz X Y ans 0 7468 3 二重积分数值求解使用MATLAB提供的dblquad函数就可以直接求出上述二重定积分的数值解 该函数的调用格式为 I dblquad f a b c d tol trace 该函数求f x y 在 a b c d 区域上的二重定积分 参数tol trace的用法与函数quad完全相同 调用函数quad8求定积分 formatlong fx inline exp x I n quad8 fx 1 2 5 1e 10 I 0 28579444254754n 33 例6 22计算二重定积分 1 建立一个函数文件fxy m functionf fxy x y globalki ki ki 1 ki用于统计被积函数的调用次数f exp x 2 2 sin x 2 y 2 调用dblquad函数求解 globalki ki 0 I dblquad fxy 2 2 1 1 kiI 1 57449318974494ki 1038 6 3离散傅立叶变换6 3 1离散傅立叶变换算法简要6 3 2离散傅立叶变换的实现一维离散傅立叶变换函数 其调用格式与功能为 1 fft X 返回向量X的离散傅立叶变换 设X的长度 即元素个数 为N 若N为2的幂次 则为以2为基数的快速傅立叶变换 否则为运算速度很慢的非2幂次的算法 对于矩阵X fft X 应用于矩阵的每一列 2 fft X N 计算N点离散傅立叶变换 它限定向量的长度为N 若X的长度小于N 则不足部分补上零 若大于N 则删去超出N的那些元素 对于矩阵X 它同样应用于矩阵的每一列 只是限定了向量的长度为N 3 fft X dim 或fft X N dim 这是对于矩阵而言的函数调用格式 前者的功能与FFT X 基本相同 而后者则与FFT X N 基本相同 只是当参数dim 1时 该函数作用于X的每一列 当dim 2时 则作用于X的每一行 值得一提的是 当已知给出的样本数N0不是2的幂次时 可以取一个N使它大于N0且是2的幂次 然后利用函数格式fft X N 或fft X N dim 便可进行快速傅立叶变换 这样 计算速度将大大加快 相应地 一维离散傅立叶逆变换函数是ifft ifft F 返回F的一维离散傅立叶逆变换 ifft F N 为N点逆变换 ifft F dim 或ifft F N dim 则由N或dim确定逆变换的点数或操作方向 例6 23给定数学函数x t 12sin 2 10t 4 5cos 2 40t 取N 128 试对t从0 1秒采样 用fft作快速傅立叶变换 绘制相应的振幅 频率图 在0 1秒时间范围内采样128点 从而可以确定采样周期和采样频率 由于离散傅立叶变换时的下标应是从0到N 1 故在实际应用时下标应该前移1 又考虑到对离散傅立叶变换来说 其振幅 F k 是关于N 2对称的 故只须使k从0到N 2即可 程序如下 N 128 采样点数T 1 采样时间终点t linspace 0 T N 给出N个采样时间ti I 1 N x 12 sin 2 pi 10 t pi 4 5 cos 2 pi 40 t 求各采样点样本值xdt t 2 t 1 采样周期f 1 dt 采样频率 Hz X fft x 计算x的快速傅立叶变换XF X 1 N 2 1 F k X k k 1 N 2 1 f f 0 N 2 N 使频率轴f从零开始plot f abs F 绘制振幅 频率图xlabel Frequency ylabel F k 6 4线性方程组求解6 4 1直接解法1 利用左除运算符的直接解法对于线性方程组Ax b 可以利用左除运算符 求解 x A b 例6 24用直接解法求解下列线性方程组 命令如下 A 2 1 5 1 1 5 0 7 0 2 1 1 1 6 1 4 b 13 9 6 0 x A b 2 利用矩阵的分解求解线性方程组矩阵分解是指根据一定的原理用某种算法将一个矩阵分解成若干个矩阵的乘积 常见的矩阵分解有LU分解 QR分解 Cholesky分解 以及Schur分解 Hessenberg分解 奇异分解等 1 LU分解矩阵的LU分解就是将一个矩阵表示为一个交换下三角矩阵和一个上三角矩阵的乘积形式 线性代数中已经证明 只要方阵A是非奇异的 LU分解总是可以进行的 MATLAB提供的lu函数用于对矩阵进行LU分解 其调用格式为 L U lu X 产生一个上三角阵U和一个变换形式的下三角阵L 行交换 使之满足X LU 注意 这里的矩阵X必须是方阵 L U P lu X 产生一个上三角阵U和一个下三角阵L以及一个置换矩阵P 使之满足PX LU 当然矩阵X同样必须是方阵 实现LU分解后 线性方程组Ax b的解x U L b 或x U L Pb 这样可以大大提高运算速度 例6 25用LU分解求解例6 24中的线性方程组 命令如下 A 2 1 5 1 1 5 0 7 0 2 1 1 1 6 1 4 b 13 9 6 0 L U lu A x U L b 或采用LU分解的第2种格式 命令如下 L U P lu A x U L P b 2 QR分解对矩阵X进行QR分解 就是把X分解为一个正交矩阵Q和一个上三角矩阵R的乘积形式 QR分解只能对方阵进行 MATLAB的函数qr可用于对矩阵进行QR分解 其调用格式为 Q R qr X 产生一个一个正交矩阵Q和一个上三角矩阵R 使之满足X QR Q R E qr X 产生一个一个正交矩阵Q 一个上三角矩阵R以及一个置换矩阵E 使之满足XE QR 实现QR分解后 线性方程组Ax b的解x R Q b 或x E R Q b 例6 26用QR分解求解例6 24中的线性方程组 命令如下 A 2 1 5 1 1 5 0 7 0 2 1 1 1 6 1 4 b 13 9 6 0 Q R qr A x R Q b 或采用QR分解的第2种格式 命令如下 Q R E qr A x E R Q b 3 Cholesky分解如果矩阵X是对称正定的 则Cholesky分解将矩阵X分解成一个下三角矩阵和上三角矩阵的乘积 设上三角矩阵为R 则下三角矩阵为其转置 即X R R MATLAB函数chol X 用于对矩阵X进行Cholesky分解 其调用格式为 R chol X 产生一个上三角阵R 使R R X 若X为非对称正定 则输出一个出错信息 R p chol X 这个命令格式将不输出出错信息 当X为对称正定的 则p 0 R与上述格式得到的结果相同 否则p为一个正整数 如果X为满秩矩阵 则R为一个阶数为q p 1的上三角阵 且满足R R X 1 q 1 q 实现Cholesky分解后 线性方程组Ax b变成R Rx b 所以x R R b 例6 27用Cholesky分解求解例6 24中的线性方程组 命令如下 A 2 1 5 1 1 5 0 7 0 2 1 1 1 6 1 4 b 13 9 6 0 R chol A Errorusing cholMatrixmustbepositivedefinite命令执行时 出现错误信息 说明A为非正定矩阵 6 4 2迭代解法迭代解法非常适合求解大型系数矩阵的方程组 在数值分析中 迭代解法主要包括Jacobi迭代法 Gauss Serdel迭代法 超松弛迭代法和两步迭代法 1 Jacobi迭代法对于线性方程组Ax b 如果A为非奇异方阵 即aii 0 i 1 2 n 则可将A分解为A D L U 其中D为对角阵 其元素为A的对角元素 L与U为A的下三角阵和上三角阵 于是Ax b化为 x D 1 L U x D 1b与之对应的迭代公式为 x k 1 D 1 L U x k D 1b这就是Jacobi迭代公式 如果序列 x k 1 收敛于x 则x必是方程Ax b的解 Jacobi迭代法的MATLAB函数文件Jacobi m如下 function y n jacobi A b x0 eps ifnargin 3eps 1 0e 6 elseifnargin epsx0 y y B x0 f n n 1 end 例6 28用Jacobi迭代法求解线性方程组 设迭代初值为0 迭代精度为10 6 在命令中调用函数文件Jacobi m 命令如下 A 10 1 0 1 10 2 0 2 10 b 9 7 6 x n jacobi A b 0 0 0 1 0e 6 2 Gauss Serdel迭代法在Jacobi迭代过程中 计算时 已经得到 不必再用 即原来的迭代公式Dx k 1 L U x k b可以改进为Dx k 1 Lx k 1 Ux k b 于是得到 x k 1 D L 1Ux k D L 1b该式即为Gauss Serdel迭代公式 和Jacobi迭代相比 Gauss Serdel迭代用新分量代替旧分量 精度会高些 Gauss Serdel迭代法的MATLAB函数文件gauseidel m如下 function y n gauseidel A b x0 eps ifnargin 3eps 1 0e 6 elseifnargin epsx0 y y G x0 f n n 1 end 例6 29用Gauss Serdel迭代法求解下列线性方程组 设迭代初值为0 迭代精度为10 6 在命令中调用函数文件gauseidel m 命令如下 A 10 1 0 1 10 2 0 2 10 b 9 7 6 x n gauseidel A b 0 0 0 1 0e 6 例6 30分别用Jacobi迭代和Gauss Serdel迭代法求解下列线性方程组 看是否收敛 命令如下 a 1 2 2 1 1 1 2 2 1 b 9 7 6 x n jacobi a b 0 0 0 x n gauseidel a b 0 0 0 6 4 3求线性方程组的通解线性方程组的求解分为两类 一类是求方程组的惟一解即特解 另一类是求方程组的无穷解即通解 这里对线性方程组Ax b的求解理论作一个归纳 1 当系数矩阵A是一个满秩方阵时 方程Ax b称为恰定方程 方程有惟一解x A 1b 这是最基本的一种情况 一般用x A b求解速度更快 2 当方程组右端向量b 0时 方程称为齐次方程组 齐次方程组总有零解 因此称解x 0为平凡解 当系数矩阵A的秩小于n n为方程组中未知变量的个数 时 齐次方程组有无穷多个非平凡解 其通解中包含n rank A 个线性无关的解向量 用MATLAB的函数null A r 可求得基础解系 3 当方程组右端向量b 0时 系数矩阵的秩rank A 与其增广矩阵的秩rank A b 是判断其是否有解的基本条件 当rank A rank A b n时 方程组有惟一解 x A b或x pinv A b 当rank A rank A b n时 方程组有无穷多个解 其通解 方程组的一个特解 对应的齐次方程组Ax 0的通解 可以用A b求得方程组的一个特解 用null A r 求得该方程组所对应的齐次方程组的基础解系 基础解系中包含n rank A 个线性无关的解向量 当rank A rank A b 时 方程组无解 有了上面这些讨论 可以设计一个求解线性方程组的函数文件line solution m 在例中可以调用line solution m文件来解线性方程组 6 5非线性方程与最优化问题求解6 5 1非线性方程数值求解1 单变量非线性方程求解在MATLAB中提供了一个fzero函数 可以用来求单变量非线性方程的根 该函数的调用格式为 z fzero fname x0 tol trace 其中fname是待求根的函数文件名 x0为搜索的起点 一个函数可能有多个根 但fzero函数只给出离x0最近的那个根 tol控制结果的相对精度 缺省时取tol eps trace 指定迭代信息是否在运算中显示 为1时显示 为0时不显示 缺省时取trace 0 例6 33求f x 在x0 5和x0 1作为迭代初值时的零点 先建立函数文件fz m functionf fz x f x 1 x 5 然后调用fzero函数求根 fzero fz 5 以 5作为迭代初值ans 5 1926fzero fz 1 以1作为迭代初值ans 0 1926 2 非线性方程组的求解对于非线性方程组F X 0 用fsolve函数求其数值解 fsolve函数的调用格式为 X fsolve fun X0 option 其中X为返回的解 fun是用于定义需求解的非线性方程组的函数文件名 X0是求根过程的初值 option为最优化工具箱的选项设定 最优化工具箱提供了20多个选项 用户可以使用optimset命令将它们显示出来 如果想改变其中某个选项 则可以调用optimset 函数来完成 例如 Display选项决定函数调用时中间结果的显示方式 其中 off 为不显示 iter 表示每步都显示 final 只显示最终结果 optimset Display off 将设定Display选项为 off 例6 34求下列方程组在 1 1 1 附近的解并对结果进行验证 首先建立函数文件myfun m functionF myfun X x X 1 y X 2 z X 3 F 1 sin x y z 2 exp x F 2 x y z F 3 x y z 在给定的初值x0 1 y0 1 z0 1下 调用fsolve函数求方程的根 X fsolve myfun 1 1 1 optimset Display off X 0 0224 0 0224 0 0000 6 5 2无约束最优化问题求解在实际应用中 许多科学研究和工程计算问题都可以归结为一个最小化问题 如能量最小 时间最短等 MATLAB提供了3个求最小值的函数 它们的调用格式为 1 x fval fminbnd filename x1 x2 option 求一元函数在 xl x2 区间中的极小值点x和最小值fval 2 x fval fminsearch filename x0 option 基于单纯形算法求多元函数的极小值点x和最小值fval 3 x fval fminunc filename x0 option 基于拟牛顿法求多元函数的极小值点x和最小值fval MATLAB没有专门提供求函数最大值的函数 但只要注意到 f x 在区间 a b 上的最小值就是f x 在 a b 的最大值 所以fminbnd f x1 x2 返回函数f x 在区间 x1 x2 上的最大值 例6 36求函数在区间 10 1 和 1 10 上的最小值点 首先建立函数文件fx m functionf f x f x 1 x 5 上述函数文件也可用一个语句函数代替 f inline x 1 x 5 再在MATLAB命令窗口 输入命令 fminbnd fx 10 1 求函数在 10 1 内的最小值点和最小值fminbnd f 1 10 求函数在 1 10 内的最小值点 注意函数名f不用加 例6 37求函数f在 0 5 0 5 0 5 附近的最小值 建立函数文件fxyz m functionf fxyz u x u 1 y u 2 z u 3 f x y 2 x 4 z 2 y 2 z 在MALAB命令窗口 输入命令 U fmin fminsearch fxyz 0 5 0 5 0 5 求函数的最小值点和最小值 6 5 3有约束最优化问题求解MATLAB最优化工具箱提供了一个fmincon函数 专门用于求解各种约束下的最优化问题 该函数的调用格式为 x fval fmincon filename x0 A b Aeq beq Lbnd Ubnd NonF option 其中x fval filename x0和option的含义与求最小值函数相同 其余参数为约束条件 参数NonF为非线性约束函数的M文件名 如果某个约束不存在 则用空矩阵来表示 例6 38求解有约束最优化问题 首先编写目标函数M文件fop m functionf fop x f 0 4 x 2 x 1 2 x 2 2 x 1 x 2 1 30 x 1 3 再设定约束条件 并调用fmincon函数求解此约束最优化问题 x0 0 5 0 5 A 1 0 5 0 5 1 b 0 4 0 5 lb 0 0 option optimset option LargeScale off option Display off x f fmincon fop x0 A b lb option 6 6常微分方程的数值求解6 6 1龙格 库塔法简介6 6 2龙格 库塔法的实现基于龙格 库塔法 MATLAB提供了求常微分方程数值解的函数 一般调用格式为 t y ode23 fname tspan y0 t y ode45 fname tspan y0 其中fname是定义f t y 的函数文件名 该函数文件必须返回一个列向量 tspan形式为 t0 tf 表示求解区间 y0是初始状态列向量 t和y分别给出时间向量和相应的状态向量 例6 39设有初值问题 试求其数值解 并与精确解相比较 精确解为y t 1 建立函数文件funt m functionyp funt t y yp y 2 t 2 4 t 1 2 求解微分方程 t0 0 tf 10 y0 2 t y ode23 funt t0 tf y0 求数值解y1 sqrt t 1 1 求精确解t y y1 y为数值解 y1为精确值 显然两者近似 例6 40已知一个二阶线性系统的微分方程 绘制系统的时间响应曲线和相平面图 函数ode23和ode45是对一阶常微分方程组设计的 因此对高阶常微