前　言 PREFACE 科学运算问题是每个理工科学生和科技工作者在课程学习、科学研究与工程实践中常常会遇到的问题，不容回避。对于非纯数学专业的学生和研究者而言，从底层全面学习相关数学问题的求解方法并非一件简单的事情，也不易得出复杂问题的解。所以，利用当前最先进的计算机工具，高效、准确、创造性地求解科学运算问题是一种行之有效的方法，尤其能够满足理工科人士的需求。 作者曾试图在同一部著作中叙述各个数学分支典型问题的直接求解方法，通过清华大学出版社出版了《高等应用数学问题的 MATLAB求解》。该书从 2004年出版之后多次重印再版，并于 2018年出版了第 4版，还配套发布了全新的 MOOC课程 ①，一直受到广泛的关注与欢迎。首次 MOOC开课的选课人数接近 14000人，教材内容也被数万篇期刊文章和学位论文引用。 从作者首次使用 MATLAB语言算起，已经有 30余年的时间了，通过相关领域的研究、思考与一线教学实践，积累了大量的实践经验资料。这些不可能在一部著作中全部介绍，所以与清华大学出版社策划与编写了这套“薛定宇教授大讲堂”丛书，系统深入地介绍基于 MATLAB语言与工具的科学运算问题的求解方法。 本丛书不是原来版本的简单改版，而是作者通过十余年的经验和资料积累，全面贯穿“再认识”的思想写作此书，深度融合科学运算数学知识与基于 MATLAB的直接求解方法与技巧，力图更好地诠释计算机工具在每个数学分支的作用，帮助读者以不同的思维与视角了解工程数学问题的求解方法，创造性地得出问题的解。 本丛书卷 I可以作为学习 MATLAB入门知识的教材与参考书，也为读者深入学习与熟练掌握 MATLAB语言编程技巧，深度理解科学运算领域 MATLAB的应用奠定一个坚实的基础。后续每一卷试图对应于一个数学专题或一门数学课程进行展开。全丛书的写作贯穿“计算思维”的思想，深度探讨该数学专题的问题求解方法。本丛书既适合于学完相应的数学课程之后，深入学习利用计算机工具的科学 . MOOC网址：https://www.icourse163.org/learn/NEU-1002660001 ·ii·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 运算问题求解方法与技巧，也可作为相应数学课程同步学习的伴侣，在学习相应课程理论知识的同时，侧重于学习基于计算机的数学问题求解方法，从另一个角度观察、审视数学课程所学的内容，扩大知识面，更好地学习相应的数学课程。 本书是丛书的卷 IV。本书系统地介绍两大主题非线性代数方程求解与最优化技术，主要解决这两个领域的数值计算问题。本书首先介绍各种非线性代数方程的解析解方法与数值解方法，并介绍多解方程的求解问题。后续各章将介绍无约束最优化、线性规划与二次型规划、非线性规划、混合整数规划、多目标规划与动态规划的基本概念与求解方法，侧重于求取最优化问题全局最优解的探讨与实践。本书还将介绍一些常用的智能优化方法，并通过一些具体的例子，对智能优化方法的效果作了必要的对比研究，得出有益的结论。 值此丛书付梓之际，衷心感谢相濡以沫的妻子杨军教授，她数十年如一日的无私关怀是我坚持研究、教学与写作工作的巨大动力。 薛定宇 2019年 9月 目　录 CONTENTS 第 1章方程求解与最优化技术 ································ 1 1.1方程与方程求解 ····································· 1 1.2最优化问题的起源与发展 ······························ 2 1.3本书框架 ··········································· 4本章习题 ················································ 5 第 2章代数方程的求解······································ 6 2.1多项式方程的求解···································· 6 2.1.1一次方程与二次方程 ···························· 7 2.1.2三次方程的解析解 ······························ 8 2.1.3四次方程的解析解 ······························ 9 2.1.4高次代数方程与 Abel–Ruffini定理·················· 11 2.2非线性方程的图解法 ·································· 11 2.2.1光滑隐函数曲线的绘制 ·························· 11 2.2.2一元方程的图解法 ······························ 12 2.2.3二元方程的图解法 ······························ 14 2.2.4方程的孤立解·································· 16 2.3代数方程的数值求解 ·································· 16 2.3.1 Newton–Raphson迭代方法 ······················· 16 2.3.2 MATLAB的直接求解函数························ 21 2.3.3求解精度的设置································ 23 2.3.4方程的复域求解································ 24 2.4联立方程组的精确求解 ································ 25 2.4.1低阶多项式方程的解析求解······················· 26 2.4.2多项式型方程的准解析解 ························ 28 2.4.3高次多项式矩阵方程的准解析解 ··················· 30 2.4.4非线性代数方程的准解析解······················· 32 ·iv·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 2.5多解矩阵方程的求解 ·································· 33 2.5.1方程求解思路与一般求解函数····················· 33 2.5.2伪多项式方程的求解 ···························· 37 2.5.3高精度求解函数································ 38 2.6欠定方程的求解 ····································· 40本章习题 ················································ 41第 3章无约束最优化 ······································· 44 3.1无约束最优化问题简介 ································ 44 3.1.1无约束最优化问题的数学模型····················· 45 3.1.2无约束最优化问题的解析解求解 ··················· 45 3.1.3无约束最优化问题的图解法······················· 45 3.1.4局部最优解与全局最优解 ························ 46 3.1.5数值求解算法的 MATLAB实现···················· 48 3.2无约束最优化问题的 MATLAB直接求解 ·················· 50 3.2.1直接求解方法·································· 50 3.2.2最优化控制选项································ 52 3.2.3附加参数的传递································ 56 3.2.4最优搜索的中间过程 ···························· 57 3.2.5最优化问题的结构体描述 ························ 58 3.2.6梯度信息与求解精度 ···························· 59 3.2.7离散点最优化问题的求解 ························ 62 3.2.8最优化问题的并行求解 ·························· 63 3.3全局最优解的尝试···································· 64 3.4带有决策变量边界的最优化问题························· 67 3.4.1单变量最优化问题 ······························ 67 3.4.2多变量最优化问题 ······························ 68 3.4.3边界问题全局最优解的尝试······················· 70 3.5最优化问题应用举例 ·································· 70 3.5.1线性回归问题的求解 ···························· 70 3.5.2曲线的最小二乘拟合 ···························· 71 3.5.3边值微分方程的打靶求解 ························ 75 3.5.4方程求解问题转换为最优化问题 ··················· 77本章习题 ················································ 78 第 4章线性规划与二次型规划 ································ 82 4.1线性规划问题简介···································· 83 4.1.1线性规划问题的数学模型 ························ 83 4.1.2二元线性规划的图解法 ·························· 84 4.1.3单纯形法简介·································· 85 4.2线性规划问题的直接求解 ······························ 88 4.2.1线性规划问题的求解函数 ························ 88 4.2.2多决策变量向量的线性规划问题 ··················· 93 4.2.3双下标的线性规划问题 ·························· 94 4.2.4线性规划的应用举例运输问题 ················· 95 4.3基于问题的线性规划描述与求解························· 98 4.3.1线性规划的 MPS文件描述························ 98 4.3.2基于问题的线性规划描述 ························ 100 4.3.3线性规划问题的转换 ···························· 104 4.4二次型规划问题的求解 ································ 106 4.4.1二次型规划的数学模型 ·························· 106 4.4.2二次型规划的直接求解 ·························· 106 4.4.3基于问题的二次型规划描述······················· 107 4.4.4双下标二次型规划 ······························ 111 4.5线性矩阵不等式问题 ·································· 112 4.5.1线性矩阵不等式的一般描述······················· 112 4.5.2 Lyapunov不等式 ······························· 113 4.5.3线性矩阵不等式问题分类 ························ 114 4.5.4线性矩阵不等式问题的 MATLAB求解 ·············· 115 4.5.5基于 YALMIP工具箱的最优化求解方法 ············· 117 4.5.6非凸最优化问题求解的尝试······················· 119 4.5.7带有二次型约束条件问题的求解 ··················· 120本章习题 ················································ 121 第 5章非线性规划 ········································· 126 5.1非线性规划简介 ····································· 127 5.1.1一般非线性规划问题的数学模型 ··················· 127 5.1.2可行解区域与图解法 ···························· 127 5.1.3数值求解方法举例 ······························ 129 ·vi·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 5.2非线性规划问题的直接求解 ···························· 131 5.2.1 MATLAB的直接求解函数························ 131 5.2.2搜索过程提前结束的处理 ························ 136 5.2.3梯度信息的利用································ 137 5.2.4多决策变量问题的求解 ·························· 138 5.2.5复杂非线性规划问题 ···························· 140 5.3非线性规划的全局最优解探讨 ·························· 141 5.3.1全局最优解的尝试 ······························ 142 5.3.2非凸二次型规划问题的全局寻优 ··················· 143 5.3.3凹费用运输问题的全局寻优······················· 146 5.3.4全局最优化求解程序的测试······················· 147 5.3.5分段目标函数的处理 ···························· 148 5.4双层规划问题 ······································· 150 5.4.1双层线性规划问题的求解 ························ 151 5.4.2双层二次型规划问题 ···························· 151 5.4.3基于 YALMIP工具箱的双层规划问题直接求解········ 152 5.5非线性规划应用举例 ·································· 154 5.5.1圆内最大面积的多边形 ·························· 154 5.5.2半无限规划问题································ 157 5.5.3混合池最优化问题 ······························ 160 5.5.4热交换网络的优化计算 ·························· 162 5.5.5基于最优化技术的非线性方程求解 ················· 165本章习题 ················································ 166 第 6章混合整数规划 ······································· 171 6.1整数规划简介 ······································· 171 6.1.1整数规划与混合整数规划 ························ 171 6.1.2整数规划问题的计算复杂度······················· 172 6.2穷举方法 ··········································· 173 6.2.1整数规划的穷举方法 ···························· 173 6.2.2离散规划问题·································· 176 6.2.3 0.1规划的穷举方法 ···························· 176 6.2.4混合整数规划的尝试 ···························· 178 6.3混合整数规划问题的求解 ······························ 181 6.3.1混合整数线性规划 ······························ 181 6.3.2整数规划问题的 LMI求解方法 ···················· 183 6.3.3混合整数非线性规划 ···························· 184 6.3.4一类离散规划问题的求解 ························ 186 6.3.5一般离散规划问题的求解 ························ 187 6.4 0.1混合整数规划的求解 ······························ 189 6.4.1 0.1线性规划问题的求解························· 189 6.4.2 0.1非线性规划问题的求解 ······················· 192 6.5混合整数规划应用···································· 194 6.5.1最优用料问题·································· 194 6.5.2指派问题 ····································· 195 6.5.3旅行商问题 ··································· 196 6.5.4背包问题 ····································· 200 6.5.5数独的填写 ··································· 201本章习题 ················································ 204第 7章多目标规划 ········································· 208 7.1多目标规划简介 ····································· 208 7.1.1多目标规划的背景介绍 ·························· 208 7.1.2多目标规划的数学模型 ·························· 209 7.1.3多目标规划问题的图解举例······················· 209 7.2多目标规划转换成单目标规划问题 ······················· 212 7.2.1无约束多目标函数的最小二乘求解 ················· 212 7.2.2线性加权变换及求解 ···························· 213 7.2.3线性规划问题的最佳妥协解······················· 215 7.2.4线性规划问题的最小二乘解······················· 216 7.3 Pareto最优解 ······································· 217 7.3.1多目标规划解的不唯一性 ························ 217 7.3.2解的占优性与 Pareto前沿 ························ 218 7.3.3 Pareto前沿的计算 ······························ 219 7.4极小极大问题求解···································· 220本章习题 ················································ 226第 8章动态规划与最优路径 ·································· 228 8.1动态规划简介 ······································· 228 8.1.1动态规划的基本概念与数学模型 ··················· 228 8.1.2线性规划问题的动态规划求解演示 ················· 229 ·viii·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 8.2有向图的路径寻优···································· 230 8.2.1有向图应用举例································ 230 8.2.2有向图最短路径问题的手工求解 ··················· 231 8.2.3逆序递推问题的动态规划表示····················· 232 8.2.4图的矩阵表示方法 ······························ 233 8.2.5有向图搜索及图示 ······························ 234 8.2.6 Dijkstra最短路径算法及实现 ····················· 237 8.3无向图的路径最优搜索 ································ 239 8.3.1无向图的矩阵描述 ······························ 239 8.3.2绝对坐标节点的最优路径规划算法与应用············ 240本章习题 ················································ 242 第 9章智能优化方法 ······································· 244 9.1智能优化方法简介···································· 244 9.1.1遗传算法简介·································· 245 9.1.2粒子群优化算法································ 246 9.2 MATLAB全局优化工具箱 ····························· 246 9.3最优化问题求解举例与对比研究························· 248 9.3.1无约束最优化问题 ······························ 248 9.3.2有约束最优化问题 ······························ 251 9.3.3混合整数规划问题求解 ·························· 257 9.3.4基于遗传算法的离散规划问题····················· 259本章习题 ················································ 261参考文献··················································· 262 MATLAB函数名索引 ········································· 265 术语索引··················································· 269 方程求解与最优化技术 第 1章 MATLAB最优化计算 最优化技术是科学与工程领域重要的数学工具，也是解决科学与工程问题的有效手段。毫不夸张地说，学会了最优化问题的理念与求解方法，可以将科研的水平提高一个档次，因为原来解决问题得到一个解就满足了，学会了最优化的思想后，很自然地将追求问题最好的解。 最优化问题与方程求解是密不可分的，所以这里首先回顾方程求解问题的发展简史，再回顾最优化与数学规划问题领域的发展过程，最后将简要地介绍本书各章的主要内容。 1.1方程与方程求解 方程（equation）是无处不在的数学模型，是在工程、科学与人们的日常生活中随时都能看到的数学模型。 　　定义1-1方程是含有一个或多个变量的等式。这些变量又称为未知变量，而这些满足等式的未知变量的值又称为方程的解。 　　定义1-2如果同时给出若干个方程，这些方程含有多个不同的变量，并要求这些方程同时成立，则这些方程称为联立方程（simultaneous equations）。 现代数学是用表达式和等号描述方程的，等号的左边有一个表达式，右侧有一个表达式，用等号连接这两个表达式。威尔士物理学家、数学家 Robert Recorde（约 1512.1558，图 1-1（a） 在 1557年发明了等号，并用数学符号描述了方程。 方程分为代数方程与微分方程，代数方程中变量与变量之间的关系是静态的，也就是说，方程的根是常数；而微分方程中，变量之间的关系是动态的。微分方程将在卷 V专门介绍，本书暂不涉及。代数方程分为线性方程、多项式方程和非线性方程。除此之外，还有参数方程、隐式方程等。线性代数方程在卷 III中已经给出了全面 介绍，本书将介绍其他方程的求解方法。 ·2·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 其实早在 Recorde使用等号描述方程之前，对诸多类型方程的研究就已经开始了。古巴比伦人在大约公元前 2000年就开始研究一元二次方程，公元 628年，古印度数学家 Brahmagupta（约 598.约 668）用语言而不是用数学公式描述了一元二次方程的求根方法。中国古代数学家刘徽（约 225.295）、王孝通（580.640）研究了一元三次方程。1554年，意大利数学家 Gerolamo Cardano（1501.1576，图 1-1 （b） 出版了一部数学著作，给出了意大利数学家 Scipione del Ferro（1465.1526）的一元三次方程的求根公式和意大利数学家 Lodovico de Ferrari（1522.1565）的一元四次方程的求根公式，Cardano还是第一个使用负数的数学家。挪威数学家 Niels Henrik Abel（1802.1829，图 1-1（c） 在 1824年证明了五次或五次以上的多项式方程是没有一般代数解法的。 （a）Robert Recorde（b）Gerolamo Cardano（c）Niels Henrik Abel 图 1-1 Recorde、Cardano和 Abel画像注：图像均来源于维基百科 还有一类特殊的方程，对所有的未知数都成立，这类方程称为恒等式（identity），例如 x 2 . y 2 =(x . y)(x + y), sin2 γ + cos2 γ =1（1-1-1） 这类方程一般都可以直接通过符号运算方式证明，所以本书不再探讨恒等式问题。 1.2最优化问题的起源与发展 最优化的理念起源于微积分研究的初期，法国数学家 Pierre de Fermat（1607. 1665，图 1-2（a） 与法国数学家 Joseph-Louis Lagrange（1736.1813，图 1-2（b） 分别提出了基于微积分的公式求解最优值的方法。除了简单的函数最优化问题之外，一般又统称最优化问题为数学规划（mathematical programming）问题。相关的历史回顾可以参见文献 [1]。 苏联学者 Leonid Vitaliyevich Kantorovich（1912.1986，图 1-2（c））在最优 （a）Pierre de Fermat（b）Joseph-Louis Lagrange（c）Leonid Kantorovich图 1-2 Fermat、Lagrange和 Kantorovich画像或照片注：图像均来源于维基百科 化领域尤其是线性规划领域做了大量的奠基工作，美国数学家 George Bernard Dantzig（1914.2005，图 1-3（a） 提出了著名的单纯形法（simplex method）[1]，求解线性规划问题。Dantzig提出单纯形法的背后有一个有趣的故事 [2,3]。1939年加州大学 Berkeley分校的博士生 Dantzig有一次上课迟到，错把老师 Jerzy Neyman （1894.1981）在黑板上写的两个世界数学难题当成课后作业，给出了问题的解，为线性规划问题提出了一种高效的求解方法。美国数学家、计算机学家 John von Neumann（1903.1957，图 1-3（b） 提出了对偶理论与计算方法，进一步提高了线性规划问题的求解效率。 美国应用数学家 Richard Ernest Bellman（1920.1984，图 1-3（c） 开创了一个最优化问题的新领域动态规划，实现了多级决策的规划问题。 （a）George Dantzig（b）John von Neumann（c）Richard Bellman 图 1-3 Dantzig、von Neumann和 Bellman的照片注：图像均来源于维基百科 最优化理念与技术为许多科学与工程领域奠定了数学基础，“最优”一词可以 ·4·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 和任何一个领域联用，为其注入新的活力。例如，在搜索引擎上搜索“最优”可以搜索到很多相关的领域，如最优控制、最优设计、最优系统、资源最优配置、最优停止理论、最优资本结构等，这些领域都是和最优化密切相关的。 1.3本书框架 最优化问题是与代数方程密切相关的，所以本书在第 2章中深入探讨了各类代数静态方程的求解方法，包括多项式方程的解析解法、非线性方程组的图解方法与基于搜索的复杂非线性方程数值求解方法。特别地，还探讨了多解非线性矩阵方程全部数值解与准解析解的方法，理论上可以求解任意复杂的非线性方程组。 第 3章介绍了最简单的一类最优化问题无约束最优化问题的求解方法，包括最优化问题的解析求解规则、简单最优化问题的图解法、基于梯度信息的求解方法等，着重介绍基于 MATLAB求解函数的直接求解方法，给出了局部最优解与全局最优解的概念，并编写了试图求解全局最优解的 MATLAB通用工具。本章还探讨了最优化技术在线性回归问题、曲线的最小二乘拟合与边值微分方程的打靶求解方面的应用。 一般数学规划问题分为线性规划问题、非线性规划问题、混合整数规划问题、多目标规划问题与动态规划问题等，本书后续内容也按照这样的分类分别介绍各种数学规划问题的求解方法。 第 4章侧重于介绍凸优化问题线性规划与二次型规划问题的求解方法，主要介绍基于 MATLAB现成工具的直接求解方法，还介绍了新版 MATLAB支持的基于问题的描述与求解方法，使得复杂线性规划与二次型规划问题的描述与求解更直接、更容易。除此之外，还介绍了线性矩阵不等式问题求解方法。 第 5章主要介绍了非线性规划问题的求解方法。首先介绍简单问题的图解法，然后介绍基于 MATLAB的非线性规划问题求解函数与复杂问题的描述与求解方法，特别地，作者编写了求解非线性非凸优化问题全局最优解的通用工具。本章还探讨了圆内最大面积的多边形、半无限规划问题、热交换网络的优化计算等应用问题的求解方法。 第 6章介绍了混合整数规划问题的求解。探讨了小规模整数规划问题的穷举方法，还介绍了线性混合整数规划问题、非线性混合整数规划问题及混合 0.1规划问题的求解方法，并探讨了最优用料问题、指派问题、背包问题等应用问题的求解方法，还介绍了基于整数 0.1规划的旅行商问题、数独问题的建模与求解方法。 第 7章侧重于多目标规划问题的求解方法，给出了多目标规划的数学模型并探讨多目标规划问题的图解方法，另外，侧重于介绍如何将多目标规划问题转换成普通最优化问题直接求解的方法。本章还给出了 Pareto解集的概念，并介绍了极小极大问题的求解方法。 第 8章简单探讨了动态规划问题的建模与求解方法，侧重介绍有向图最短路径问题的求解方法，还探讨了无向图的路径最优问题求解方法。 传统的最优化求解方法主要是基于搜索的方法，有时可能得出问题的局部最优解。第 9章将简要地介绍基于 MATLAB的智能优化方法，如遗传算法、粒子群优化算法、模式搜索算法与模拟退火方法等，并通过算例对比研究智能优化方法与传统优化方法，得出有意义的结论。 本章习题 1.1看看能否利用 MATLAB求解下面的方程。 { （1） x1 + x2 = 35 2x1 +4x2 = 94 { 2e.xy /2 （2） x2+ e.x/2 sin(xy)=0 y2 cos(x + y2)+ x2ex+y =0 1.2试在 .2 . x . 11区间内找出下面函数的最小值 [4]。 52397179 1 54 f(x)= x 6 . x + x + x 3 . x 2 . x + 25801020 10 1.3绘制下面函数的曲面，试旋转观察曲面，找出 x和 y取何值时曲面能达到谷底。 （1）f(x, y)= .(y + 47) sin √ x 2 +(y + 47) . x sin √ x . (y + 47) ()2 ()2 [()()] xy xy （2）f(x, y) = 20+30 . 1+20 . 1.10 cos 30 . 1巾 + cos 20 . 1巾 代数方程的求解 第 2章 MATLAB最优化计算 方程在人们日常生活、科学研究与工程实践中都是经常遇到的数学模型，所谓方程就是含有未知数的等式。方程是用来描述变量与变量之间数学关系的。方程分为代数方程、微分方程等，本章主要探讨代数方程的求解方法，兼顾代数方程的解析解方法与数值解方法，并试图得出多解方程全部的解。 本书卷 III侧重于探讨多元一次线性方程的求解方法，不但能求解 AX = B类简单线性代数方程的唯一解、无穷解与最小二乘解，还可以求解 XA = B， AXB = C及其多项型线性代数方程的解，此外还给出了一般 Sylvester方程及多项 Sylvester方程的求解方法。上述方程均可以利用 MATLAB的强大功能求取出数值解与解析解。有关线性方程的求解方法可以参见卷 III的相关内容。 本章侧重于介绍多项式方程与一般非线性方程的求解方法。 2.1节主要探讨低阶多项式方程的求解公式，并给出底层的 MATLAB实现程序，从数值运算角度看，该程序尤其适合于含有重根的低阶多项式方程的求解。 2.2节介绍一般一元与二元方程的图解方法，并介绍方程的实际求解方法。 2.3节介绍一般非线性方程组的数值求解方法。首先介绍经典的 Newton–Raphson迭代方法及 MATLAB实现，然后介绍 MATLAB提供的非线性代数方程与矩阵方程的求解方法。 2.4节介绍基于符号运算的低次代数方程解析解方法，然后介绍高次代数方程与非线性矩阵方程的准解析解方法。 2.5节介绍多解矩阵方程的求解方法、伪多项式方程的求解方法，并介绍高精度求解方法与实现。2.6节还将探讨欠定方程的求解方法。 2.1多项式方程的求解 多项式方程是实际应用中经常遇到的方程，本节先介绍多项式方程的数学形式，再介绍多项式方程的求解方法。　　定义 2-1多项式方程的一般形式为 x n + a1x n.1 + a2x n.2 + ··· + an.1x + an =0（2-1-1） . . . .. . 多项式方程根与系数的关系满足如下的 Viète定理，又称 Viète公式，该定理是以法国数学家 Fran.ois Viète（1540.1603）命名的。　　定理 2-1（Viète定理）假设多项式方程的根为 x1，x2，···，xn，则有 x1 + x2 + ··· + xn = .a1 (x1.x2 +x1x3 +···+x1xn)+(x2x1 +x2x3 +···+x2xn)+··· = a2 （2-1-2） . . 一次与二次方程都有很简单的求解公式，这里将通过例子演示一元方程与二{ x1x2 ··· xn =(.1)nan 本节侧重介绍低阶多项式方程的求根公式及其 MATLAB实现，并给出高阶多项式方程的 Abel–Ruffini定理。 2.1.1一次方程与二次方程 元方程的直接求解方法。 　　例 2-1一次多项式方程 x + c =0的解是什么？ 　　解显然，一次多项式方程的解是 x = .c，不论 c是何值。　　例 2-2公元四至五世纪的中国古代著名的数学著作《孙子算经》曾给出了鸡兔同笼问题：“今有雉兔同笼，上有三十五头，下有九十四足，问雉兔各几何？”　　解古典数学著作中有各种各样的方法求解鸡兔同笼问题。若引入代数方程的思维，则假设鸡的个数为 x1，兔的个数为 x2，这样容易地列出下面的二元一次方程组： x1 + x2 = 35 2x1 +4x2 = 94由第一个方程，令 x1 = 35 . x2，将其带入第二个方程，有 2(35 . x2)+4x2 = 70 . 2x2 +4x2 = 70+2x2 = 94 立即得出 x2 = 12，代入则可以得出 x1 = 35 . x2 = 23。将根代入原始方程，则可以看出误差都为零，说明得出的根是正确的。当然，求解这类方程有许多种方法，这里暂时不去探讨了。 　　例 2-3一元二次方程 ax2 + bx + c =0的求解。　　解古巴比伦人早在公元前 1800年就开始研究这类问题，后来诸多数学家也在研究二次方程的求解方法，直到 1615年出版的法国数学家 Fran.ois Viète的著作中，给出了一般多项式方程根与系数之间的关系，即前面给出的 Viète定理。 提取方程两端的系数 a，则可以得出 a(x2 + bx/a + c/a)=0，如果 a不等于零，则探讨方程 x2 + bx/a + c/a =0的根就可以了。如果采用配方法，可以从底层推导一元二次 ·8·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 方程的求解公式。 2 bc x + x + aa ()2 ()2 bb bc 2 = x + x + . + a 2a 2aa ()2 bb2 . 4ac = x + . =0 2a 4a2 由最后一个方程显然可以得出 ()2 bb2 . 4ac x += 2a 4a2 两端开方，求代数平方根，再经过简单处理，则可以得出原方程的两个根： ∞ .b ± b2 . 4ac x1,2 = 2a 2.1.2三次方程的解析解 三次或三次以上代数方程的求解就没有这么简单了。古巴比伦人研究过三次多项式方程。中国三国时期著名数学家刘徽在 265年出版了《九章算术注》，其中描述了三次方程（ cubic equations）问题。中国唐代数学家王孝通在其著作《辑古算经》中建立并求解了 25个特殊三次方程。 　　定义 2-2三次方程的一般形式为 x3 + c2x2 + c3x + c4 =0。 为简单起见，该方程进行了首一化处理。令 x = t . c2/3，则可以将三次方程变换成 t3 + pt + q =0的形式，其中 p = . 1 c22 + c3,q =2 c 23 . 1 c2c3 + c4（2-1-3） 3 273 　　例 2-4试利用 MATLAB证实式（2-1-3）的叙述。　　解利用 MATLAB进行变量替换，则可以给出下面的语句： >> syms x c2 c3 c4; f=x^3+c2*x^2+c3*x+c4; syms t; f1=subs(f,x,t-c2/3); f2=collect(f1,t) 则可以立即得出结果如下，由此可以证实式（2-1-3）中的结果。 () 1 21 32 3 f2 = t+ . c2 + c3t + c2 . c3c2 + c4 =0 3 273 这样，变量 t三个根的闭式求解公式为 t = .k u + .2k v, k =0, 1, 2 （2-1-4） 其中，. =(.1+ ∞ 3j)/2，且 √ √ √ √ 23 23 qqpqqp u = 3 . 2+ 4+27 ,v = 3 . 2 . 4 +27（2-1-5） 不过在 MATLAB实际计算中，由上式独立计算 u和 v的值有时会出现错误，因为开方运算解是不唯一的。为准确计算 u和 v值，可以考虑先计算 u，再用 Viète定理计算 v的值，v = .p/(3u)。 当然，有一种特殊情况必须考虑，就是当 u =0时，v =0，这时方程的三重根为 x = .c2/3。 将得出的 t代入 x = t . c2/3，可以得出方程的三个根。这个闭式解公式是意大利数学家 Gerolamo Cardano在 1545年提出的。 2.1.3四次方程的解析解 一般四次方程（ quartic equation）可以先经过特殊的变量替换转换成特殊形式的四次方程，最后给出闭式求解公式。 　　定义 2-3四次方程的一般形式为 x4 + c2x3 + c3x2 + c4x + c5 =0。 　　例 2-5对一般四次方程，若令 x = y . c2/4，则原方程会变换成什么形式？　　解利用下面的语句可以进行方程的变换。 >> syms x c2 c3 c4 c5 y; f=x^4+c2*x^3+c3*x^2+c4*x+c5; f1=subs(f,x,y-c2/4); f2=collect(f1,y) 可以将方程变换成 y4 + py2 + qy + s =0的特殊四次方程形式，其中， 3 2 p = c3 . c 2 8 q =1 c 32 . 1 c3c2 + c4 （2-1-6） 82 31 1 42 s = . c2 + c3c2 . c4c2 + c5 25616 4 由给出的特殊形式，可以将 y的方程变换成下面的形式。 ∞ ∞ (y 2 + p 2+m. 2my+2 ∞q 2m)(y 2 + p 2+m+2my. 2 ∞q 2m)=0（2-1-7） 其中，m为三次方程 8m3 +8pm2 + (2p2 . 8s)m . q2 =0的根。这个算法是意大利数学家 Lodovico de Ferrari提出的。当然，.0。0， 这样做的前提是 m = 如果 m = 则意味着 q =0，这样 y的方程可以写成 y4 + py2 + s =0，易于求解 y2，再求出 y。有了 y，则可以由变换表达式求出相应的 x。 MATLAB提供的 roots()函数是基于矩阵特征值计算的多项式方程求解函数，具体的调用格式为 r=roots(p)，其中，p为降幂排列的多项式系数向量，r为方程的数值解，为列向量。 该函数在求解有重根方程时经常被诟病，因为得出的根误差比较大。综合考虑上面给出的低次方程求解算法，可以编写出 roots()函数的替代函数，该函数有 ·10·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 望得出精确的低次方程的数值解，而方程阶次高于 4时自动嵌入 MATLAB的原始 roots()函数，得出方程根的数值解。 function r=roots1(c) i=find(c~=0); c=c(i(1):end); n=length(c)-1; c=c/c(1); switch n case 1, r=-c(2); case 2 d=sqrt(c(2)^2-4*c(3)); r=[-c(2)+d; -c(2)-d]/2; case 3 p=c(3)-c(2)^2/3; q=2*c(2)^3/27-c(2)*c(3)/3+c(4); v=0; u=(-q/2+sqrt(q^2/4+p^3/27))^(1/3); if u~=0, v=-p/3/u; end xi=(-1+sqrt(3)*1i)/2; w=xi.^[0,1,2]'; r=w*u+w.^2*v-c(2)/3; case 4 p=-3*c(2)*c(2)/8+c(3); q=c(2)^3/8-0.5*c(2)*c(3)+c(4); s=-3*c(2)^4/256+c(2)*c(2)*c(3)/16-0.25*c(2)*c(4)+c(5); rr=roots1([8,8*p,2*p^2-8*s,-q^2]); m=rr(1); d=sqrt(2*m); if q==0, r1=roots1([1,p,s]); r=[sqrt(r1); -sqrt(r1)]; else r=[roots1([1,-d,p/2+m+q/2/d]); roots1([1,d,p/2+m-q/2/d])]; end, r=r-c(2)/4; otherwise, r=roots(c); end 　　例 2-6试求解方程 (s . 5)4 =0。　　解当然，如果以这种形式给出方程，则可以立即得出方程的四个重根，都是 5。在实际应用时经常已知方程的展开形式，不知道分解的形式，如何求解呢？可以尝试 MATLAB提供的 roots()函数，也可以尝试新编写的 roots1()函数。 >> P=conv([1,-5],conv([1,-5],conv([1,-5],[1,-5]))) d1=roots(P), d2=roots1(P), d3=roots(sym(P)) e1=norm(polyval(P,d1)), e2=norm(polyval(P,d2)) 可以看出， roots1()函数得出所期望的四重根 s =5，而 roots()函数得出的结果为 5.0010，5 ± 0.001j，4.9990，可以看出，方程有重根时 MATLAB函数数值求解有较大的误差，在双精度数据结构下，用其他数值算法也有同样的问题，在实际应用中可能引起麻烦。代入原方程后， e1 =2.1690×10.12，e2 =0。通过上面演示还可以看出，利用符号运算总能得出正确的结果。 　　例 2-7试在双精度框架下求解复系数多项式方程并评价精度。 f(s)= s 4 + (5 + 3j)s 3 + (6 + 12j)s 2 . (2 . 14j)s . (4 . 4j)=0 　　解到现在为止介绍了两个多项式方程的数值求解函数 roots()和 roots1()，其 实，构造原方程时，是假设方程的三重根为 .1+1j，还有一个实数根 .2。现在可以由下面的语句直接求解，得出结果后与已知根的解析解相比，可见， roots()函数得出的根的误差范数为 3.5296×10.5，而 roots1()函数的误差范数为 0，说明这里给出的方程的解是精确的，roots1()函数同样适合求解复系数多项式方程。 >> p=[1,5+3j,6+12i,-2+14j,-4+4j]; %复数标记 i与 j是相同的 r1=roots(p), err1=norm(r1-[-2; -1-1j; -1-1j; -1-1j]) r2=roots1(p), err2=norm(r2-[-1-1j; -1-1j; -1-1j; -2]) 2.1.4高次代数方程与 Abel–Ruffini定理 　　定义2-4方程的代数解法是利用有限次加减乘除、整数次乘方、开方运算可以构造出来的闭式求解公式。 　　定理 2-2（Abel–Ruffini定理）任意系数的五次或五次以上的多项式方程是没有代数解法的。该定理又称 Abel不可能性定理。 意大利数学家 Paolo Ruffini（1765.1822）在 1799年给出了该定理的不完全证明，挪威数学家 Niels Henrik Abel在 1824年给出了定理的证明。所以对一般高次多项式方程而言，只能采用数值解方法。 2.2非线性方程的图解法 如果方程含有一个或两个未知数，则可以通过图解法求解方程。如果未知数过多，则不适合使用图解法，而应该尝试其他方法。本节将介绍一元与二元方程的图解方法，并分析总结图解法的优势与劣势。 2.2.1光滑隐函数曲线的绘制 MATLAB提供了一般双变量隐函数模型的绘制函数 fimplicit()，正常情况下该函数的默认设置可以得出比较光滑的隐函数曲线，不过在一些特定的场合下，需要手工调节该函数的控制参数，以获得光滑的函数曲线。下面通过例子演示光滑隐函数曲线的绘制方法。 　　例 2-8试绘制隐函数 y2 cos(x + y2)+ x2ex+y =0在 (.2巾, 2巾)求解的光滑曲线。　　解下面的函数可以直接绘制出该隐函数的曲线，如图 2-1所示。可以看出，该函数可以直接绘制出隐函数曲线，不过曲线的某些地方比较粗糙，例如曲线左上角与右下角区域出现不光滑的毛刺现象。 >> syms x y; p=2*pi; f=y^2*cos(x+y^2)+x^2*exp(x+y); fimplicit(f,[-p,p])隐函数曲线的光滑度受网格密度属性 'Meshdensity'直接影响，其默认值为 151，如果发现得出的曲线不光滑，不妨将该值设置为大一些的值，如 500，这时得出的隐函 ·12·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 数曲线如图 2-2所示。可见，这样得出的曲线光滑度是令人满意的，即使作局部放大，曲线仍然是光滑的。 >> fimplicit(f,[-p,p],'Meshdensity',500) 图 2-1默认的隐函数曲线 图 2-2光滑的隐函数曲线 2.2.2一元方程的图解法　　定义 2-5一元方程的一般数学形式为 f(x)=0（2-2-1）对于任意的单变量方程 f(x)=0，可以考虑将方程用符号表达式描述，然后调用 fplot()函数绘制出方程的曲线，这样，就可以用图解方法求出方程曲线与横轴的交点，这些交点就是方程的解。 　　例 2-9根式方程的解析求解是有诸多条件的，如果条件不满足则无法解析求解。试用图解方法求解下面的根式方程。 √√√ √ 2x2 +3+ x2 +3x +2 . 2x2 . 3x +5 . x2 . 5x +2=0 　　解先用符号表达式表示方程的左端，则可以调用 fplot()函数，并叠印上横轴，则可以得出如图 2-3所示的曲线。从得出的结果看，方程与横轴只有一个交点，该交点就是方程的解。 >>symsx %声明符号变量 f=sqrt(2*x^2+3)+sqrt(x^2+3*x+2)-... sqrt(2*x^2-3*x+5)-sqrt(x^2-5*x+2); fplot(f), line([-5 5],[0,0]) %绘制曲线并叠印横轴 图 2-3一元方程的曲线与求解 如果想得到方程的解，则需要单击图形坐标轴工具栏中的 图标，对交点附近 的区域作局部放大，用户可以反复地使用放大功能，直至 x的标度都大致一样，这时，可以认为得到了方程的解。对这个具体的方程而言，通过局部放大得出方程的解为 x =0.13809878，如图 2-4所示，代入方程则可以得出误差为 6.018×10.9。 图 2-4局部放大的结果与方程的解 　　例 2-10试求解下面的一元超越方程。 e.0.2x sin(3x + 2) . cos x =0,x → (0, 20) 　　解这个超越方程是没有解析解的，必须用数值解的方法求解。图解法是一种数值 ·14·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 求解方法。前面的例子曾经使用符号表达式来描述原始方程，这里采用匿名函数来描述方程，这两种方法是等效的。由 fplot()函数可以直接绘制出方程的曲线，并同时叠印出横轴，如图 2-5所示。这样，方程曲线与横轴的交点都是方程的解。 >> f=@(x)exp(-0.2*x).*sin(3*x+2)-cos(x); %注意应该采用点运算 fplot(f,[0,20]), line([0,20],[0,0]) %绘制曲线并叠印横轴 图 2-5超越方程的图解法 可见，在给定的区间内，该方程曲线与实轴有 6个交点，这些交点处的 x都是方程的根。与前面的叙述一致，可以通过局部放大的方法逐一求出方程的根，不过求解过程还是很麻烦的。例如，可以求出方程的一个根为 x = 10.9601289577，代入方程可以得出误差为 .3.4739×10.11。后面将探讨更好的方法。 2.2.3二元方程的图解法 二元方程联立组的数学形式与定义在下面给出，本节将探讨利用图解的方式求解相应的二元联立方程组，并指出图解法存在的问题。 　　定义 2-6二元联立方程的一般形式为 { f(x, y)=0 （2-2-2）g(x, y)=0从给出方程的数学形式看， f(x, y)=0可以看成关于自变量 x和 y的隐函数表达式，故使用 fimplicit()函数即可以直接绘制该隐函数的曲线，而曲线上的所有点都满足该方程。同样， g(x, y)=0也是隐函数的数学表达式，由 fimplicit()函数可以求解该方程。如果将这两个函数在同一坐标系下绘制出来，得出的曲线交点则为联立方程的解。 　　例 2-11试求解二元方程在 .2巾 . x, y . 2巾区域内的解。 { 2e.xy /2 x2+ e.x/2 sin(xy)=0 y2 cos(x + y2)+ x2ex+y =0 　　解要想求解联立方程，可以声明符号变量 x和 y，然后将两个方程用符号表达式分别表示出来，再调用 fimplicit()就可以绘制出两个方程的解曲线了，如图 2-6所示。图中给出的每条曲线都满足一个方程，而联立方程的解为曲线的交点。可以看出，该方程在给定区域内有很多解。由得出的图形可见，隐函数曲线的光滑度不影响曲线交点的求解，所以可以按默认形式绘制曲线。 >> syms x y f1=x^2*exp(-x*y^2/2)+exp(-x/2)*sin(x*y); f2=y^2*cos(x+y^2)+x^2*exp(x+y); fimplicit([f1 f2],[-2*pi,2*pi]) 图 2-6联立方程图解法示意图 如果想得出某个具体交点的信息，则可以对该点作局部放大，大致地得出交点处的 x与 y值，不过可以预计，这样的解不会太精确，此外由于这个联立方程存在太多交点，所以一个一个局部放大去求解的方式显然不适用，应该考虑引入能一次性求出所有交点的全新方法。 　　例 2-12试用图解法求解下面的联立方程。 { 2 x + y2 =5 x +4y3 +3y2 =2 　　解先用符号表达式表示这两个方程，然后将这两个隐函数绘制出来，得出的曲线如图 2-7所示。可见，图中显示这两组曲线有两个交点。能因此得出结论，说原方程有两个根吗？ >> syms x y f1=x^2+y^2-5; f2=x+4*y^3+3*y^2-2; fimplicit([f1,f2],[-pi,pi]) 如果将第二个表达式稍加变换，则 x = .4y3 . 3y2 +2，代入第一个方程，有 16y 6 + 24y 5 +9y 4 . 16y 3 . 11y 2 . 1=0 ·16·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 图 2-7联立方程图解法示意图可见，这是一个关于 y的六次多项式方程，很可能该方程有六个根，而不是图 2-7中所示的两个根。为什么图中只给出两个根呢？因为原方程有两个实根，其他四个根应该是两对共轭复数根。在图解法中只能表示出方程的实数根，而不能显示、求取复数根。 2.2.4方程的孤立解 观察例 2-11中给出的方程，不难发现，将 x =0，y =0这个点代入两个方程，这两个方程都是满足的。如果从给出的曲线看，第一个方程的曲线似乎有意回避了这个点，而第二个方程也不经过这个点。这个点不是由曲线其他点演化而来的，这样的解称为孤立解（isolated solution）。 目前没有任何方法可以求取方程的孤立解，只能由用户的经验观察与判断某些点是不是方程的孤立解。例如，通过观察可见， x =0，y =0点是该方程的孤立解，代入原方程则可以验证其正确性。 2.3代数方程的数值求解 前面介绍的图解方法只是非线性方程组众多求解方法的一种，图解法有其明显的优势，但也有劣势。图解法只能用于求解一元或二元方程的实数根，对多元方程是不能采用图解法求解的。本节将探讨一般方程的求解思路与实用求解函数。 2.3.1 Newton–Raphson迭代方法 为简单起见，可以先探讨一元方程的求解方法。 Newton–Raphson迭代方法是以英国科学家 Isaac Newton与英国数学家 Joseph Raphson（约 1648.约 1715）命名的一般方程的迭代方法。 假设一元方程是由 f(x)=0描述的，且在 x = x0点处函数的值 f(x0)为已知的。这样，可以在 (x0,f(x0))点作函数曲线的一条切线，如图 2-8所示，则该切线与横轴的交点 x1可以认为是找到的方程的第一个近似的根。由图 2-8给出的斜率为 f ′ (x0)的切线方程，可以得出 x1的位置为 x1 = x0 . ff′((xx00))（2-3-1） 式中， f ′ (x0)为 f(x)关于 x的导函数在 x0点的值。再由 x1出发作切线，则可以得出 x2，以 x2出发找到 x3，···。若已经找到了 xk，则从该点出发可以搜索到下一个点。 xk+1 = xk . ff′((xxkk)) ,k =0, 1, 2, ···（2-3-2） 如果 |xk+1 . xk| . ξ1或 |f(xk+1)| . ξ2，其中，ξ1与 ξ2为预先选定的误差容限，则可以认为 xk为原方程的一个解。 图 2-8 Newton–Raphson迭代法示意图　　定义 2-7对多元向量函数 f(x)，其 Jacobi矩阵的定义为 . .... .... df1(x)/dx1 df1(x)/dx2 ··· df1(x)/dxn J[f(x)] =df2(x.. )/dx1 df2(x.. )/dx2 · ..·· df2(x.. )/dxn （2-3-3）.. . . dfm(x)/dx1 dfm(x)/dx2 ··· dfm(x)/dxn 　　定理 2-3更一般地，对多元方程 f(x)= 0，其中 x为向量或矩阵，而非线性函数 f(x)也是同维数的函数，则可以由下式迭代求出: xk+1 如果 ||xk+1 . xk|| . ξ1了 Jacobi矩阵。 = xk . J.1 [ 或 ||f(xk+1)|| <ξ2，则 xk f(xk) ] f(xk) （2-3-4） 为方程的根。在定义中使用 根据给出的算法，编写出 MATLAB的通用求解函数: function x=nr_sols(f,df,x0,epsilon,key) ·18·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 if nargin<=4, key=0; end, x1=x0; if nargin==3, epsilon=eps*10; end while (1), x=x0-df(x0)\f(x0); if norm(x-x0)> syms x; f=exp(-0.2*x)*sin(3*x+2)+cos(x), diff(f) 可以得出函数的导数为 f ′ (x)=3e.x/5 cos (3x + 2) . sin x . 1 e.x/5 sin (3x + 2) 5 有了函数及其导数，则可以用匿名函数表示它们，再设定搜索的初值为 x0 = 10，这时调用 Newton–Raphson求解算法，则可以得出方程的解搜索的中间点为 [10, 10.8809, 11.0700, 11.0593, 11.0593]T，方程的解为 x = 11.0593，将其代入原方程，则可以得出误差为 .5.1348×10.16。可以看出，利用这样的方法求解精度还是比较高的。整个求解过程的示意图如图 2-9所示，可以看出经过几步迭代就可以得出方程的解。 >> f=@(x)exp(-0.2*x).*sin(3*x+2)+cos(x); d=@(x)3*exp(-x/5).*cos(3*x+2)-sin(x)-(exp(-x/5).*sin(3*x+2))/5; x=nr_sols(f,d,10,1e-15,1), x1=x(end), f(x1) fplot(f,[9,12]), hold on, plot(x,zeros(size(x)),'o'), hold off 图 2-9一元方程求解的中间过程 　　例 2-14试用正割法重新求解例 2-13中方程的一个根。　　解选择两个初始点 x0 =9，x1 = 10，调用求解函数则可以得出求解过程的中间结果为 x = [9, 10, 12.9816, 11.3635, 10.7250, 11.0222, 11.0633, 11.0592, 11.0593, 11.0593]，其中搜索过程如图 2-10所示。可以看出，虽然这里给出的求解函数效率不如 Newton– Raphson算法，但其优势是无须用户提供函数的导函数，所以该算法是有意义的。 >> f=@(x)exp(-0.2*x).*sin(3*x+2)+cos(x); x=sec_sols(f,9,10,1e-15,1), x1=x(end), f(x1), z=zeros(size(x)); fplot(f,[8.5,13.5]), hold on, plot(x,z,'o'), hold off 图 2-10一元方程求解的中间过程 　　例 2-15试求解例 2-11中的二元方程，以 (1, 1)点为初值搜索到一个解。　　解由于同时含有自变量 x和 y，这样的方程是不能直接求解的，需要将其改写成 x的方程，最简单的方法是令 x1 = x，x2 = y，这样，原始的方程可以改写成 [ ] x21e.x1x 22/2 + e.x1/2 sin(x1x2) f(x1,x2)= 2 22 x2 cos(x1 + x2)+ x1ex1+x2 ·20·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 Jacobi矩阵不易用手工的方法推导出来，是需要通过解析推导求解的，所以应该在符号运算框架下输入原始函数，并通过 jacobian()函数计算出该矩阵。 >> syms x1 x2 f=[x1^2*exp(-x1*x2^2/2)+exp(-x1/2)*sin(x1*x2); x2^2*cos(x1+x2^2)+x1^2*exp(x1+x2)]; J=jacobian(f,[x1,x2]) 可以推导出函数的 Jacobi矩阵为 [ 2 2 22 2x1e.x1x2/2 . e.x1/2 sin (x1x2) /2 . x1x2e.x1x2/2/2 + x2e.x1/2 cos (x1x2) J = () 22 2ex1+x2 2x1ex1+x2 . xsin x2 + x1+ x1 2 ] 2 x1e.x1/2 cos (x1x2) . x13x2e.x1x2/2 () () 2 322ex1+x2 2x2 cos x2 + x1. 2x2 sin x2 + x1+ x1 有了原函数与 Jacobi矩阵，则可以手工写出这两个函数的匿名函数，然后调用基于 Newton–Raphson算法的求解函数。 >> f=@(x)[x(1)^2*exp(-x(1)*x(2)^2/2)+exp(-x(1)/2)*sin(x(1)*x(2)); x(2)^2*cos(x(1)+x(2)^2)+x(1)^2*exp(x(1)+x(2))]; J=@(x)[2*x(1)*exp(-x(1)*x(2)^2/2)-... exp(-x(1)/2)*sin(x(1)*x(2))/2-... x(1)^2*x(2)^2*exp(-x(1)*x(2)^2/2)/2+... x(2)*exp(-x(1)/2)*cos(x(1)*x(2)),... x(1)*exp(-x(1)/2)*cos(x(1)*x(2))-... x(1)^3*x(2)*exp(-x(1)*x(2)^2)/2; 2*x(1)*exp(x(1)+x(2))-x(2)^2*sin(x(2)^2+x(1))+... x(1)^2*exp(x(1)+x(2)),... 2*x(2)*cos(x(2)^2+x(1))-2*x(2)^3*sin(x(2)^2+x(1))+... x(1)^2*exp(x(1)+x(2))]; x=nr_sols(f,J,[1; 1],1e-15,1), length(x), norm(f(x(:,end))) 由该初值点出发得出的方程的解为 x = [5.1236, .12.2632]T，中间点的个数为 18。代入原方程后得出的误差矩阵范数为 3.9323×10.12。从求解的结果看，确实通过该函数可以求出方程的一个根，不过从实际操作看，这样的方法未免过于复杂，由于需要推导 Jacobi矩阵，并将其矩阵用匿名函数的形式手工表示出来，该过程比较容易出错。 若采用正割方法，则可以给出下面的语句，得出方程的根为 x = [1.0825, .1.1737]T，代入方程得出误差为 2.2841×10.14。 >> x=sec_sols(f,[1;1],[1.5; -3],1e-14,1) length(x), norm(f(x(:,end))) 虽然这时函数调用无须用户提供 Jacobi矩阵，但所需的迭代步数为 265，求解效率较低，所以对代数方程求解问题而言，需要更好的方法。 2.3.2 MATLAB的直接求解函数 由于上面给出的求解算法比较麻烦，需要提供的信息又比较难获取，所以在实际方程求解中应该考虑采用更好的方法。 MATLAB提供了更实用的求解函数 fsolve()，无须提供 Jacobi矩阵的句柄，只需给出方程函数的句柄和初值，则可以直接求解任意复杂的非线性方程组，由给出的初值搜索出方程的一个根。该函数的调用格式为 x=fsolve(f,x0)式中， f为方程函数的句柄； x0为初始向量或矩阵； f函数的维数与 x0完全一致，正常情况下得出的 x为方程的数值解。 该函数可以至多返回四个变量，这时完整的调用格式为 [x,F ,flag,out]=fsolve(f,x0,opts)式中， x为方程的解； F为 x处方程函数的值矩阵； flag如果为正说明求解成功， out变量还将返回一些中间信息。用户还可以增加输入选项 opts控制求解的算法与精度，后面将通过例子演示。 　　例 2-16试利用这里介绍的求解函数直接求解例 2-15中的方程。　　解仍然可以使用匿名函数描述原始的方程，且无须提供 Jacobi矩阵的函数句柄，直接调用求解函数即可以得出方程的解为 x1 = [0, 2.1708]T，将解代入方程则得出误差向量的范数为 3.2618×10.14。虽然这样得出的解不是例 2-15中得出的解，但也是方程的一个解。此外，还可以看出，迭代步数为 14次，f函数的调用次数为 38，与例 2-15中的结果相仿。不过该函数的优势是无须用户提供方程的导函数，使用该函数更适合于实际应用，建议采用该方法直接求解方程。 >> f=@(x)[x(1)^2*exp(-x(1)*x(2)^2/2)+exp(-x(1)/2)*sin(x(1)*x(2)); x(2)^2*cos(x(1)+x(2)^2)+x(1)^2*exp(x(1)+x(2))]; x0=[1; 1]; [x1,f1,key,cc]=fsolve(f,x0) 如果将初值修改为 x0 = [2, 1]T，则搜索到的解为 x1 =[.0.7038, 1.6617]T，该解对应的误差为 f1 =2.0242×10.12。 >> x0=[2; 1]; x1=fsolve(f,x0), f1=norm(f(x1)) 与前面介绍的 Newton–Raphson迭代法相比，求解方程的过程已经极大地简化了，由于只需描述方程本身，无须描述更复杂的 Jacobi矩阵，使得方程的求解变得轻而易举，所以在实际方程求解问题中可以放心使用这样的方法。 在前面的例子中，未知自变量 x与方程函数 f(x)都是同维数的向量，编写出匿名函数就可以描述方程函数，就可以求解方程从而得出方程的解。如果方程 f(x)= 0中，x与 f为同维数的矩阵，也可以直接利用 fsolve()函数搜索出方程 ·22·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 的数值解。下面以 Riccati方程为例介绍矩阵型方程的求解方法。　　定义 2-8 Riccati代数方程的数学形式为 ATX + XA . XBX + C = 0（2-3-6）式中，每个矩阵都是 n × n矩阵。 Riccati方程是以意大利数学家 Jacopo Francesco Riccati（1676.1754）命名的，本意是对应于一类一阶微分方程，其原型要求 B为正定矩阵，C为对称矩阵。后来因为微分方程难以求解，所以将其简化成上面的 Riccati代数方程。 从数学角度看，这几个矩阵可以为任意矩阵。在控制科学领域，可以考虑采用控制系统工具箱提供的 are()函数直接求取方程的数值解，不过该函数只能求出 Riccati方程的一个根。如果想获得方程全部的根，则可以考虑调用 vpasolve()函数直接求解。下面将给出例子演示 Riccati方程与多项式矩阵方程的方法。 　　例 2-17 试求解下面的 Riccati代数方程。 . . . .1 1 1 2 1 1 0 .2 .3 A = . 1 0 2 . , B = . .1 1 .1 . , C = . 1 3 3 . .1 .1 .3 .1 .1 0 .2 .2 .1 　　解可以先输入这几个矩阵，然后调用控制系统工具箱中的 are()函数求解 Riccati方程，并计算得出的误差矩阵范数。 >> A=[-1,1,1; 1,0,2; -1,-1,-3]; B=[2,1,1; -1,1,-1; -1,-1,0]; C=[0,-2,-3; 1,3,3; -2,-2,-1]; X=are(A,B,C) %解方程 norm(A'*X+X*A-X*B*X+C) %将得出的解代入原方程求误差矩阵的范数 由控制系统工具箱的 are()函数可以直接得出方程的一个根如下，代入原方程可见，该解导致的误差为 1.3980×10.14，说明该解的精度是比较高的。 . 0.21411546 .0.30404517 .0.57431474 X = . 0.83601813 1.60743230 1.397651600 ..0.004386346 0.20982900 0.24656718 由于该方程是二次型方程，人们很自然就可以想到，该方程可能有多个根，如何求出其他的根呢？显然可以尝试选择一个不同的初始矩阵，例如幺矩阵，从该矩阵求解方程的根。 >> f=@(x)A.'*x+x*A-x*B*x+C; x0=ones(3); %尝试另一个初值重新求解 [X1,f1,flag,cc]=fsolve(f,x0), norm(f1) 得出方程的解如下，对应的误差为 6.3513×10.9。 . 2.1509892 2.9806867 2.4175971 X1 = . .0.9005114 .1.3375371 .1.2847861 .0.95594506 1.8384489 1.7300025 显然，这时得出的解与 are()函数得出的解是不一致的。该函数返回的 f1矩阵就 是方程解处的函数值，即 f(X1)。另外，在这个例子中返回的 flag值为 1，因为它为正数，所以表示求解成功。另外返回的 cc信息包括如下内容，表明迭代步数为 11，函数调用次数为 102，说明求解效率还是很高的。 iterations: 11 funcCount: 102 algorithm: 'trust-region-dogleg' firstorderopt: 4.9891e-08 message: Equation solved. 2.3.3求解精度的设置 从上面例子得出的解看，误差偏大。有没有什么办法控制误差的大小呢？ 前面在定理 2-3中描述了一般迭代方法收敛的条件，给出了两个误差限 ξ1和 ξ2，这样的参数是可以人为设定的。可以由 opts=optimset命令得出方程求解的控制选项 opts，该变量是一个结构体型的数据，其中很多成员变量是可以修改的，例如，AbsTol是绝对误差限，与前面介绍的常数 ξ1是一致的，更常用的是相对误差限 RelTol，另外，函数值误差限 FunTol与常数 ξ2是一致的，所以可以通过设置这些误差限来控制求解的精度。 除了这两个选项之外，迭代步数 MaxIter和方程函数调用次数 MaxFunEvals也经常需要重新设置，否则，在求解方程过程完成之后可能会给出提示，指明求解次数超限。在这种情况下，一方面可以将这两个选项设置为更大的值，另一方面，也可以将得出的结果作为初值，重新调用 fsolve()函数继续求解，直至找出方程的解为止。这样的方法还可以与循环结构配合使用。 下面将通过例子演示求解精度的设定与控制方法。 　　例 2-18试选择幺矩阵作初值，重新求解例 2-17中的 Riccati代数方程。　　解如果重新设置求解精度控制选项 ff，在调用语句中可以直接使用该控制变量，得出更精确的解，这时，误差矩阵的范数为 2.5020×10.15，比默认精度有了明显的提高。 >> A=[-1,1,1; 1,0,2; -1,-1,-3]; B=[2,1,1; -1,1,-1; -1,-1,0]; C=[0,-2,-3; 1,3,3; -2,-2,-1]; f=@(x)A.'*x+x*A-x*B*x+C; ff=optimset; ff.TolX=eps; ff.TolFun=1e-20; x0=ones(3); X1=fsolve(f,x0,ff), norm(f(X1)) 　　例 2-19 试求解下面的 Riccati代数方程。 . . . .2 1 .3 2 1 .1 5 .4 4 A = . .1 0 .2 . , B = . 1 2 0 ., C = . 1 0 4 . 0 .1 .2 .1 0 .4 1 .1 5 　　解由下面的语句可以尝试求解 Riccati方程。 ·24·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 >> A=[-2,1,-3; -1,0,-2; 0,-1,-2]; B=[2,1,-1; 1,2,0; -1,0,-4]; C=[5 -4 4; 1 0 4; 1 -1 5]; X=are(A,B,C) 不过求解过程中会得出 “No solution: (A,B) may be uncontrollable or no solution exists”（(A,B)可能不可控或方程无解）错误信息。尽管 are()函数失效，仍然可以尝试求解原方程的数值解方法。例如，可以由下面的语句直接求解。 >> f=@(x)A.'*x+x*A-x*B*x+C; x0=-ones(3); [X1,f1,flag,cc]=fsolve(f,x0), norm(f(X1)) 得出的一个解如下，该解的误差范数为 1.2515×10.14。 . 5.5119052 .4.2335285 0.18011933 X1 = . 3.8739991 .3.6090612 .0.36369470 . 6.6867568 .4.9302000 1.07355930 2.3.4方程的复域求解 使用 fsolve()的另一个优势是，当初值选作复数时，有可能得出方程的复数根。另外，该函数还可以求取复数系数方程的数值解。 　　例 2-20试求例 2-17方程的复数根。　　解如果选择复数初值，则也有可能得出方程的复数根，同时可以验证，该根的共轭复数矩阵也满足原始方程。 >> A=[-1,1,1; 1,0,2; -1,-1,-3]; B=[2,1,1; -1,1,-1; -1,-1,0]; C=[0,-2,-3; 1,3,3; -2,-2,-1]; f=@(x)A.'*x+x*A-x*B*x+C; ff=optimset; ff.TolX=1e-15; ff.TolFun=1e-20; x0=eye(3)+ones(3)*1j; X1=fsolve(f,x0,ff), norm(f(X1)) X2=conj(X1), norm(f(X2)) 由该初值可以搜索出的解如下，相应的误差为 1.1928×10.14。对这个例子还可以同步求出方程根的共轭矩阵，经检验该矩阵也满足于原方程。 . 1.0979 + 2.6874j 1.1947 + 4.5576j 0.7909 + 4.1513j X1 = . .3.5784 + 1.3112j .5.8789 + 2.2236j .5.4213 + 2.0254j ..4.7771 . 0.4365j .7.8841 . 0.7403j .7.1258 . 0.6743j 对这个具体问题而言，如果初始值选择成实部为幺矩阵，虚部为单位阵的矩阵，则搜索出的将是实数解，这也说明由复数初始搜索点出发也能搜索出实数根。不过值得注意的是，这样搜索出的结果可能带有微小的虚部，其范数在该例子中为 6.4366×10.19，应该由 real()函数提取出来。 >> x0=ones(3)+eye(3)*1i; X1=fsolve(f,x0,ff), norm(f(X1)) norm(imag(X1)), X2=real(X1), norm(f(X2)) 　　例 2-21如果 Riccati代数方程的系数矩阵 A变成复数矩阵，试重新求解该方程。 . .1+8j 1+ j 1+6j A = . 1+3j 5j 2+7j . .1+4j .1+9j .3+2j 　　解可以将各个矩阵输入到 MATLAB的工作空间，然后使用匿名函数描述原方程。注意，原方程使用的是矩阵 A的直接转置 AT，不是 Hermite转置 AH，所以在匿名函数中应该使用 A.'，而不能使用 A'，否则方程描述是错误的，求解就没有意义了。可以使用下面命令直接求解并检验原方程，并求出解的共轭复数矩阵。 >> A=[-1+8j,1+1j,1+6j; 1+3j,5j,2+7j; -1+4j,-1+9j,-3+2j]; B=[2,1,1; -1,1,-1; -1,-1,0]; C=[0,-2,-3; 1,3,3; -2,-2,-1]; f=@(x)A.'*x+x*A-x*B*x+C; X0=ones(3); X1=fsolve(f,X0), norm(f(X1)), X2=conj(X1), norm(f(X2)) 这样可以得出方程的解为 . 0.0727 + 0.1015j 0.2811 . 0.1621j .0.3475 . 0.0273j X1 = . .0.0103 + 0.0259j .0.2078 + 0.2136j 0.2940 + 0.2208j ..0.0853 . 0.1432j .0.0877 . 0.0997j .0.0161 . 0.1826j 如果将解代入原方程，则得出误差矩阵的范数为 5.4565×10.13。对这个特定的例子而言，即使使用实数起始搜索矩阵，也能搜索出方程的复数解。此外，虽然能直接求出解矩阵的共轭复数矩阵，但显然该矩阵不满足原方程，这与实数矩阵方程不同，因为在实数矩阵方程中，如果找到了方程的一个解，其共轭矩阵一般会满足原方程。 如果采用 MATLAB控制系统工具箱的 are()函数求解方程的解： >> X=are(A,B,C), norm(f(X))得出矩阵方程的解如下，在求解过程中也没有任何警告或错误信息，不过试图将该解代入原方程后得出的误差过大，为 10.5210，说明得出的解根本不满足原始方程。 . 0.512140 .0.165860 0.020736 X = . .0.390440 0.322260 .0.051222 . .0.093655 .0.075522 0.110550 2.4联立方程组的精确求解 前面介绍过，利用图解方法只能求出给定方程的实数根，并不能求出方程的复数根，具体例子可以参见例 2-12。另外，如果联立方程有多个实数根，则只能用图形方法绘制出根所在的位置，并不能直接得出根的具体值，需要逐个根进行局部放大求解，求解过程比较烦琐。此外，前面介绍的数值求解方法由于每次只能求出方程的一个根，使用起来有时也不方便。 MATLAB工具箱提供了代数方程的解析求解函数 solve()，可以直接用于求解解析解存在的代数方程。如果方程的解析解不存在，则可以采用 vpasolve()函数求取方程的高精度数值解，解的误差可能达到 10.30或更高精度，远大于双精度数据结构下的数值解。本书称这类解为准解析解，以区别于方程的解析解与双精度意义下的数值解。 ·26·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 2.4.1低阶多项式方程的解析求解 一元一次和一元二次方程可以利用 solve()函数直接求解，该函数还可以用于含有其他参数的方程求解。不过如果有其他参数的存在，则可能利用现有的 solve()函数不易得出三次或四次方程的解析解，尽管这些方程是存在代数解法的，除非给出特别的控制选项，后面将通过例子演示。 solve()函数的调用格式为 S=solve(eqn1,eqn2,··· ,eqnn), %求解方程 S=solve(eqn1,eqn2,··· ,eqnn,x1,x2,··· ,xn), %指定未知变量求解方程式中，待求解的方程由符号表达式 eqni表示，自变量由 xi表示。返回的解是一个结构体型变量，其解由 S.xi直接提取。在调用格式中 eqni可以是单个的方程也可以是向量、矩阵描述的一组方程，还可以将所有的方程描述成一个向量与矩阵符号表达式 eqn1，直接求解这些方程。 当然，用下面的调用格式还可以直接获得方程的解。 [x1,x2,··· ,xn]=solve(...), %输入变量表示与前面一致 从函数的调用和使用方面看，这种直接返回变量的调用格式比返回结构体变量的格式更实用，所以本书尽量使用这样的格式。 　　例 2-22试重新求解例 2-2中的鸡兔同笼问题。　　解声明符号变量，并将方程用符号表达式表示出来，则可以调用 solve()函数直接求解给定的方程。注意，方程中的等号应该由双等号表示。 >> syms x1 x2 [x1,x2]=solve(x1+x2==35,2*x1+4*x2==94) 得出方程的解为 x1 = 23，x2 = 12。返回变量的名字还可以选择为其他的变量，此外，如果方程右边为零，则可以省略等号。例如，上面的求解语句还可以修改成 >> syms x1 x2; [x0,y0]=solve(x1+x2-35,2*x1+4*x2-94) 　　例 2-23中国唐代数学家王孝通所著《缉古算经》中有一道应用题，翻译成现代数学符号可以写出下面的三次方程，试求解该方程。 11 x 3 + 5004x 2 + 1169953 x = 41107188 33 　　解利用符号运算方法可以直接求解这个三次方程。 >> syms x solve(x^3+5004*x^2+(1169953+sym(1/3))*x==(41107188+sym(1/3))) 得出方程的三个解如下。王孝通只得到了 31这个解，由于它为整数，是原应用题的解，其他两个是负数，不是应用题的解。 ∞∞∞ ∞ 5035 3 5 29 414767 31, .± 26 　　例 2-24试求解下面的联立方程组，并给出方程有实数根的条件。 { ax + cy =2 bx2 + cx + ay2 . 4xy = .3 　　解可以看出，如果对方程进行简单变换，则可以得到关于 x或 y一元二次方程，利用相应的求根公式，则可以写出方程的解析解。不过从给出的表达式看，如果想手工求解这样的方程并不是一件简单的事情，所以应该利用计算机来完成这样烦琐的推导。 先声明必要的符号变量，则可以将两个方程用符号表达式描述出来，再给出下面的直接求解命令。 >> syms x y a b c f1=x*a+c*y-2; f2=b*x^2+c*x+a*y^2-4*x*y+3; [x0,y0]=solve(f1,f2) 则可以得出方程的一对解为 (∞ )c 8a +4bc + ac2 ± ac4 .48ac.8a2c.12bc2 .12a3 .16c2 .16ab+642 x0 = . . 2a (a3 +4ac + bc2) a ∞ 8a +4bc + ac2 ± ac4 .48ac.8a2c.12bc2 .12a3 .16c2 .16ab+64 y0 = 2(a3 +4ac + bc2) 由得出的结果看，如果期望方程有实数根，则根号下的表达式应该大于等于零，即 2 c 4 . 48ac . 8ac . 12bc2 . 12a 3 . 16c 2 . 16ab + 64 . 0 　　例 2-25试在符号运算的框架下求取例 2-7中复系数多项式方程的解析解。　　解由于已知方程的解析解存在，所以可以调用 solve()函数解方程，最终得出方程的解析解为 [.2, .1 . j, .1 . j, .1 . j]T。 >> syms s; f=s^4+(5+3i)*s^3+(6+12i)*s^2-(2-14i)*s-(4-4i); s0=solve(f) 　　例 2-26试求解四次方程的解析解。 7s 4 + 119s 3 + 756s 2 + 2128s + 2240 = 0 　　解由于四次方程有代数解求根公式，所以可以将其输入 MATLAB环境，则可以调用 solve()函数尝试求解。 >> syms s; f=7*s^4+119*s^3+756*s^2+2128*s+2240; x=solve(f) 得出方程的解为 x = .5，.4，.4，.4。事实上，虽然四次方程有自己的代数解求根公式，得出的解也可能是无理数，换句话说，真正意义下的解析解也可能不存在。例如，如果原方程左侧加 1，则方程仍然是四次多项式方程，这时方程的解析解可能不存在。 >> f=f+1; x=solve(f)这样的方程由 solve()函数是不能直接求解的，因为直接得出的结果如下： ·28·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 root(z^4 + 17*z^3 + 108*z^2 + 304*z + 2241/7, z, 1) root(z^4 + 17*z^3 + 108*z^2 + 304*z + 2241/7, z, 2) root(z^4 + 17*z^3 + 108*z^2 + 304*z + 2241/7, z, 3) root(z^4 + 17*z^3 + 108*z^2 + 304*z + 2241/7, z, 4) 表明方程的解是 z4 + 17z3 + 108z2 + 304z + 2241/7 = 0的多项式方程的四个根。如果实在想得出方程的数值解，则可以调用 vpa()函数或直接调用 vpasolve()函数直接求解方程，得出的误差向量范数为 4.7877×10.36。 >> x1=vpa(x), norm(subs(f,s,x1)) 其实，若真的想求出三次或四次方程的解析解，还是可以通过设定'MaxDegree'选项实现的，不过得出的解很冗长。下面将通过例子演示得出的解。 　　例 2-27试求出例 2-26改变后四次方程的解析解。　　解可以重新求解该方程，由于是四次方程，所以将'MaxDegree'选项设置为 4。 >> syms s; f=7*s^4+119*s^3+756*s^2+2128*s+2241; x=solve(f,'MaxDegree',4), x(1) 得出的结果还是很烦琐的，经过细心的手工变量替换与化简可得 √ 2 22 2 9r2r/2 . 9r4r. 12r/7 + 27r3/4 r22217 x1 = . 2 . . 6r 6rr2 4 式中， √ ∞∞∞ √ 1 3767i 9r2 12 4 r = 6 + ,r2 = +9r4 + 14882 4 7 有了这样的求解公式，则可以得出任意精度的方程的解。例如，由下面的语句可以得出误差为 7.9004×10.104。 >> norm(vpa(subs(f,s,x),100)) 2.4.2多项式型方程的准解析解 对于一般的多项式型代数方程而言，采用 vpasolve()函数有望得出方程全部的准解析解，而一般非线性代数方程只能得到一个准解析解。本节将通过例子演示各类方程的准解析解方法。 　　例 2-28试求解例 2-12中给出的联立方程。　　解由图解法并不能有效求解例子中的联立方程，因为原方程既含有实数根，也含有复数根。可以先将方程用符号表达式的方式输入给计算机，然后调用 vpasolve()函数，则可以直接求解该方程组。 >> syms x y; [x0,y0]=vpasolve(x^2+y^2==5,x+4*y^3+3*y^2==2) ..... ..... 得出的结果在 和向量中返回，所以要检验得出的误差并不是一件容易的事，xy00因为要重新输入方程的表达式，再求值。另一种求解的方法是将两个方程用符号变量表示出来，再直接求解，这样得出的结果与前面是完全一致的。将解代入原方程，可以得出f1=x^2+y^2-5;f2=x+4*y^3+3*y^2-2; >> norm([subs(f1,{x,y},{x0,y0}),subs(f2,{x,y},{x0,y0})])从这里给出的例子看，求解这样复杂的高阶代数方程组的难度对用户而言，与求解鸡兔同笼问题是一样的，用户需要做的就是将方程用符号表达式表示出来送给求解函vpasolve()更简单地，还可以用向量的形式表示两个方程，这样，再调用函数则可 得出的解为 .2.0844946216518881587217059536735 2.25174086438823135598559542792.0.00566796162874471821186269437616j x0 = .2.2559104938211695667214623634±0.2906418591129883375435021776769j 2.0928338805177645801934398246011 0.80924790534443253378482794183867 0.0473834464122358325607839974582±0.26935104521931155183432307891j y0 = .0.808292289977322744563366509227.0.81116945942300081793886298527j .0.7874302182142587097796629183005 误差为 1.0196 × 10.38。 [x0,y0]=vpasolve(f1,f2) 数，然后等待得出的解就行了。 以重新求解相应的方程，得出的结果与前面语句完全一致。 >> F=[x^2+y^2-5; x+4*y^3+3*y^2-2]; [x0,y0]=vpasolve(F) norm(subs(F,{x,y},{x0,y0})) 　　例 2-29试求解下面看起来很复杂的代数方程。 . ..... . ..... . .. .. 1 325 3 2 x + x +++ +=0 22 y 2y2 x3 y 31 + + +5y 4 =0 22xx4 　　解这个方程与例 2-12中给出的方程不一样，至少例 2-12中的方程可以将一个方程代入另一个方程，最终手工得出一个高阶的多项式方程，而这个方程就没那么容易变换了。不用说求解方程，就是判定方程有多少个根，不借助计算机也不是一件容易的事。 不过不必考虑或担心这类底层问题，只需用下面的语句将原始的方程用规范的语句原原本本表示出来，就可以直接得出原方程的准解析解。 >> syms x y; f1(x,y)=x^2/2+x+3/2+2/y+5/(2*y^2)+3/x^3; f2(x,y)=y/2+3/(2*x)+1/x^4+5*y^4; ·30·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 [x0,y0]=vpasolve(f1,f2), size(x0) %解方程并得到解的个数 e1=norm(f1(x0,y0)), e2=norm(f1(x0,y0)) %检验误差将得出的全部 26个根代入原始方程，则能得出很小的计算误差，达到 10.33级，说明该方程的各个解都是非常精确的。求解这样看起来难度极高的代数方程，对用户而言，求解的难度也与鸡兔同笼问题是一样的。 即使著名的 Abel–Ruffini定理已经指明这类高阶多项式方程没有解析解，还是可以通过代数方程求解方法得出高精度的准解析解，得出解的精度是非常高的。 2.4.3高次多项式矩阵方程的准解析解 定义 2-8中描述了代数 Riccati方程，该方程是关于 X的二次型方程，是多项式型方程。如何用符号表达式表示 Riccati方程是求解该方程的关键一步。这里先通过例子探讨使用 vpasolve()函数求解该方程的方法。 　　例 2-30试求解例 2-17中给出的 Riccati代数方程全部的根。　　解如果想得出方程全部的根，则应该尝试 vpasolve()函数。若想表示未知的 X矩阵，则需要将其设置为符号变量，再构成符号矩阵，这样就可以由简单的符号表达式描述 Riccati方程本身，再调用 vpasolve()函数求解方程。 >> A=[-1,1,1; 1,0,2; -1,-1,-3]; B=[2,1,1; -1,1,-1; -1,-1,0]; C=[0,-2,-3; 1,3,3; -2,-2,-1]; X=sym('x%d%d',3); F=A.'*X+X*A-X*B*X+C tic, X0=vpasolve(F), toc .由符号运算可以推导出 Riccati方程对应的联立方程为 x12 .x13 .x31(1+x11 .x12)+x11(x12 .2x11 +x13 .2).x21(1+x11 +x12 .x13)=0 x11 .x13 .x32(1+x11 .x12)+x12(x12 .2x11 +x13 .1).x22(x11 +x12 .x13 .1)=2 x11 +2x12 .x33(1+x11 .x12)+x13(x12 .2x11 +x13 .4).x23(x11 +x12 .x13 .1) = 3 . x22 .x23 .x31(1+x21 .x22)+x11(1+x22 .2x21 +x23).x21(1+x21 +x22 .x23)= .1 . x21 .x23 .x32(1+x21 .x22)+x12(x22 .2x21 +x23).x22(x21 +x22 .x23)=.3 x21 +2x22 .x33(1+x21 .x22)+x13(1+x22 .2x21 +x23).x23(3+x21 +x22 .x23)=.3 x32 .x33 .x31(4+x31 .x32)+x11(1+x32 .2x31 +x33).x21(x31 +x32 .x33 .2)=2 x31 .x33 .x32(3+x31 .x32)+x12(1+x32 .2x31 +x33).x22(x31 +x32 .x33 .2)=2 . . x31 +2x32 .x33(6+x31 .x32)+x13(1+x32 .2x31 +x33).x23(x31 +x32 .x33 .2)=1 经过 23.04s的等待，可以得出方程全部的 20组根，其中，第 5、10、15、18这四组根为实数矩阵，其余为共轭复数矩阵。得出的四个实数矩阵分别为 . 1.9062670985148 2.6695228037028 4.1090269897993 X5 = . .4.3719461205750 .3.2277001659457 .5.7232367559600 . .8.1493167800300 .2.8535676463847 .12.90505951546 . 8.69508738924215 .8.369677652147 15.08579547886 X10 = . .15.265261644743 14.485759397031 .23.336518219659 . .20.7167885668890 17.582212700878 .33.22526573531 . 0.21411545933325 .0.3040451651414 .0.5743147449581 X15 = . 0.83601813100313 1.60743227422054 1.397651628726543 . .0.0043863464229 0.209828998159396 0.246567175609337 . 2.1509892346834 2.980686747114 2.4175971297531 X18 = . .0.900511398366 .1.3375371059663 .1.284786114934488 .0.95594505835086 1.83844891740592 1.7300024774690319 例 2-17中得出的实数根是原方程的第 15个根，可以提取该根，代入方程检验，得出 的误差为 1.8441×10.39。 >> k=15; %选择提取根的序号 X1=[X0.x11(k) X0.x12(k) X0.x13(k); X0.x21(k) X0.x22(k) X0.x23(k); X0.x31(k) X0.x32(k) X0.x33(k)]; norm(A.'*X1+X1*A-X1*B*X1+C) 　　例 2-31试用准解析解方法求解例 2-19中方程的全部根。 　　解给出下面的命令即可求出方程的全部 20个根，其中有 8个实根，其余为共轭复 数根，求解过程耗时为 36.75s。 >> A=[-2,1,-3; -1,0,-2; 0,-1,-2]; B=[2,1,-1; 1,2,0; -1,0,-4]; C=[5 -4 4; 1 0 4; 1 -1 5]; X=sym('x%d%d',3); F=A.'*X+X*A-X*B*X+C tic, X0=vpasolve(F), toc 前面介绍了 Riccati方程的求解方法，如果将 Riccati方程中 AT项替换成一个 自由矩阵 D，则可以引出广义 Riccati方程。 　　定义 2-9广义 Riccati代数方程的数学形式为 DX + XA . XBX + C = 0 （2-4-1） 式中，各个矩阵都是 n × n矩阵。 广义 Riccati方程没有现有的 MATLAB函数直接求解，仍可尝试 vpasolve() 函数直接求解，得出方程全部的根。 　　例 2-32若已知如下矩阵，试求解广义 Riccati方程。 . ... .111 211 0 .2 .32 .1 .1 . .. .. ... A= 102 , B =.11 .1 , C =133 , D=1 1 .1 .1 .1 .3 .1 .10 .2 .2 .11 .1 　　解和以前一样，先输入几个矩阵，然后由符号表达式的方式描述方程，再给出求 解命令就可以直接求解方程了。 >> A=[-1,1,1; 1,0,2; -1,-1,-3]; B=[2,1,1; -1,1,-1; -1,-1,0]; C=[0,-2,-3; 1,3,3; -2,-2,-1]; D=[2,-1,-1; 1,1,-1; 1,-1,0]; X=sym('x%d%d',3); F=D*X+X*A-X*B*X+C tic, X0=vpasolve(F), toc ·32·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 经过 25s的等待，由上述求解语句可以直接得出方程的 20个根，其中， 8个根为实数矩阵，其余的为共轭复数矩阵。 　　例 2-33本丛书卷 III例 5-42曾探讨了一个高次矩阵方程的求解问题。 AX3 + X4D . X2BX + CX . I = 0 其中，已知的矩阵为 . . . . 2 1 9 0 3 6 7 0 3 3 9 5 A = . 9 7 9 . , B = . 8 2 0 ., C = . 5 6 4 . , D = . 1 2 9 . 6 5 3 8 2 8 1 4 4 3 3 0 试用这里给出的方法求取该方程的解矩阵。　　解很自然地，可以给出下面的求解命令，不幸的是，经过 1843.8s的长时间等待，只得出了方程的一个根，并给出警告信息 “Solutions might be lost”（根可能丢失了），说明用这样的方法也不能确保得出方程所有的根。 >> A=[2 1 9; 9 7 9; 6 5 3]; B=[0 3 6; 8 2 0; 8 2 8]; C=[7 0 3; 5 6 4; 1 4 4]; D=[3 9 5; 1 2 9; 3 3 0]; X=sym('x%d%d',3); f=A*X^3+X^4*D-X^2*B*X+C*X-eye(3); tic, X0=vpasolve(f), toc 　　例 2-34试求解含有复数矩阵的 Riccati方程的准解析解。　　解很自然地可以考虑使用下面的语句直接求解需要的复系数 Riccati方程，同样的步骤和语句，由于涉及复系数，求解过程极其耗时，大约 379s才能得出结果（实系数方程耗时 23.04s，相差 16倍）。可以得出方程的 20个根，全部为复数根。 >> A=[-1+8i,1+1i,1+6i; 1+3i,5i,2+7i; -1+4i,-1+9i,-3+2i]; B=[2,1,1; -1,1,-1; -1,-1,0]; C=[0,-2,-3; 1,3,3; -2,-2,-1]; X=sym('x%d%d',3); F=A.'*X+X*A-X*B*X+C tic, X0=vpasolve(F), toc 2.4.4非线性代数方程的准解析解 如果用符号表达式可以描述出非线性联立方程组，则可以由 vpasolve()函数直接求解方程，得出方程的解。与多项式类方程不同，这样得出的解很可能是众多解中的一个，如果原始方程含有多个解，则用户可以自己选择一个搜索的初值，从该初值 x0出发直接搜索准解析解。相应的函数调用格式为 x=vpasolve(eqn1,eqn2,··· ,eqn,x1,x2,··· ,xn,x0) n 式中， x0为未知量向量的初始搜索点。如果原方程是多项式型的联立方程，则 x0对求解结果没有影响。该函数的另一种调用格式为 x=vpasolve(eqn1,eqn2,··· ,eqn,x1,x2,··· ,xn,xm,xM) n 式中，xm与 xM为未知数向量 x的下界与上界向量。 　　例 2-35试求解例 2-11给出的代数方程的根，可以选择初始搜索点 (2, 2)。　　解利用下面自然的方式可以直接求解非线性代数方程。 >> syms x1 x2 f=[x1^2*exp(-x1*x2^2/2)+exp(-x1/2)*sin(x1*x2); x2^2*cos(x1+x2^2)+x1^2*exp(x1+x2)]; [x0,y0]=vpasolve(f,[2; 2]), norm(subs(f,{x1,x2},{x0,y0})) 得出方程的解如下，代入原方程则误差的范数为 2.0755×10.37。由此可见，方程解的精度远高于前面介绍的数值解法精度。 x0 =3.0029898235869693047992458712192 y0 = .6.2769296697194948789764344182923 从给出的例子可见，尽管这里给出的方法能得出方程的高精度解，但仍然有一个问题尚未解决，就是如何一次性地求出如图 2-6所示的所有交点坐标，即联立方程在感兴趣区域内所有的解，这将是 2.5节需要解决的问题。 2.5多解矩阵方程的求解 尽管前面介绍的 vpasolve()函数能够一次性求出某些方程全部的解，但对一般矩阵方程，尤其是非线性矩阵方程而言是无能为力的，也没有其他 MATLAB程序能够求解任意的非线性矩阵方程。 前面介绍了由给定初值求解函数的几种方法，但这些方法很难一次性求解多解非线性方程所有的解，所以应该构建更简单的函数完成这样的任务，本节给出一种求解的思路，并依照该思路编写出 MATLAB通用程序，试图得出方程感兴趣区域内全部的解，并再扩展一步该方法，试图得出方程的全部准解析解。 2.5.1方程求解思路与一般求解函数 由前面给出的求解函数可见，如果选定了一个初值，则可以通过随机数矩阵生成的方式产生一个初始搜索点，得出方程的解。更一般地，可以建立一个循环结构实现这样的操作。由初始搜索点，调用 fsolve()函数得出方程的一个解，如果这个解已经被记录，则可以比较这个解和已记录解的精度，如果新的解更精确，则用这个解取代已记录的解，否则舍弃。如果这个解是新的解，则记录该解。 如果这个循环结构设计成死循环，则有望得出方程全部的解。根据这样的思路可以编写出一个通用的求解函数，其实这个函数以前曾经发布了多个版本，这个版本中，特别增加了一些处理，例如，可以尝试零矩阵是不是方程的孤立解；如果找到的根比以前存储的更精确，则替换该根；如果找到的新解为复数，则检验其共轭复数是不是方程的根。基于这样的考虑，编写了下面的求解函数。这个函数有一个特 ·34·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 点，运行的时间越长，得出的结果可能越精确。 function more_sols(f,X0,varargin) [A,tol,tlim,ff]=default_vals({1000,eps,30,optimset},varargin{:}); if length(A)==1, a=-0.5*A; b=0.5*A; else, a=A(1); b=A(2); end ar=real(a); br=real(b); ai=imag(a); bi=imag(b); ff.Display='off'; ff.TolX=tol; ff.TolFun=1e-20; [n,m,i]=size(X0); X=X0; tic if i==0, X0=zeros(n,m); %判定零矩阵是不是方程的孤立解 if norm(f(X0))1e-5, x0=x0+(ai+(bi-ai)*rand(n,m))*1i; end [x,aa,key]=fsolve(f,x0,ff); t=toc; if t>tlim, break; end if key>0, N=size(X,3); %读出已记录根的个数，若找到的根已记录，则放弃 for j=1:N, if norm(X(:,:,j)-x)<1e-4; key=0; break; end, end if key==0 %如果找到的解比存储的更精确，则替换 if norm(f(x))0, X(:,:,i+1)=x; %记录找到的根 if norm(imag(x))>1e-5 && norm(f(conj(x)))<1e-8 i=i+1; X(:,:,i+1)=conj(x); %若找到复根，则测试其共轭复数 end assignin('base','X',X); i=i+1, tic %更新信息 end, assignin('base','X',X); end, end该函数调用底层公用函数 default_vals()，其内容在其他卷中有过描述，为 方便起见，这里重新给出。 function varargout=default_vals(vals,varargin) %读默认参数通用子函数 if nargout=length(vals), error('number of arguments mismatch'); else, nn=length(varargin)+1; %用循环结构指派各个默认参数 varargout=varargin; for i=nn:nargout, varargout{i}=vals{i}; end, end, end方程求解函数 more_sols()的调用格式为 more_sols(f,X0,a,.,tlim,opts)式中， f为方程的函数句柄，可以由匿名函数与 M函数描述原代数方程； X0为三维数组，用于描述解的初值，如果首次求解方程，建议将其设置为 zeros(n,m,0)，即空白三维数组， n和 m为解矩阵的维数；方程的解被自动存在 MATLAB工作空间中的三维数组 X中，如果想继续搜索方程的解，则应该在 X0的位置填写 X；a的 默认值为 1000，表示在 [.500, 500]区间内大范围搜索方程的解； .的默认值为 eps； tlim的默认值为 30，表示 30s没有找到新的解就自动终止程序；还可以指定求解的控制选项 opts，默认值为 optimset。a还可以取为复数，表示需要求取方程的复数根。另外，a还可以给定为求解区间 [a, b]。 　　例 2-36试重新求解例 2-10中的一元超越方程。　　解从图 2-5给出的曲线可见，该方程在期望的区域内有 6个交点，利用编写的 more_sols()函数可以求出这 6个交点，并在图 2-11上直接标注出来。得出方程的解为 x1 = [1.4720, 4.6349, 7.7990, 10.9601, 14.1159, 17.2666]。 >> f=@(x)exp(-0.2*x).*sin(3*x+2)-cos(x); more_sols(f,zeros(1,1,0),[0,20]) x0=X(:); x1=x0(x0>=0 & x0<=20), fplot(f,[0,20]), hold on plot(x1,f(x1),'o',[0,20],[0,0],'--'), hold off 图 2-11超越方程的全部实根 　　例 2-37试求解例 2-11给出的代数方程在 .2巾 . x1,x2 . 2巾区域内全部的根。　　解利用下面自然的方式可以直接求解非线性代数方程，先用匿名函数的形式描述联立方程，这样就可以调用 more_sols()函数求取方程在感兴趣区域内全部的数值解。这里将 A可以选作 13，比 4巾略大的值，得出的解个数多于在感兴趣区域内的解。 >> f=@(x)[x(1)^2*exp(-x(1)*x(2)^2/2)+exp(-x(1)/2)*sin(x(1)*x(2)); x(2)^2*cos(x(1)+x(2)^2)+x(1)^2*exp(x(1)+x(2))]; A=13; more_sols(f,zeros(2,1,0),A) 可以提取出感兴趣区域内全部的根，可以发现，总的根的个数为 110个，可以将得到的解叠印到图解法的图形上，如图 2-12所示。 >> ii=find(abs(X(1,1,:))<=2*pi && abs(X(2,1,:))<=2*pi); X1=X(:,:,ii); size(ii) %提取感兴趣区域内的根 x=X1(1,1,:); x=x(:); y=X1(2,1,:); y=y(:); plot(x,y,'o') syms x y; f1=x^2*exp(-x*y^2/2)+exp(-x/2)*sin(x*y); ·36·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 图 2-12方程的图解法与得出的结果 f2=y^2*cos(x+y^2)+x^2*exp(x+y); hold on; fimplicit([f1,f2],[-2*pi,2*pi],'Meshdensity',800) hold off, axis(2*pi*[-1,1,-1,1]) 　　例 2-38试用数值方法重新求解例 2-32中的广义 Riccati方程。　　解和以前一样，先输入几个矩阵，然后由匿名函数的方式描述方程，再给出求解命令就可以直接求解方程，得到方程的 8个实数根。 >> A=[-1,1,1; 1,0,2; -1,-1,-3]; B=[2,1,1; -1,1,-1; -1,-1,0]; C=[0,-2,-3; 1,3,3; -2,-2,-1]; D=[2,-1,-1; 1,1,-1; 1,-1,0]; F=@(X)D*X+X*A-X*B*X+C, more_sols(F,zeros(3,3,0)) 继续求解方程则可以得出方程全部的复数根，如果将 A设置成复数量。这样总共可以得到方程的 20个根，与例 2-32得出的结果完全一致。 >> more_sols(F,X,1000+1000i) 　　例 2-39如果例 2-32中的广义 Riccati方程变换成 DX + XA . BX2 + C = 0，试用不同的方法求解该方程。　　解如果用 more_sols()函数直接求解，则可以得到 19个实数根。 >> A=[-1,1,1; 1,0,2; -1,-1,-3]; B=[2,1,1; -1,1,-1; -1,-1,0]; C=[0,-2,-3; 1,3,3; -2,-2,-1]; D=[2,-1,-1; 1,1,-1; 1,-1,0]; F=@(X)D*X+X*A-B*X^2+C; more_sols(F,zeros(3,3,0)) 如果在复数范围内搜索方程的根，总共可以找到 57个根。现在尝试 vpasolve()函数的准解析解求解方法，经过 5873.8s的等待，可以得到 60个根。 >> X=sym('x%d%d',[3,3]); F=D*X+X*A-B*X^2+C, tic, X0=vpasolve(F), toc 两种方法得出的根的数目差不多，不过仔细对比得出的根，如 x11的值，可以发现，由后者得出的根有 11个位于无穷远处，其幅值大于 1010，而前者得出的都是在合理范围内的数值。两种方法得出的 x11对比在表 2-1中给出。此外，经检验， x11 = 38的根是第 21个根，将其提取并代入原方程，得出误差矩阵的范数为 18981334.363，不满足原方 程。x11 = .1.5146的根是第 15个根，也不满足原始方程。由此可见，采用 vpasolve()函数求解，即使处理多项式矩阵方程，得出的结果也是不可靠的。 >> length(find(abs(X0.x11)>1e10)), k=21; x0=[X0.x11(k) X0.x12(k) X0.x13(k); X0.x21(k) X0.x22(k) X0.x23(k); X0.x31(k) X0.x32(k) X0.x33(k)] norm(F(x0)) 表 2-1两种方法得到的 x11值对比 准解析解数值解 .3.006466262658×1063 . ..71.6679956 .1.5146484375 . .0.393316966940371 .0.3933169669 0.19217475249284 0.1921747525 1.10128812201771 1.101288122 1.83357133024232 1.83357133 3.1272528908906671 3.127252891 7.0244296789169 7.02442972 . 15.92032583 . 52.986452 79369526682545291264.0 . 5.7961230691484403×1059 . 准解析解数值解 ..431.5070709 .3.4447376962310 .3.444737696 ..1.450204069 0.1563285970559477 0.1563285971 0.925385242867366 0.9253852429 1.712893235210277 1.712893235 3.121670783661601 3.121670784 5.803147834137380 5.803147834 7.387972365540733 7.387972366 38.0 . . 577.0832185 169534813891083632640.0 . 4.692574987609119×1067 . 　　例 2-40试求解非线性矩阵方程。 eAX sin BX . CX + D = 0其中 A、B、C和 D矩阵在例 2-32中给出，试求出该方程的全部实根。　　解可以用下面的语句直接求解这里给出的复杂非线性矩阵方程，已经找到 122个实根。用户还可以自己尝试，看看能不能得出更多的实根（根的存储文件为 data2_36.mat）。 >> A=[2 1 9; 9 7 9; 6 5 3]; B=[0 3 6; 8 2 0; 8 2 8]; C=[7 0 3; 5 6 4; 1 4 4]; D=[3 9 5; 1 2 9; 3 3 0]; f=@(X)expm(A*X)*funm(B*X,@sin)-C*X+D; more_sols(f,zeros(3,3,0),10); X 2.5.2伪多项式方程的求解 伪多项式方程是非线性矩阵方程的一个特例，因为未知矩阵实际上是一个标量。这里将给出伪多项式方程的定义，并给出求解方法。 　　定义 2-10伪多项式（pseudo-polynomial）方程的一般数学形式为 p(s)= c1sα1 + c2sα2 + ··· + cn.1sαn.1 + cnsαn =0（2-5-1） ·38·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 式中，向i为实数。 可见，伪多项式方程是常规多项式方程的扩展，求解方法可能远难于普通多项式方程的求解方法。这里将探讨不同方法的可行性。 　　例 2-41试求解伪多项式方程 [5] x 2.3 +5x 1.6 +6x 1.3 . 5x 0.4 +7=0。　　解一种容易想到的方法是引入新变量 z = x0.1，这样原方程可以映射成关于 z的多项式方程如下： z 23 +5z 16 +6z 13 . 5z 4 +7=0 可以求出，该方程有 23个根，再用 x = z10就可以求出方程全部的根。这样的思想可以由下面的 MATLAB语句实现。 >> syms x z; f1=z^23+5*z^16+6*z^13-5*z^4+7; p=sym2poly(f1); r=roots(p); f=x^2.3+5*x^1.6+6*x^1.3-5*x^0.4+7; r1=r.^10, double(subs(f,x,r1)) 不过，把这样得出的解代回到原来的方程可以发现，绝大部分的根都不满足原有的伪多项式方程。原方程到底有多少个根呢？上面得到的 x只有两个根满足原方程，即 x = .0.1076 ± 0.5562j，其余的 21个根都是增根。由下面语句也可以得出同样两个根。 >> f=@(x)x.^2.3+5*x.^1.6+6*x.^1.3-5*x.^0.4+7; more_sols(f,zeros(1,1,0),100+100i), x0=X(:) 从数学角度看，这对真实的根位于第一 Riemann叶上。其余的解位于其他 Rie-mann叶上，都是方程的增根，但不满足原方程。 　　例 2-42试求解无理阶次伪多项式方程 s5 + 25s3 + 16s2 . 3s0.4 +7=0。 　　解由于这里的阶次是无理数，无法将其转换成前面介绍的普通多项式方程，所以 more_sols()函数就成了求解这类方程的唯一方法，可以由下面的语句直接求解该方程。可见，该无理阶伪多项式方程只有两个根，位置为 s = .0.0812 ± 0.2880j。 >> f=@(s)s^sqrt(5)+25*s^sqrt(3)+16*s^sqrt(2)-3*s^0.4+7; more_sols(f,zeros(1,1,0),100+100i); x0=X(:) err=norm(f(x0(1))) 将得出的解代回原方程，则可见误差为 9.1551×10.16，该解在数值意义下足够精确。从这个例子还可以看出，即使阶次变成了无理数，并未给求解过程增加任何麻烦，求解过程和计算复杂度与前面的例子完全一致。 2.5.3高精度求解函数 在这里给出的 more_sols()函数中，核心求解工具是 fsolve()函数，若将其替换为高精度的 vpasolve()函数，则可以编写出高精度的非线性函数求解程序。 function more_vpasols(f,X0,varargin) [A,tlim]=default_vals({1000,60},varargin{:}); X=X0; %读入默认参数 if length(A)==1, a=-0.5*A; b=0.5*A; else, a=A(1); b=A(2); end ar=real(a); br=real(b); ai=imag(a); bi=imag(b); [i,n]=size(X0); tic while (1), %死循环结构，可以按 Ctrl+C组合键中断，也可以等待 x0=ar+(br-ar)*rand(1,n); %生成初始随机实矩阵 if abs(imag(A))>1e-5, x0=x0+(ai+(bi-ai)*rand(1,n))*1i; end V=vpasolve(f,x0); N=size(X,1); key=1; x=sol2vec(V); %求解方程 if length(x)>0 %如果找到的解非空，则继续判定，否则直接放弃 t=toc; if t>tlim, break; end %若一段时间没找到新根，则终止整个程序 for j=1:N, if norm(X(j,:)-x)<1e-5; key=0; break; end, end if key>0, i=i+1; X=[X; x]; %若找到新根，则记录该根 disp(['i=',int2str(i)]); assignin('base','X',X); tic end, end, end function v=sol2vec(A) %子函数，将根转换成行向量 v=[]; A=struct2cell(A); for i=1:length(A), v=[v, Ai]; end %转换成行向量 该函数的调用格式为 more_vpasols(f,X0,A,tlim)，其中还嵌入了为其设计的底层支持子函数 sol2vec()，将得出的解转换成行向量。输入变量 f可以为符号型的行向量来描述联立方程，初始矩阵 X0指定为 zeros(0,n)，其中， n为未知数的个数。其他的输入变元与前面介绍的 more_sols()函数是一致的。返回的变量 X(i,:)存储找到的第 i个解。值得指出的是， more_vdpsols()函数的运行速度比 more_sols()函数慢得多，但精度也高得多。 　　例 2-43考虑例 2-11中的联立方程，试找出 .2巾.x, y .2巾范围内所有准解析解。　　解可以用下面的命令直接求解联立方程： >> syms x y; t=cputime; F=[x^2*exp(-x*y^2/2)+exp(-x/2)*sin(x*y); y^2*cos(x+y^2)+x^2*exp(x+y)]; more_vpasols(F,zeros(0,2),4*pi); cputime-t %求根并计时要检验得出方程根的精度，则首先应该提取出感兴趣区域内的根 x0和 y0，并对其进行排序，得出的根代入原方程后的误差范数为 7.79×10.32，比例 2-37中得出的精度要高得多，所需的时间在半个小时左右，也远高于 more_sols()函数，用这样的方法只找到区域内的 105个根。得出的图形类似于图 2-12，不过有几个点缺失。 >> x0=X(:,1); y0=X(:,2); ii=find(abs(x0)<2*pi & abs(y0)<2*pi); x0=x0(ii); y0=y0(ii); [x0 ii]=sort(x0); y0=y0(ii); double(norm(subs(F,{x,y},{x0,y0}))), size(x0) fimplicit(F,[-2*pi,2*pi]), hold on; plot(x0,y0,'o') 　　例 2-44试求解例 2-40的高精度数值解。　　解可以试用下面语句求解方程，不过在用符号表达式描述方程时，并不能计算出方程的表达式 f，所以不能调用后续的 more_vpasols()函数，不能得出方程的高精度 ·40·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 数值解，求解这类方程只能采用例 2-40的普通数值解方法。 >> A=[2 1 9; 9 7 9; 6 5 3]; B=[0 3 6; 8 2 0; 8 2 8]; C=[7 0 3; 5 6 4; 1 4 4]; D=[3 9 5; 1 2 9; 3 3 0]; X=sym('x%d%d',3); f=expm(A*X)*funm(B*X,@sin)-C*X+D; more_vpasols(f,zeros(3,3,0),20); 2.6欠定方程的求解 前面介绍方程时，一直在假设方程的个数与未知数的个数是一致的，这些方程都是正常的方程，本节将探讨异常的方程类型欠定方程的概念与求解方法。 　　定义 2-11如果方程的个数等于未知数的个数，则方程称为适定方程（ well-posed equation，又称恰定方程）；如果方程的个数少于未知数的个数，则方程称为欠定方程（ underdetermined equation）；如果方程的个数大于未知数的个数，则方程称为超定方程（overdetermined equation）。 前面演示的隐式方程 f(x, y)=0就是一个常见的欠定方程，如果由 ezplot()或 fimplicit()函数用图解法求解，则得出的曲线上所有的点都满足原欠定方程，这时，欠定方程有无穷多解。 在一些特殊的场合下，用隐函数绘制函数不能绘制出任何曲线，这时方程可能有个别孤立解。这种情况下也可以考虑采用 fsolve()函数直接求解，不过在默认的设置下， fsolve()函数并不能求解方程与未知数个数不同的代数方程，需要将求解算法设置成 levenberg-marquardt，即采用 Levenberg–Marquardt算法求解欠定方程。如果采用 more_sols()函数，也应该作相应的算法设置。本节将通过例子演示具有孤立解的欠定方程求解方法。 　　例 2-45试求解下面的欠定方程。 32 32 (4x1 +4x1x2 +2x2 . 42x1 . 14)2 + (4x2 +2x1 +4x1x2 . 26x2 . 22)2 =0 　　解如果手工求解可以发现，原欠定方程可能分拆成两个独立方程，这样，方程的个数与未知数的个数一致，就可以调用 more_sols()类函数直接求解方程。手工转换的方法带有很多的人为性，因为并不是所有欠定方程都是可以手工拆分的。这里不作这种手工转换，试图直接求解欠定方程。如果尝试用下面的语句绘制隐函数曲线，在调用过程中没有任何警告信息，但最终不能得到任何曲线，说明方程只可能存在有限个孤立解。 >> f=@(x1,x2)(4*x1.^3+4*x1.*x2+2*x2.^2-42*x1-14).^2+... (4*x2.^3+2*x1.^2+4*x1.*x2-26*x2-22).^2; fimplicit(f) 得出该方程精确的解？ Levenberg–Marquardt现在可以人为地选择算法，再求解欠定方程，经过一段时间9的运行，有可能找出该欠定方程所有的个根。 ff=optimset;ff.Algorithm='levenberg-marquardt';>> F=@(x)f(x(1),x(2));%用匿名函数描述欠定方程 more_sols(F,zeros(2,1,0),1000,eps,300,ff)9得出的方程的 个解为 、、、、.(2805131313)(32)(0086728843)(3385200739).,.,.,..,.、、、和.......(3584418481)(3779332832)(0128019537)(3073000813).,..,..,..,.其中，搜索第四个解比较耗时。..(0270809230)，.,.本章习题 2.12-27试验证例中手工化简的是正确的。x12.2432试求解多项式方程。能否用数值方法+14+735+1715+1500625=0xxxx...2.32-39试验证例得出的一些无穷远处的解都不满足原方程。 2.42-35由于涉及复数矩阵，使用准解析解方法求解例 是很耗时的，试用数值解方法重新求解该方程，并体验复数矩阵是不是为方程求解带来额外的麻烦。 . .. 12222（） ...24=13xzxzxz.. 2.5求解能转换成多项式方程的联立方程，并检验得出的高精度数值解的精度。 2 24xy . x2 . y2 . xy2 = 13 2 24yz . y2 . z2 . yz2 = 13 . .. .. 2222 xy2 . zxy . 4xyz= xz （2） xy3 . 2yz2 =3x3z2 +4xzy2 224 yx . 7xy2 +3xz= xzy . .. .. x +3y3 +2z2 = 1/2 （3） x2 +3y + z3 =2 x3 +2z +2y2 = 2/4 2.6试求解下面的联立方程 [4]。 x1x2 + x1 . 3x5 =0 22 2x1x2 + x1 +3R10x2 + x2x3 + R7x2x3 + R9x2x4 + R8x2 . Rx5 =0 22 2x2x3 + R7x2x3 +2R5x3 + R6x3 . 8x5 =0 2 R9x2x4 +2x4 . 4Rx5 =0 22 22 x1x2 + x1 + R10x2 + x2x3 + R7x2x3 + R9x2x4 + R8x2 + R5x3 + R6x3 + x4 =1 其中，0.0001 . xi . 100，i =1, 2, 3, 4, 5，且已知常数 R = 10，R5 =0.193， R6 =4.10622×10.4，R7 =5.45177×10.4，R8 =4.4975×10.7，R9 =3.40735×10.5， R10 =9.615×10.7。 . . . . . ·42·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 2.7试求解下面的方程 [4]。 bb(1 + aT0) T T ec/T . ec/T + . 1=0 T0 aT0 T0 其中，100 . T . 1000，并已知常数 a = .1000/(3.H)，b =1.344 × 109， c = .7548.1193，T0 = 298，且 .H有三个取值，分别为 .50000，.35958和 .35510.3。 2.8试求解下面的联立超越方程 [4]。 { 0.5 sin x1x2 . 0.25x2/巾 . 0.5x1 =0 [] (1 . 0.25/巾)e2x1 . e+ ex2/巾 . 2ex1 =0 其中，0.25 . x1 . 1，1.5 . x2 . 2巾。 （1）文献 [4]给出了方程的两个解。若求解区间增大到 x1 → (.5, 5)，x2 →(.10, 10)，试求出方程全部的根，并用图解法验证得出的解，观察有没有没有找到的实根。 （2）该方程有实数根吗？在虚部约束为 (.10, 10)内总共可以找到多少根？ 2.9试求解下面方程中的 t并验证结果 [6]。 { t31 + t23y + t17x + t112 y+ t5xy + t2x2 =0 t37 + t29y + t19x + t132 y+ t7xy + t3x2 =0 2.10试求解带有参数的方程。{ 2 x+ ax2 +6b +3y2 =0 y = a + x +3 2.11试用图解法求解下面的一元和二元方程，并验证得出的结果。 （1）f(x)= e.(x+1)2+巾/2 sin(5x + 2) { 2 2.y (x+ y2 + 10xy)e.x 2.xy =0 （2） x3 +2y =4x +5 { 2.12用图解法和数值解法求出下面联立方程在 .2巾 . x, y . 2巾区域内全部的根 [7]。 2e.xy /2 x2+ e.x/2 sin(xy)=0 y2 cos(x + y2)+ x2ex+y =0 2. 13用数值求解函数求解习题 2.11中方程的根，并对得出的结果进行检验。 2. 14试求解下面的机器人动力学方程，看看总共可以找到多少实根 [4]。 . 4.731×10.3x1x3 .0.3578x2x3 .0.1238x1 +x7 .1.637×10.3x2 .0.9338x4 =0.3571 0.2238x1x3 +0.7623x2x3 +0.2638x1 . x7 . 0.07745x2 . 0.6734x4 . 0.6022 = 0 x6x8 +0.3578x1 +4.731×10.3x2 =0 . . .0.7623x1 +0.2238x2 +0.3461 = 0 22 x1 + x2 . 1=0 22 x3 + x4 . 1=0 2 x+ x6 . 1=0 52 . . 22 x7 + x8 . 1=0 且 .1 . xi . 1，i =1, 2, ··· , 8。试验证得出的方程的根。如果扩大求解区间，能否找到其他的根？该方程有复数根吗？ 2. 15试求出伪多项式方程 x 7 +2x 3 +3x 2.1 +4=0所有的根，并检验结果。 2.16试找出下面 Riccati变形方程全部的解矩阵，并验证得出的结果。 AX + XD . XBX + C = 0 其中， . . . . 2 1 9 0 3 6 7 0 3 3 9 5 A = . 9 7 9 . , B = . 8 2 0 . , C = . 5 6 4 . , D = . 1 2 9 . 6 5 3 8 2 8 1 4 4 3 3 0 2. 17已知上题给出的矩阵，试求解 AX + XD + CX2 . XBX + X2C + I = 0。 2. 18试求下面线性代数方程的解析解，并检验解的正确性。 2. 19已知下面的联立线性方程，试用 solve()函数得出并验证方程的解。 2 .9 3 .2 .1 . .1 .4 0 . . . 10 8 .1 .2 10 .4 5 .6 0 3 . . X = . . .3 0 .8 3 .4 3 . . .5 .6 .6 .8 .4 9 .5 3 . x1 + x2 + x3 + x4 + x5 =1. 3x1 +2x2 + x3 + x4 . 3x5 =2 . x2 +2x3 +2x4 +6x5 =3 5x1 +4x2 +3x3 +3x4 . x5 =4 . . 4x2 +3x3 . 5x4 = 12 2.20假设非线性方程为 AX3 + X4D . X2BX + CX . I = 0，且 A、B、C和 D矩阵在习题 2.16中给出，试求出该方程的全部实根。假设已经求出了方程的 77个实根，总共 3351个复数根，在 data2ex1.mat文件给出，试接着求解该方程，看看能不能找到新的解。 无约束最优化 第 3章 MATLAB最优化计算 最优化技术是当前科学研究中一类重要的手段。所谓最优化就是找出使得目标函数值达到最小或最大的自变量值的方法。毫不夸张地说，学会了最优化问题的思想与求解方法，可以将科研的水平提高一个档次，因为原来解决问题得到一个解就满足了，学会了最优化的思想后，很自然地将追求问题最好的解。最优化问题从其分类看分为无约束最优化问题和有约束最优化问题。 本章侧重于介绍无约束最优化问题以及 MATLAB求解方法，在 3.1节中先给出无约束最优化问题的定义与标准数学模型，然后介绍无约束最优化问题的解析解方法与图解法，并给出全局最优解与局部最优解的定义与判定方法，最后以简单的一元函数最优化问题为例，介绍最优化问题的算法与 MATLAB实现。 3.2节侧重于基于 MATLAB最优化工具箱函数的最优化问题求解方法，并通过例子演示相关求解函数的使用格式与应用技巧，演示梯度信息在最优化问题求解中的应用与效果，还将给出基于并行计算的最优化问题求解方法。 3.3节将探讨一般最优化问题的全局最优解方法，给出一个尝试求解问题全局最优解的思路和其 MATLAB实现，并通过一个改进的测试函数检验了该算法，验证了该算法的有效性。 3.4节介绍带有决策变量边界限制的最优化问题及其求解方法，并试图得到单变量与多变量最优化问题的全局最优解。 3.5节还将探讨如何使用最优化技术求解一些实际应用问题。 3.1无约束最优化问题简介 无约束最优化（unconstrained optimization）问题是最常见也是最简单的一类最优化问题。本节首先介绍无约束最优化问题的标准数学模型，并探讨最优化问题的解析解方法与图解方法，还将给出全局最优解与局部最优解的概念，以及简单一元最优化问题的求解算法并通过例子演示其 MATLAB实现。 3.1.1无约束最优化问题的数学模型　　定义 3-1无约束最优化问题的一般数学描述为 min f(x) （3-1-1） x 式中，x =[x1,x2, ··· ,xn]T称为决策变量；标量函数 f(·)称为目标函数。 上述数学描述的物理含义是如何找出求取一组 x向量，使得目标函数 f(x)的值为最小，故该问题又称为最小化问题。由于这里的 x向量的值可以任取，所以这类最优化问题又称为无约束最优化问题。 其实，这里给出的最小化是最优化问题的通用描述，它不失普遍性。若要想求解最大化问题，那么只需给目标函数 f(x)乘以 .1就能立即将其转换成最小化问题。所以本章及后续介绍中描述的全部问题都只考虑最小化问题，非最小化问题需要事先转换成最小化标准型。 3.1.2无约束最优化问题的解析解求解 f 无约束最优化问题的最优点 x处，目标函数 f(x)对 x各个分量的一阶导数为 0，从而可以列出下面的方程。 df =0, df =0, ··· , df =0（3-1-2） dx1 x=x扩 dx2 x=x扩 dxn x=x扩 求解这些方程构成的联立方程可以得出极值点。其实，解出的一阶导数均为 0的极值点不一定都是极小值的点，其中有的还可能是极大值点。极小值问题还应该有正的二阶导数。对于单变量的最优化问题，可以考虑采用解析解的方法进行求解。然而多变量最优化问题因为需要将其转换成求解多元非线性方程，其难度甚至高于直接最优化问题，所以没有必要用解析解方法求解。 3.1.3无约束最优化问题的图解法 如果一元方程或二元方程组已知，则可以将方程用匿名函数或符号表达式描述出来，然后采用 fplot()或 fimplicit()函数将方程曲线绘制出来，再通过读取曲线的交点信息即可以获得方程的解。 　　例 3-1考虑一元函数 f(t)= e.3t sin(4t + 2) + 4e.0.5t cos(2t) . 0.5，试用解析求解和图形求解的方法研究该函数的最优性。　　解可以先表示该函数，并解析地求解该函数的一阶导数，用 fplot()函数可以绘制出 t → [0, 4]区间内原函数与一阶导函数的曲线，如图 3-1所示。 >> syms t; y=exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)*cos(2*t)-0.5; %描述目标函数 ·46·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 y1=diff(y,t) %求取目标函数的一阶导函数 fplot([y,y1],[0,4]), line([0,4],[0,0]) %绘制一阶导函数曲线 图 3-1一元函数的导数和方程图解法 通过下面的计算可以得出导数函数为 f ′ (t)= .e.3t(3 sin(4t + 2) . 4 cos(4t + 2)) . 2e.t/2(cos 2t +4 sin 2t) 求解方程 f ′ (t)=0，则可以得出方程的两个解为 x1 =1.4528，x2 =3.0190。对原函数求二阶导数，则这两个点的导函数值分别为 z1 =7.8553，z2 = .3.646。 >> x1=vpasolve(y1,2), x2=vpasolve(y1,3) y2=diff(y,2); z1=subs(y2,t,x1), z2=subs(y2,t,x2) 其实，求解导函数等于 0的方程不比直接求解其最优值简单。用图解法可以看出，在这个区间内有两个点，即 A1和 A2，使得它们的一阶导函数为 0，但从其一阶导数走向看， A2点对应负的二阶导数值，所以该点对应于极大值点，而 A1点对应于正的二阶导数值，故为极小值点。 然而因为给定的函数是非线性函数，所以用解析法或类似的方法求解最小值问题一点都不比直接求解最优化问题简单。因此，除演示之外，不建议用这样的方法求解该问题，而直接采用最优化问题求解程序得出问题的解。 3.1.4局部最优解与全局最优解 一般的数值最优解方法都采用搜索方法，用户可以给出一个初始搜索点，从这个搜索点出发，采用不同的数值解算法，根据目标函数的实际情况找到下一步的一个点，再根据该点决策变量的目标函数值搜索下一个点。很显然，可以考虑迭代的方法一步一步地搜索出最优的点，使得目标函数的值最小。 对一元问题而言，由于目标函数可以表示为曲线形式，所以一般搜索方法可以形象地理解成，在初始搜索点处放置一个小球，让小球沿曲线滚下，这样最终小球将在某一个点处停下来，这时，小球的速度为零，即在这点处目标函数的导数为零，这样的点就是人们期望的最优点。下面将给出一个例子，介绍数值最优化中局部最优值与全局最优值的概念。 　　例 3-2假设目标函数为 y(t)= e.2t cos 10t + e.3t.6 sin 2t，0 . t . 2.5，试观察不同的初值能得出的最小值，并讨论局部最小值与全局最小值的概念。　　解可以由下面的语句直接绘制出目标函数在感兴趣区域的曲线，如图 3-2所示。 >> f=@(t)exp(-2*t).*cos(10*t)+exp(-3*(t+2)).*sin(2*t); %目标函数 fplot(f,[0,2.5]); 图 3-2函数定义域为 t . 0 如果将 A点设为初始搜索点，则小球滚下后经过反复滚动，最终可能停止在 t1点。如果初始点数值为 B点，则小球经过滚动后最终收敛到 t2点，该点目标函数的一阶导数也是零，所以该点也是原问题的最优解。如果初始值选择在 t较大的位置，可能找到的最优值还可能是 t3或 t4。 从这里的例子可以看出，一个已知的目标函数可能有多个最优值，但其目标函数的值可能是不同的。在这个例子中，从给定的区间看， t1点称为原问题的全局最优解，因为该点对应的目标函数值是最小的。其他各点， t2、t3和 t4又称为问题的局部最优解。 　　定义 3-2如果在一个邻域 R内存在 xf点，使得 f(xf) . f(x)对所有 x → R 成立，则 xf称为局部最优解。　　定义 3-3如果在实数域内存在一个 xf点，使得 f(xf) . f(x)成立，则决策变量 xf称为无约束最优化问题的全局最优解。 从实际问题求解的角度看，人们追求的是问题的全局最优解，而不是局部最优解，如何获得问题的全局最优解是研究者普遍感兴趣的问题。事实上，目前所有的最优化算法没有哪一种能保证能求出最优化问题的全局最优解，只能说“更可能”获得全局最优解。在实际应用中，局部最优解没有什么价值，因为其他很多 ·48·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 非最优解处的函数值都小于局部最优解处的目标函数值。例如，在图 3-2曲线上看， t → (0.174, 0.434)处的函数值都优于局部最优解 t2。 　　例 3-3重新考虑例 3-2中的问题，如果感兴趣的区域变成了 t → (.0.5, 2.5)，试重新分析问题的全局最优解与局部最优解。　　解现在再考虑更大些的定义域，即 t . .0.5，则用下面的语句能绘制出该函数在新定义域内的曲线，如图 3-3所示。由得出的曲线看，可能的最优解为 t0，t1，t2，t3和 t4。这时的 t1点已经不再是问题的全局最优值，而是一个局部最优值，新的全局最优值变成了 t0。从这个例子可以看出，同样的目标函数随着感兴趣区域的不同，可能有不同的全局最优解。 >> f=@(t)exp(-2*t).*cos(10*t)+exp(-3*(t+2)).*sin(2*t); fplot(f,[-0.5,2.5]); ylim([-2,1.2]) 图 3-3函数定义域为 t . .0.5 若将定义域扩展到 t → (.～, ～)区间，则原问题没有真正意义的全局最优解。 3.1.5数值求解算法的 MATLAB实现 求解最优化问题有许多种数值算法，本书给出一种简单的方法演示无约束最优化问题的数值求解，并给出 MATLAB的实现。 通常用到的数值最优化方法都是先选定初值 x0，然后用迭代方法得出原始问题的数值解。假设已知第 k步迭代的决策变量值为 xk，依据 Taylor级数展开技术，可以用一个二次函数 g(x)去逼近 x处的函数值 [8]。 g(x)= f(xk)+ f ′ (xk)(x . xk)+ 1 f ′′ (xk)(x . xk)2（3-1-3） 2 且已知 f(xk)= g(xk)，f ′ (xk)= g ′ (xk)，f ′′ (xk)= g ′′ (xk)。现在不去优化 f(x)函数，而考虑简单二次函数 g(x)的优化问题，该函数最优解存在的条件是 g ′ (x)=0，由其 推导出 0= g ′ (x)= f ′ (xk)+ f ′′ (xk)(x . xk)（3-1-4） 令 x = xk+1，则可以推导出 xk+1 = xk . ff ′′′ ((xxkk)) （3-1-5） 如果用后向差分数值导数计算函数的二阶导数，则可以得出迭代公式。 xk+1 = xk . xk . xk.1 f ′ (xk)（3-1-6） f′(xk) . f′(xk.1) 其实，仔细观察这里的求解方法，不难看出该方法事实上就是求解 f ′ (x)=0代数方程的方法，递推公式类似 Newton–Raphson迭代法，因为函数的二阶导数采用了后向差分算法取代，从而得出该方程的数值解。下面将通过例子演示如何用 MATLAB来求解一元最优化问题最优值的方法与过程。 　　例 3-4试利用这里给出的算法求解例 3-1中的问题。　　解如果想求解这样问题的最优解，则需要用匿名函数描述目标函数的导数，而例 3-1已经给出了其导函数的解析表达式，所以这里给出简单的匿名函数实现。再选择两个初始参考点 t0 =1，t1 =0.5，这样就可以由循环结构实现迭代过程。可以将循环条件设置为 |tk+1 . tk| <.，并假设 . = 10.5，意即两次搜索点之间的距离足够小，则停止循环。运行下面的语句，则可以得出函数的最优解。 >> df=@(t)-exp(-3*t).*(3*sin(4*t+2)-4*cos(4*t+2)) ... -2*exp(-t/2).*(cos(2*t)+4*sin(2*t)); t0=1; t1=0.5; t=[t0,t1]; while abs(t1-t0)>1e-5 t2=t1-(t1-t0)/(df(t1)-df(t0))*df(t1); t0=t1; t1=t2; t=[t, t2]; end得出的中间结果如下，最优解为 tf =1.4528。可以看出，对这里给出的简单问题而言，仅仅经过几步迭代，就可以得出原始问题的数值最优解。 t = [1, 0.5, 1.7388, 1.4503, 1.4534, 1.4528, 1.4528] 有了中间结果，则可以利用下面的语句描述目标函数和导函数曲线，标注出中间的搜索点并同时标注辅助线，如图 3-4所示。从这些点可以大致了解搜索的中间过程。 >> f=@(t)exp(-3*t).*sin(4*t+2)+4*exp(-0.5*t).*cos(2*t)-0.5; fplot(f,[0,4]), hold on fplot(df,[0,4]), plot(t,df(t),'o'), hold off 当然，这里仅给出了一元函数最优解的一种常用的搜索方法，对多元问题而言，搜索过程可能要麻烦得多。另外，由原始问题可能不能得出函数的一阶导数（对单元问题而言应该为梯度）信息，只知道目标函数本身的信息，求解全局最优化问题可能更加麻烦。本章后续内容将不再介绍底层的求解算法，只介绍基于 MATLAB最优化工具箱的直接求解方法。 ·50·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 图 3-4最优值搜索的中间结果 3.2无约束最优化问题的 MATLAB直接求解 MATLAB提供了两个高效的无约束最优化问题求解函数，本节将介绍如何利用这些函数直接求解无约束最优化问题的方法，并探讨如何提高求解的精度与效率，还将尝试利用并行计算的方法求解最优化问题。 3.2.1直接求解方法 MATLAB语言中提供了求解无约束最优化的函数 fminsearch()，其最优化工具箱中还提供了函数 fminunc()，二者的调用格式完全一致，为 x=fminsearch(Fun,x0), %最简求解语句 [x,f0,flag,out]=fminsearch(Fun,x0,opt), %带有控制选项 [x,f0,flag,out]=fminsearch(Fun,x0,opt,p1,p2,··· ), %带有附加参数式中， Fun为描述目标函数的函数句柄，它可以是 MATLAB函数的文件名，也可以是匿名函数或 inline函数； x0为搜索初值向量； opt为控制选项。除此之外，该函数还允许使用附加参数 p1，p2，···，但不建议使用附加参数，后面将介绍一种替代附加参数的方法。在返回的变量中， x为决策变量的最优解； f0为最优的传递函数值； flag为计算结果的标志，若 flag为 0，则说明求解过程不成功， flag为 1则表示求解成功；返回量 out包含一些求解的中间信息，如迭代步数等。 从数值算法看， fminsearch()函数采用了文献 [9]中提出的改进单纯形算法，而 fminunc()函数可以使用拟 Newton算法，也可以使用置信域（trust region）求解算法。从采用的算法看， fminsearch()函数无须目标函数的梯度信息，拟 Newton算法的 fminunc()函数也无须目标函数的梯度信息，而 fminunc()采用置信域算法时是需要用户提供梯度信息的。 从求解效果上作者发现，对很多算例而言，较新版本的 fminunc()函数效率得到了大幅度提升，所以建议尽可能使用新版本的求解函数。下面将通过例子来演示无约束最优化问题的数值解法。 　　例 3-5已知二元函数 z = f(x, y)=(x 2 . 2x)e.x .y .xy，试用 MATLAB提供的 22 求解函数求出其最小值，并解释解的几何意义。 　　解因为函数中给出的自变量是 x和 y，而最优化函数需要求取的是自变量向量 x，故在求解前应该先进行变量替换，如令 x1 = x，x2 = y，这样就可以将目标函数手工地修改成 z = f(x)=(x12 . 2x1)e.x.x1x2 2.x122 如果想求解最优化问题，首先需要将目标函数用 MATLAB表示出来。通常可以采用两种方法描述目标函数，其一是采用 MATLAB函数的方法，另一种是匿名函数的方法。先看一下 MATLAB函数的编程方法，在这种方法下需要由决策向量 x直接计算目标函数的值 y。 function y=c3mopt(x) y=(x(1)^2-2*x(1))*exp(-x(1)^2-x(2)^2-x(1)*x(2)); 将上述的函数存成 c3mopt.m文件后，则可以由下面的语句直接求解最优化问题，得出最优解为 x = [0.6111, .0.3056]，这时返回的标志 flag为 1,表示求解成功，另外可以看出，在求解过程中调用了 90次目标函数，总的迭代步数为 46。 >> [x,b,flag,c]=fminsearch(@c3mopt,[1;1]) %给出初值并求解最优化问题 另一种描述目标函数的方法是采用匿名函数的方法，这样做的一个好处是无须建立一个实体的 MATLAB函数文件，只需给出动态命令即可；另一个好处是可以直接使用 MATLAB工作空间中的变量。用下面的语句可以先用匿名函数定义出目标函数，然 后求解最优化问题，得出的解与前面的方法完全一致。 >> f=@(x)(x(1)^2-2*x(1))*exp(-x(1)^2-x(2)^2-x(1)*x(2)); %目标函数 x0=[2; 1]; [x,b,flag,d]=fminsearch(f,x0) %由初值并求解最优化问题同样的问题用 fminunc()函数求解，则可以得出同样的结果。这时，目标函数的调 用次数为 66，总的迭代步数为 7。 >> [x,b,flag,d]=fminunc(f,[2; 1]) %另一个求解函数比较两种方法，显然可以看出，用 fminunc()函数的效率明显高于 fminsearch()， 因为对目标函数调用的次数与迭代的步数都明显少于后者。所以在无约束最优化问题求解时，如果安装了最优化工具箱，则建议使用 fminunc()函数。其实，利用 surf()函数可以绘制出目标函数的曲面，如图 3-5所示。可见，得出的最小值点为得出曲面的谷底。 >> syms x y; f=(x^2-2*x)*exp(-x^2-y^2-x*y); fsurf(f,[-3 3, -2 2]) ·52·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 图 3-5给定函数的三维曲面图 　　例 3-6试求解 De Jong基准测试问题。 20 ∑ 2 min x i x i=1 　　解很显然，问题的解析解为 xi =0，i =1, 2, ··· , 20。现在可以用匿名函数描述目标函数，然后给出 fminsearch()直接求解最优化问题。其中，在描述目标函数时， x(:)命令将 x向量转换为列向量，所以这时不论给出的初值 x0为行向量还是列向量都能正确地计算目标函数的值。 >> f=@(x)x(:)'*x(:); x0=ones(20,1); [x0 f0 flag cc]=fminsearch(f,x0) 调用上述的求解语句，得出目标函数为 f0 =0.0603，距离期望的理论值有较大的差异，且给出警告信息“正在退出 :超过了函数计算的最大数目，请增大 MaxFunEvals选项。当前函数值： 0.060339”。从给出的提示看，搜索过程非正常结束，原因是目标函数计算的次数大于预设的 MaxFunEvals选项。同时，应该注意，这时返回的 flag为 0，说明求解不成功。 3.2.2最优化控制选项 MATLAB提供的最优化函数的搜索过程可以这样理解：由用户给定的初值 x0进行搜索，每搜索一步，需要计算 ||xk+1 . xk||与 |f(xk+1) . f(xk)|的值，一般情况下，当下面两组条件之一满足，则搜索过程终止。 ||xk+1 . xk|| <.1, |f(xk+1) . f(xk)| <.2（3-2-1） k1 >k1 max,k2 >k2 max （3-2-2） 式中， .1和 .2为用户指定的误差限，应该选择为很小的正数，这两个参数决定搜索结果的精度； k1是目标函数调用次数； k1 max为用户选定的目标函数计算的最大允许次数； k2是实际迭代步数； k2 max为最大允许的迭代步数， k1和 k2这两个选项的默认值均为 200n，其中，n为决策变量的个数。 如何修改这些控制选项呢？与前面介绍的解方程问题一样，需要调用 optimset读入控制选项的模板，然后修改有关的选项。模板 optimset支持的主要控制选项在表 3-1中给出，这些选项是以结构体成员变量的形式给出的，可以直接修改。 表 3-1最优化求解函数的控制选项表 选项名 选项说明 Display 中间结果显示方式，其值可以取 'off'表示不显示中间值， 'iter'表示逐步显示， 'notify'表示在求解不收敛时给出提示（默认选项），'final'只显示最终值 FunValCheck 检验目标函数的值是否有效，可选值为 'on'与'off'，如果设置为前者，则目标函 数的值为复数或 NaN时给出错误信息，默认选项为'off' GradObj 求解最优化问题时使用，表示目标函数的梯度是否已知，可以选择为'off'或'on' LargeScale 表示是否使用大规模问题算法，取值为'on'或'off'，一般几个变量的问题不必采 用该算法 MaxIter 方程求解和优化过程最大允许的迭代步数，若方程未求出解，可以适当增加该值， 即 k2 max MaxFunEvals 方程函数或目标函数的最大调用次数，即 k1 max OutputFcn 常用于中间结果的处理，用户可以编写一个函数，描述在每步迭代中如何处理中 间结果，后面将给出演示例子 TolFun 误差函数误差限控制量，即 .2，当目标函数满足式（ 3-2-1）即终止求解 TolX 解的误差限控制量，即 .1，当解满足式（ 3-2-1）即终止求解 在一般求解格式的调用语句中，正常情况下，当 .1与 .2条件同时满足时搜索过程正常结束，返回的标志 flag为正数，而 k1 max与 k2 max满足时为非正常结束，返回的标志 flag为 0。所以可以考虑依赖 flag的值搜索出有意义的最优值。 在较新版本的最优化工具箱中，还提供了另一套控制选项可以由 opt=optimoptions('fminunc')命令列出所有的选项，其中常用的选项在表 3-2中给出。该表中合并了一些收敛条件，使得控制选项设置更简洁，所以在较新版的最优化问题求解中可以考虑使用新的控制选项。值得指出的是，该控制选项不能用于 fminsearch()函数。后面将通过例子演示选项的设置方法。 　　例 3-7试用不同的算法求解例 3-5中的无约束最优化问题，并显示中间结果。　　解先考虑采用 fminsearch()函数直接求解，该函数采用了改进的单纯形法。另外，为比较算法，可以设置中间结果的显示，得出的目标函数值为 .0.641423723509499。 >> f=@(x)(x(1)^2-2*x(1))*exp(-x(1)^2-x(2)^2-x(1)*x(2)); ff=optimset; ff.Display='iter'; x0=[2; 1]; [x,b,c,d]=fminsearch(f,x0,ff) %单纯形法 得出的中间结果如下： ·54·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 表 3-2新版本最优化求解函数的控制选项表 选项名 选项说明 Algorithm 最优化搜索算法的选择与设置， 'quasi-newton'为拟 Newton算 法，而'trust-region'为置信域算法，前者为默认的算法。注意， 置信域算法需要已知目标函数的梯度，否则不能使用 Display 与表 3-1中的选项是一致的 SpecifyObjectiveGradient 求解最优化问题时使用，表示目标函数的梯度是否已知，可以选择 为 0或 1 MaxIterations 优化过程最大允许的迭代步数，默认值为 400 MaxFunctionEvaluations 目标函数的最大调用次数，默认值为 100倍决策变量个数 OutputFcn 常用于中间结果的处理，用户可以编写一个函数，描述在每步迭代 中如何处理中间结果，后面将给出演示例子 OptimalityTolerance 统一的最优化误差限控制量 Iteration Func-count min f(x) Procedure 01 0 1 3 0 initial simplex 2 5 -0.000641131 expand 3 7 -0.00181849 expand 4 9 -0.0132889 expand 5 11 -0.0654528 expand 6 13 -0.0835065 reflect ...略去了中间结果 45 88 -0.641424 contract inside 46 90 -0.641424 contract inside 现在再调用 fminunc()函数求解最优化问题，该函数采用的是拟 Newton算法，调用该函数也可以同样求解最优化问题，得出的目标函数值为 .0.641423726326004，略优于前面的寻优方法。 >> ff=optimoptions('fminunc','Display','iter'); [x,b,c,d]=fminunc(f,x0,ff) %拟 Newton法求解得出的中间结果如下，可以看出，迭代步数更少，求解算法更高效。 Iteration Func-count f(x) Step-size optimality 0 3 0 0.00182 1 24 -0.134531 872.458 0.324 2 36 -0.134533 0.001 0.324 3 48 -0.623732 172 0.205 4 54 -0.641232 0.311866 0.0357 5 60 -0.641416 0.329315 0.00433 6 63 -0.641424 1 0.000218 7 66 -0.641424 1 1.49e-08 　　例 3-8试求出例 3-6精确的解。 　　解前面介绍过，由于在默认调用格式下给出提示，目标函数的计算次数大于最大允许的次数，所以一种很自然的方法是增大 MaxFunEvals，不过这并非是一种好的方法。另一种方法是将得出的 x0作为初值重新搜索，再判定得出的 flag值是否为正，如 果为正则可以终止循环过程，否则继续搜索，直到得出所需的决策变量为止。为得到精确的结果，可以给出更严格的误差限。 >> f=@(x)x(:)'*x(:); x0=ones(20,1); flag=0; k=0; ff=optimset; ff.TolX=eps; ff.TolFun=eps; while flag==0, k=k+1 [x0 f0 flag cc]=fminsearch(f,x0,ff); end, norm(x0) 可以看出，经过 19次循环（k = 19），得到的 ||x0|| =4.0405×10.14，目标函数的值低至 f0 =1.6326×10.27，圆满地求解了原始问题。 　　例 3-9试利用 fminunc()函数重新求解例 3-6中的最优化问题。　　解假设初始搜索点为位于 (.1000, 1000)区间的均匀分布随机数，则可以给出下面的命令，由 fminunc()函数直接进行寻优计算。 >> f=@(x)x(:)'*x(:); x0=-1000+2000*rand(20,1); ff=optimset; ff.TolX=eps; ff.TolFun=eps; [x0 f0 flag cc]=fminunc(f,x0,ff), norm(x0) 得出 x0向量的范数为 f0 =3.0924×10.10，可以看出，得出的结果比较接近理论值。不过，如果想使用循环结构，则不会得出更好的求解结果。 从这里给出的例子可见， fminunc()函数可以从给定的初值迅速得出比较精确的结果，然后利用 fminsearch()函数求精确解，将 k1 max与 k2 max设置成大值即可。例如可以尝试下面的语句，这时目标函数的总调用次数为 27169，总耗时为 0.44s，误差的范 数为 ||x0|| =5.7339×10.13。 >> x1=-1000+2000*rand(20,1); x0=x1; k=0; tic [x0 f0 flag cc]=fminunc(f,x0,ff); ff.MaxFunEvals=100000; ff.MaxIter=100000; M=cc.funcCount; [x0 f0 flag cc]=fminsearch(f,x0,ff); toc, M=M+cc.funcCount, norm(x0) 如果单纯采用 fminsearch()函数，则目标函数的调用总次数为 73760，明显高于上述的方法，误差的范数为 ||x0|| =9.8226×10.16，总耗时为 0.56s，在求解时间和效率上低于上述的方法，但精度提高了很多。 >> x0=x1; tic, [x0 f0 flag cc]=fminsearch(f,x0,ff); toc, cc.funcCount, norm(x0) ·56·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 3.2.3附加参数的传递 在实际的最优化问题中，目标函数除了与 x有直接关系外，还可以带有其他附加参数，本节将演示附加参数的使用方法，还将探讨避免使用附加参数的方法。 　　例 3-10已知扩展的 Dixon问题。 . n/10 10∑i.1 ∑ 2 min .(1 . x10i.9)2 + (1 . x10i)2 +(xj . xj+1)2. x i=1 j=10i.9 如果 n = 50，试求解该最优化问题。　　解显然，原问题的解析解为 xi =1，i =1, 2, ··· ,n。如果 n是 10的倍数，则可以将原始的 x向量转换成 10行，n/10列的 X矩阵，这样，利用向量化方法就可以编写出下面给出的 MATLAB函数来描述目标函数。 function y=c3mdixon(x,n) m=n/10; X=reshape(x,10,m); y=0; for i=1:m y=y+(1-X(1,i))^2+(1-X(10,i))^2+sum((X(1:8,i).^2-X(2:9,i)).^2); end 注意，由于标准目标函数传递的格式下， n的值并不能直接传入目标函数，所以应该将 n用作附加参数，这样，该函数的输入变元除了 x之外，还有第二输入变元 n。 描述了目标函数，则可以考虑由下面的命令直接求解最优化问题，注意，在求解语句调用时，应该给出附加参数 n的值，与描述目标函数的 m函数必须保持一致，否则不能直接求解。经过 38步循环，可以得到较高精度的数值解，误差的范数为 1.2318×10.13。 >> ff=optimset; ff.TolX=eps; ff.TolFun=eps; n=50; x0=10*rand(n,1); flag=0; k=0; while flag==0, k=k+1 [x0,f0,flag,cc]=fminsearch(@c3mdixon,x0,ff,n);norm(x0-1) end, norm(x0-1)如果用 fminunc()替换求解程序 fminsearch()，则得出的解的精度并不能这么高，即使使用循环结构也不能得出精确的结果。 在较新版本的最优化工具箱中，很多描述似乎不适合使用附加参数，所以这里考虑一种替代的方法，即为 MATLAB文件编写一个匿名函数接口，直接传入附加参数的方法，因为匿名函数可以直接使用 MATLAB工作空间中的变量。这里将通过例子演示这样的方法。 　　例 3-11重新考虑例 3-10中的问题，试不采用附加参数重新求解最优化问题。　　解例 3-10中曾经给出了描述目标函数的 MATLAB函数 c3mdixon()，该函数使用了附加参数 n。在后面将介绍的结构体调用格式下并不允许使用附加参数，这时应该为其建立一个匿名函数接口，将原来的附加参数在 MATLAB工作空间内赋值，这样就可以调用下面的语句重新求解上述的无约束最优化问题。 >> n=50; f=@(x)c3mdixon(x,n); x0=rand(n,1); [x0,f0,flag,cc]=fminunc(f,x0) 3.2.4最优搜索的中间过程 表 3-1中提及了 OutputFcn选项，可以将该选项设置为响应函数的函数句柄，并编写出响应函数，以便在每一步迭代中显示或处理中间结果。响应函数的格式将通过例子演示。 　　例 3-12试找出并显示例 3-5最优化过程的各个中间点。　　解为截取寻优过程的中间点，可以用开关结构编写一个如下的输出处理函数。注意，响应函数的输入、返回变元都应该写成这里给出的固定格式，在寻优过程中 MATLAB的内在机制会自动生成这些输入变元的值。 function stop=c3myout(x,optimValues,state) stop=false; switch state %开关结构，监视中间结果 case 'init', hold on %初始化响应：设置坐标系保护 case 'iter', plot(x(1),x(2),'o'), %迭代响应：将中间结果用圆圈表示 text(x(1)+0.1,x(2),int2str(optimValues.iteration)); %迭代步数 case 'done', hold off %结束监控过程：取消坐标系保护 end 这样就可以在每步迭代中将中间结果标识出来了。要启动这样的监控过程，需要将 OutputFcn选项设置为 @c3myout。要演示整个优化过程，可以先绘制出原目标函数曲面的等高线图，选择初始搜索点 x0 = [2, 1]T，则可以用下面的语句开始带有监控的优化过程，这样就可以在等高线上叠印出中间搜索点，如图 3-6所示，中间搜索点做了编号与标记处理，如果两个中间点的距离特别小，发生重叠（为排版需要已手工移动标记， 图 3-6求解过程示意图 ·58·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 避免重叠），则说明这时的计算步长很小，接近于收敛值。 >> [x,y]=meshgrid(-3:0.1:3, -2:0.1:2); %生成网格矩阵 z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y); contour(x,y,z,30); %等高线图 f=@(x)(x(1)^2-2*x(1))*exp(-x(1)^2-x(2)^2-x(1)*x(2)); %目标函数 ff=optimset; ff.OutputFcn=@c3myout; x0=[2 1]; x=fminunc(f,x0,ff) %解最优化问题 3.2.5最优化问题的结构体描述 MATLAB最优化工具箱还支持用结构体变量来描述最优化问题，这样可以使最优化问题的描述更规范。可以建立一个结构体变量 problem，该结构体的各个成员变量在表 3-3给出。完成结构体的描述，就可以由下面的语句直接求解最优化问题的结构体变量 problem。 [x,fm,flag,out]=fminunc(problem) 注意， fminsearch()并不支持这样的调用格式，只能采用前面给出的一般调用方法使用，求解无约束最优化问题。 表 3-3无约束最优化结构体成员变量表 成员变量名 成员变量说明 objective 目标函数句柄 x0 初始搜索向量 options 控制选项的设置，可以由 problem.options=optimset语句设置成默认的选项，也 可像前面例子介绍的那样，自己修改控制选项，然后将选项赋给 options成员变量 solver 应该将其设置为'fminunc' 　　例 3-13试用结构体的方式重新描述并求解例 3-5中的无约束最优化问题。　　解可以用下面命令建立起最优化问题的结构体变量 P，然后调用 fminunc()函数即可以直接求解原始问题，得出的结果与例 3-5中的结果完全一致。 >> P.solver='fminunc'; P.options=optimset; %用结构体描述整个问题 P.objective=@(x)(x(1)^2-2*x(1))*exp(-x(1)^2-x(2)^2-x(1)*x(2)); P.x0=[2; 1]; [x,b,c,d]=fminunc(P) %直接求解最优化问题 　　例 3-14试用结构体格式描述例 3-11给出的最优化问题，并重新求解该问题。　　解由于结构体的格式并不支持附加参数的使用，所以应该仿照例 3-11介绍的方法，为目标函数 c3mdixon()设置一个匿名函数的接口，将 n的值传入匿名函数，这样就可以用结构体的格式描述原始的无约束最优化问题。可以考虑给出下面的语句重新求解最优化问题，得出的结果与前面的方法完全一致。 >> P.solver='fminunc'; n=50; P.options=optimset; P.objective=@(x)c3mdixon(x,n); %建立匿名函数句柄 P.x0=rand(n,1); [x,b,c,d]=fminunc(P) %直接求解最优化问题 3.2.6梯度信息与求解精度 有时最优化问题求解速度较慢，甚至无法搜索到较精确的最优点，尤其是变量较多的最优化问题，所以需要引入目标函数梯度，以加快计算速度，改进搜索精度。然而，有时计算梯度也是需要时间的，也会影响整个运算速度，所以实际求解时应该考虑是不是值得引入梯度的概念。 在利用 MATLAB最优化工具箱求解最优化问题时，也应该和目标函数在同一函数中描述梯度函数，亦即这时 MATLAB的目标函数应该返回两个变量，第一个变量仍然表示目标函数，第二个变量可以返回梯度函数。同时，还应该将求解控制变量的 GradObj属性设置成'on'，或将新的控制选项 SpecifyObjectiveGradient设置为 1，这样就可以利用梯度来求解最优化问题。 ()2 ()2 　　例 3-15试求解 Rosenbrock函数 f(x1,x2) = 100x2 . x12 +1 . x1的无约束最小值问题。　　解这个问题是英国著名学者 Howard Harry Rosenbrock（1920.2010）构造的用于测试最优化问题求解算法的测试问题。从目标函数可以看出，由于它为两个平方数的和，所以当 x2 = x1 =1时，整个目标函数有最小值 0。用下面语句可以绘制出目标函数的三维等高线图，如图 3-7所示。 >> [x,y]=meshgrid(0.5:0.01:1.5); %生成网格数据 z=100*(y.^2-x).^2+(1-x).^2; %计算目标函数 contour3(x,y,z,100), zlim([0,310]) %绘制三维等高线图 图 3-7 Rosenbrock目标函数的三维等高线图 由得出的曲线看，其最小值点在图中的一个很窄的白色带状区域内，故 Rosenbrock目标函数又称为香蕉函数，而在这个区域内的函数值变化较平缓，这就给最优化求值带 ·60·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 来很多麻烦。该函数经常用来测试最优化算法的优劣。现在观察用下面语句求解最优化问题，可以尝试用 fminunc()函数求解最优化问题，再将结果作为初值继续搜索，经过 100次循环，即可求解最优化问题。 >> f=@(x)100*(x(2)-x(1)^2)^2+(1-x(1))^2; %目标函数的匿名函数描述 ff=optimset; ff.TolX=eps; ff.TolFun=eps; x=[0;0]; for k=1:100, k [x,f0,flag,cc]=fminunc(f,x,ff) %最优化问题的数值求解 end 这时得出的最优解为 x = [0.999995588079578, 0.999991166415393]T。可见，即使经过了长时间的运算，该算法也无法精确搜索到真值 (1, 1)，用传统的最速下降法更无法搜索到真值，所以这时需要引入梯度的概念。对给定的 Rosenbrock函数，利用符号运算工具箱即可以求出其梯度向量。 >> syms x1 x2; f0=100*(x2-x1^2)^2+(1-x1)^2; J=jacobian(f0,[x1,x2]) %梯度计算 可以求出梯度向量为 [ ] 2 2 J = . 400(x2 . x1)x1 . 2+2x1, 200x2 . 200x 1 这时，可以在目标函数中描述其梯度，故需要重新编写目标函数为 function [y,Gy]=c3fun3(x) y=100*(x(2)-x(1)^2)^2+(1-x(1))^2; %需要返回两个变量，不能用匿名函数 Gy=[-400*(x(2)-x(1)^2)*x(1)-2+2*x(1); 200*x(2)-200*x(1)^2]; 这样，就应该给出如下命令得出 x = [1.000000000000012, 1.000000000000023]T。 >> ff.GradObj='on'; x=fminunc(@c3fun3,[0;0],ff) %利用梯度重新求解可见，引入了梯度则可以明显加快搜索的进度，且最优解也基本上逼近真值，这是不使用梯度不可能得到的，所以从本例可以看出梯度在搜索中的作用。然而，在有些例子中引入梯度也不是很必要，因为梯度本身的计算和编程需要更多的时间。如果用结构体的方式描述最优化问题，则可以得出与前面一致的解。 >> problem.solver='fminunc'; ff=optimset; problem.x0=[2; 1]; problem.objective=@c3fun3; ff.GradObj='on'; ff.TolX=eps; ff.TolFun=eps; problem.options=ff; [x,b,c,d]=fminunc(problem) %用结构体描述整个问题，直接求解 如果不使用梯度信息，但采用 fminsearch()函数进行求解，则可以得出原问题精确的解。可见，不基于梯度信息的求解方法也可能较好地求解这类问题。 >> ff.GradObj='off'; x=fminsearch(f,[0;0],ff) %直接求解 值得指出的是， Rosenbrock函数是为检测寻优算法优劣而建立起来的人造函数，解决该问题的有效方法需要引入目标函数的梯度。实际应用中，很多寻优算法都是无须梯度信息的，利用目标函数本身的信息即可成功地解决数值寻优的问题。 　　例 3-16重新考虑例 3-6中介绍的最优化问题，由于采用的最优化算法未利用梯度信息，所以得出的结果还是比较慢的，试利用梯度信息重新求解该问题。　　解由给定的目标函数，可以立即看出目标函数的梯度为 2[x1,x2, ··· ,x20]=2x，所以可以利用 MATLAB函数的形式将目标函数与梯度信息表示出来。注意，由于该函数需要返回两个变元，所以不能采用匿名函数的形式描述目标函数与梯度。 function [y,G]=c3mdej1(x) y=x(:).'*x(:); G=2*x; 有了相应的信息，则可以由下面的语句直接最优化问题。可以看出，利用梯度信息，则用 3步迭代即可以得出问题的最优解， x向量的范数低至 9.0317×10.31，目标函数的值低至 f0 =8.1572×10.61，精度远高于其他的数值算法。 >> clear problem; problem.solver='fminunc'; ff=optimset; problem.objective=@c3mdej1; ff.GradObj='on'; ff.TolX=eps; ff.TolFun=eps; problem.options=ff; problem.x0=-100+200*rand(20,1); [x,f0]=fminunc(problem), norm(x) %用结构体描述问题并求解 　　例 3-17试利用置信域算法重新求解例 3-5的最优化问题。　　解如果采用置信域算法，则需要首先获得目标函数的梯度。可以使用符号运算的方式，直接推导目标函数的梯度表达式。 >> syms x1 x2; f=exp(-x1^2-x1*x2-x2^2)*(x1^2-2*x1) G=simplify(jacobian(f,[x1,x2])) 这样可以直接推导出目标函数的梯度为 22 21 [ ] 2 23 2x1 +2x1x2 . x1x2 +4x1 . 2x1 . 2 .x1x2.x G(x)= e.x 2 (.x1 +2x1)(x1 +2x2) 有了目标函数与目标函数的梯度数学表达式，则可以建立下面的 MATLAB函数，直接描述这两个量，通过下面的函数返回。 function [f,Gy]=c3mgrad(x) u=exp(-x(1)^2-x(1)*x(2)-x(2)^2); f=u*(x(1)^2-2*x(1)); Gy=[2*x(1)+2*x(1)*x(2)-x(1)^2*x(2)+4*x(1)^2-2*x(1)^3-2; (-x(1)^2+2*x(1))*(x(1)+2*x(2))]*u; 由相关的信息构造出最优化问题的结构体变量，可以将求解算法直接设置为置信域算法，然后给出 fminunc()函数直接求解问题。 >> clear problem; problem.solver='fminunc'; ff=optimset; problem.objective=@c3mgrad; ff.Algorithm='trust-region'; ff.GradObj='on'; ff.Display='iter'; problem.options=ff; problem.x0=[2;1]; [x,f0]=fminunc(problem) ·62·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 下面是得出的中间结果： Norm of First-order Iteration f(x) step optimality CG-iterations 0 0 0.00182 1 0 10 0.00182 1 2 0 2.5 0.00182 0 3 -0.0120683 0.625 0.0558 0 4 -0.46831 1.25 0.434 1 5 -0.46831 1.25226 0.434 1 6 -0.602878 0.313066 0.252 0 7 -0.640866 0.296542 0.0372 1 8 -0.641424 0.0310125 6.49e-05 1 9 -0.641424 6.10566e-05 6.96e-10 1 从这里给出的例子来看，执行效率显然低于例 3-7测试的拟 Newton算法，所以可以得出结论：并不是所有的无约束最优化问题都适合使用梯度信息的。 3.2.7离散点最优化问题的求解 在实际应用中，某一个需要优化的目标函数有时其原型是未知的，只有一些相应的、离散分布的样本数据点，这样，就可以采用样条插值或其他插值方法去拟合目标函数，从而优化这样的目标函数。下面分别介绍一元、二元函数甚至多元函数的插值方法，由这些插值方法即可以计算出感兴趣点的目标函数值。 （1）一元问题。如果已知样本点 x0和 y0，则对应任意的 x向量，可以由 f=@(x)interp1(x0,y0,x,'spline')匿名函数命令计算出 x向量各点的函数插值结果。 （2）二元问题。如果已知样本点向量 x0、y0和 z0，则可以令 p1 = x，p2 = y，这样，函数值的插值结果可以由函数句柄表示。 f=@(p)griddata(x0,y0,z0,p(1),p(2),'v4') （3）三元问题。如果已知样本点向量 x0、y0、z0与 v0，则可以令 p1 = x，p2 = y， p3 = z，这样，函数值的插值结果可以由函数句柄表示。 f=@(p)griddata(x0,y0,z0,v0,p(1),p(2),p(3),'v4') （4）多元函数的散点插值还可以考虑使用 griddatan()函数直接构造。 f=@(p)griddatan([x1(:),x2(:),··· ,xm(:)],z,p)有了函数值，则可以调用 fminunc()等函数直接求解最优化问题。下面将通过 例子演示这样的最优化问题求解方法。 　　例 3-18重新考虑例 3-5中的函数，试由已知的函数在 x → [.3, 3]，y → [.2, 2]区间生成均匀分布的 200个样本点，然后仅利用这些离散的散点求出对应函数的最小值，并检验所得出的结果。　　解仿照前面的例子，可以首先生成一些离散点，再由这些离散点通过插值方法构造目标函数的匿名函数，对该匿名函数进行优化，则可以得出最优解为 x =0.6069， y = .0.3085。可以看出，这样得到的最优点很接近例 3-5由目标函数表达式得出的最优解，所以这里给出的优化方法是可行的。 >> x=-3+6*rand(200,1); y=-2+4*rand(200,1); z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y); %生成散点数据 f=@(p)griddata(x,y,z,p(1),p(2),'v4'); %重建目标函数 x=fminunc(f,[0,0]) %由散点数据进行寻优 3.2.8最优化问题的并行求解 通常，大规模最优化问题可能涉及巨大的运算量，所以在实际应用中可以考虑最优化工具箱提供的自动并行计算功能，这需要如下设置控制选项 options=optimoptions('fminunc','UseParallel',true) 如果并行计算的选项开启，则求解机制会自动使用并行运算机制，利用并行计算的方法求解最优化问题。如果关闭并行计算机制，则采用普通的方法求解最优化问题。下面将给出较大规模问题的求解实例，不过从无约束最优化问题求解看，开启并行计算的意义并不是很大。 　　例 3-19考虑例 3-11中给出的最优化问题，如果 n = 2000，试比较使用并行计算与不使用并行计算的求解方法。　　解如果使用并行计算的方法，则需要给出下面的语句，这时耗时为 107.49s，得出的误差为 5.3493×10.6。MATLAB自动将计算任务分派给 4个处理器并行处理。 >> problem.solver='fminunc'; n=2000; ff=optimoptions('fminunc'); ff.OptimalityTolerance=eps; ff.MaxFunctionEvaluations=1000000; ff.MaxIterations=1000000; ff.UseParallel=true; problem.options=ff; %使用并行计算 problem.objective=@(x)c3mdixon(x,n); %建立匿名函数句柄 x0=rand(n,1); problem.x0=x0; tic, [x,b,c,d]=fminunc(problem); toc, norm(x-1) 如果不采用并行计算，则可以给出下面语句直接求解。为对比方便采用了相同的初值，这时得出的误差与前面得出的完全一致，耗时为 147.22s，略高于并行计算的耗时，但对此例而言并行计算的优势不明显。 >> ff.UseParallel=false; problem.options=ff; %关闭并行计算 tic, [x,b,c,d]=fminunc(problem); toc, norm(x-1) %重新计算 ·64·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 3.3全局最优解的尝试 前面介绍了全局最优解问题。本节将给出一个多谷底的基准测试函数，首先演示该测试函数最优解对初始值的依赖性，然后编写一个能尝试求出全局最优解的 MATLAB函数，并探讨其求解效率。 　　例 3-20考虑一个著名的无约束最优化问题的基准测试函数（ benchmark prob-lem） Rastrigin函数 [10]。 f(x1,x2) = 20+ x12 + x22 . 10(cos 2巾x1 + cos 2巾x2) 试绘制出目标函数的表面图，并用简单的最优化求解函数求解这样的问题，看看会发生什么。　　解目标函数的表面图可以由下面的语句直接得出，如图 3-8所示。可以看出，表面图凹凸不平，其中有很多波峰与波谷，全局最优解的图解法留作习题 3.20。 >> f=@(x1,x2)20+x1.^2+x2.^2-10*(cos(2*pi*x1)+cos(2*pi*x2)); fsurf(f) % 绘制目标函数的表面图 图 3-8 Rastrigin函数的表面图 还可以得出目标函数的等高线图，如图 3-9所示。从给出的图形可见，中间的点是全局最小点，另外还有很多波谷点，但它们都是局部最小点。另外，全局最优点附近的几个点可以认为是次最优（subminimum）点。 >> fcontour(f,'MeshDensity',1000) %绘制等高线图 选择几个不同的初始搜索点，则可以由下面语句得出不同的优化结果。在重新定义目标函数时，没有必要全盘改写原有的匿名函数，只需在原来函数的基础上写一个新的接口函数。有了目标函数的匿名函数描述，则可以在以下几个初值下直接寻优计算。 >> f1=@(x)f(x(1),x(2)); %可以在前面匿名函数的基础上定义新函数 x1=fminunc(f1,[2,3]), f1(x1), x2=fminunc(f1,[-1,2]), f1(x2) x3=fminunc(f1,[8 2]), f1(x3), x4=fminunc(f1,[-4,6]), f1(x4) 得到如下几个结果： 图 3-9目标函数的等高线图 x1 =[1.9602, 1.9602],f(x1)=7.8409, x2 =[.0.0000, 1.9602],f(x2)=3.9205 x3 =[7.8338, 1.9602],f(x3)=66.6213, x4 =[.3.9197, 5.8779],f(x4)=50.9570 可以看出，这样得出的最优化结果都是“最优”的，但有显著差异，大多数点为局部最小值点。可以看出，采用传统的最优化搜索方法，如果初始值选择不当，很可能陷入局部最小值。另外，上面并未得到全局最优解 x = [0, 0]。 为避免局部最小值问题，经常采用某些智能优化方法，如遗传算法（ genetic algorithm）或其他进化类算法，这类方法将在后续章节中给出概略性的介绍。不过即使遗传算法这类的方法也不能确保得到全局最优解，只不过进化类算法号称更可能得出全局最优解。 类似于前面介绍的方程求解的思路，可以采用下面的新算法来作全局寻优。首先，用随机的方式在感兴趣区域 (a, b)选择初值，则通过普通的搜索方法得出最优解 x，并得出最优目标函数 f1 = f(x)，如果得出的最优目标值比已经得到的还小，则记录该最优值。重复 N次这类求解过程，则可能得出问题的全局最优解。基于此思路，可以编写出如下的 MATLAB函数求解全局最优化问题。 function [x,f0]=fminunc_global(f,a,b,n,N,varargin) k0=0; f0=Inf; if strcmp(class(f),'struct'), k0=1; end %可以用结构体描述 for i=1:N, x0=a+(b-a).*rand(n,1); %用循环结构生成随机初始搜索点 if k0==1, f.x0=x0; [x1 f1 key]=fminunc(f); %结构体描述问题求解 else, [x1 f1 key]=fminunc(f,x0,varargin{:}); end %无约束最优化 if key>0 & f1> f=@(x)20+x(1)^2+x(2)^2-10*(cos(2*pi*x(1))+cos(2*pi*x(2))); [x,f0]=fminunc_global(f,-2*pi,2*pi,2,50) %尝试获得全局最优解为进一步演示这样的全局最优解求解过程，可以用循环调用 100次这一求解程序，可以看到，每次都能找到全局最优解。 >> F=[]; N=50; %建立一个目标函数值的存储向量，初始值为空白向量 for i=1:100 %调用求解函数 100次并评估找到全局最优解的成功率 [x,f0]=fminunc_global(f,-2*pi,2*pi,2,N); F=[F,f0]; end 当然，由于使用了均匀分布的随机数，这样的全局最优解 x1 = x2 =0很容易被找到。所以用这个例子评估全局优化算法并不公平，后面将试图给出更公平的测试函数评价全局最优算法。 　　例 3-22假设将经典的 Rastrigin函数修改成 ()2 ()2 [() ()] f(x1,x2) = 20+ x1 . 1+ x2 . 1. 10 cos 2巾 x1 .1+ cos 2巾 x2 .1 3020 30 20 可以运行 100次这个求解程序，测试一下找到全局最优解的成功率是多少。　　解显然，经过这样的改写，可以得出最优化问题全局最优解的理论值为 x1 = 30， x2 = 20。如果将感兴趣搜索区间扩展到 ±100，则可以进行如下的测试。 >> f=@(x)20+(x(1)/30-1)^2+(x(2)/20-1)^2-... 10*(cos(2*pi*(x(1)/30-1))+cos(2*pi*(x(2)/20-1))); F=[]; tic, N=100; for i=1:100 %运行该求解函数 100次 [x,f0]=fminunc_global(f,-100,100,2,N); F=[F,f0]; %记录最优值 end, toc 可以看出，这次执行 100次寻优，有 17次没有找到全局最优值 (30, 20)，其余的 83次都找到了全局最优点，成功率 83％，总耗时 57.07s。还可以看出，这里给出的函数 fminunc_global()是可信赖的，一般情况下很可能得出全局最优解。如果将 N改为50，则总耗时降至 26.23s，全局最优解成功率为 67％。 从这里给出的搜索算法看，应该很适合使用并行计算的策略加快搜索过程，不过若尝试使用并行计算，则 MATLAB执行机构会给出错误信息，指出匿名函数句柄的处理不能进行并行资源的分派，从而不能使用并行计算的方法求解问题。 3.4带有决策变量边界的最优化问题 本章前面的内容介绍了理论上的无约束最优化问题，但在很多演示例子中，采用的是指定决策变量范围的“无约束”最优化问题。例如，例 3-2给出的函数在不同感兴趣区间内的全局最优解是不同的，都不是数学意义下的无约束最优化问题。所以，本节将考虑带有决策变量边界的最优化问题的求解方法，首先介绍单变量最优化求解问题，然后给出多变量最优化的求解方法。 3.4.1单变量最优化问题 本节先给出单变量带决策变量边界的无约束最优化问题的定义，然后介绍 MATLAB提供的直接求解函数，再介绍全局优化的解决方法。　　定义 3-4单变量带决策变量边界的最优化问题的数学形式为 min f(x)（3-4-1） x s.t.xm.x.xM 式中， xm与 xM为决策变量的下界与上界；记号 s.t.是英文 subject to的缩写，表示满足后面的关系。 MATLAB最优化工具箱提供了 fminbnd()函数，可以用于直接求解带有决策变量的单变量最优化问题，该函数的调用格式为 x=fminbnd(f,xm,xM) [x,f0,flag,out]=fminbnd(f,xm,xM,options) [x,f0,flag,out]=fminbnd(problem)式中，f是单变量目标函数的句柄；xm和 xM为决策变量的下界与上界。如果采用结构体 problem描述最优化问题，则 xm和 xM相应的成员变量名为 x1和 x2。 　　例 3-23试重新求解例 3-3中的最优化问题。　　解显然，这个最优化问题带有决策变量的边界 tm = .0.5，tM =2.5，这样，可以由下面的命令直接求解最优化问题，得出的结果为 tf =1.5511。遗憾的是，这样的求解函数并不能得出感兴趣区域的全局最优解，即图 3-3中标注的 t0，而是局部最优解 t3，而 fminbnd()函数本身并不能选择初值，所以，只利用该函数并不能有效地求解全局最优化问题。 ·68·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 >> f=@(t)exp(-2*t).*cos(10*t)+exp(-3*(t+2)).*sin(2*t); tm=-0.5; tM=2.5; t=fminbnd(f,tm,tM) 由于该函数本身没有可调的参数，所以为得到全局最优解，需要将 (xm,xM)区间划分成 m个子区间，对每个子区间都求取最小值，找出其中目标函数值最小的一个，则这个值很可能是问题的全局最优解。利用这样的思路可以编写出下面的 MATLAB函数，该程序中使用的 fminsearchbnd()函数后面将介绍。 function [x,f0]=fminbnd_global(f,xm,xM,n,m,varargin) f0=Inf; M=ones(n,1); if length(xm)==1, xm=xm*M; end, if length(xM)==1, xM=xM*M; end for i=1:m if n==1, h=(xM-xm)/m; %单变量问题 [x1,f1]=fminbnd(f,xm+(i-1)*h,xm+i*h,varargin{:}); else, x0=xm+(xM-xm).*rand(n,1); %多变量问题 [x1 f1 key]=fminsearchbnd(f,x0,xm,xM,varargin{:}); end %用循环结构生成随机初始搜索点 if f1> f=@(t)exp(-2*t).*cos(10*t)+exp(-3*(t+2)).*sin(2*t); tm=-0.5; tM=2.5; x=fminbnd_global(f,tm,xM,1,10) 3.4.2多变量最优化问题 如果定义 3-4中的 x为 x向量，则相应的最优化问题可以改写成 min f(x) （3-4-2） x s.t. xm.x.xM 所以式（ 3-4-2）所描述的问题是， x在指定的范围内取多少时能使得目标函数取最优值。这样的问题由 fminsearch()和 fminbnd()函数是不能直接求解的。 John D’Errico提出了一种变换方法 [11]，引入新决策变量 zi，使得 xi = xmi + 1(xMi . xmi )(sin zi + 1)（3-4-3） 2 这样，就可以将关于 xi的决策变量边界问题巧妙地转换成关于 zi的无约束最优化问题。 John D’Errico还开发了 fminsearchbnd()函数，扩展了现有 fminsearch() 函数的功能，使其能直接求解式（ 3-4-2）中的问题，该函数的调用格式为 [x,f0,flag,out]=fminsearchbnd(f,x0,xm,xM) [x,f0,flag,out]=fminsearchbnd(f,x0,xm,xM,opt) [x,f0,flag,out]=fminsearchbnd(f,x0,xm,xM,opt,p1,p2,··· )如果上界或下界约束没有给出，则可以将其设置为空矩阵 []，这时，该函数会 自动调整式（ 3-4-3）中的变换公式，仍然能正常求解最优化问题。 　　例 3-25试求解下面函数的最小值。 ()2 ()2 2 3 f(x, y) = (1.5 . x + xy)2 +2.25 . x + xy +2.625 . x + xy 其中，.4.5 . x, y . 4.5。　　解令 x1 = x，x2 = y，则可以将目标函数改写成 ()2 ()2f(x) = (1.5 . x1 + x1x2)2 +2.25 . x1 + x1x 2 +2.625 . x1 + x1x 3 2 2 可以由匿名函数的形式直接由点运算描述目标函数，然后绘制出目标函数的曲面，如图 3-10所示。 >> f=@(x,y)(1.5-x+x.*y).^2+(2.25-x+x.*y.^2).^2+(2.625-x+x.*y.^3).^2; fsurf(f,[-4.5,4.5]) 图 3-10目标函数的曲面 如果想求解原始的最优化问题，并没有必要真采用上面 f(x)的格式重新输入目标函数，只需按照下面的形式改写目标函数，就可以构造出目标函数句柄 f0，这时再求解最优化问题，则可以得出全局最优解为 x1 = [3, 0.5]T，这时的目标函数为 f(x1)=2.6871×10.29。 >> f0=@(x)f(x(1),x(2)); opt=optimset; opt.TolX=eps; ff.TolFun=eps; [x1,f1]=fminsearchbnd(f0,[0;0],-4.5*[1;1],4.5*[1;1],opt) 如果将感兴趣区间设置为 .2 . x, y . 2，则显然刚才得到的最优解不在此区域内，所以可以利用下面的语句重新搜索，得出的最优解为 x1 = [2, 0.1701]T，这时的最优目 ·70·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 标函数值为 f1 =0.5233。 >> [x1,f1]=fminsearchbnd(f0,[0;0],-2*[1;1],2*[1;1],opt) fsurf(f,[-2,2]) 3.4.3边界问题全局最优解的尝试 对多峰的目标函数或者大范围搜索，前面介绍的 fminsearchbnd()函数可能不能确保得出问题的全局最优解，所以需要引入全局最优的方法。前面编写的 fminbnd_global()函数可以尝试求取最优化问题的全局最优解。 　　例 3-26如果求解区间变成 .500 .x, y .500，试求解例 3-25中的最优化问题。　　解如果使用 fminsearchbnd()函数，当然可以给出下面的语句，不过得出的结果可能为 x1 =[.50, 1.0195]T，同时给出警告信息， 且 f0 =0.4823，“超过了函数计算的最大数目，请增大 MaxFunEvals选项”，说明求解不成功。反复运行下面最后一条指令可以发现，大约有一半的运行得不出全局最优解 [3, 0.5]。 >> f=@(x,y)(1.5-x+x.*y).^2+(2.25-x+x.*y.^2).^2+(2.625-x+x.*y.^3).^2; f0=@(x)f(x(1),x(2)); opt=optimset; opt.TolX=eps; ff.TolFun=eps; [x1,f1]=fminsearchbnd(f0,500*rand(2,1),-500*[1;1],500*[1;1],opt) 调用前面编写的 fminbnd_global()函数重新求解原始问题，则反复运行下面的命令，可以看出 100％的运行都能得到原始问题的全局最优解。 >> [x1,f1]=fminbnd_global(f0,-500,500,2,5,opt) 3.5最优化问题应用举例 本节将给出几个最优化技术应用的实例。首先引入线性回归问题的求解方法，如果已知实验数据，通过实验数据将其对应的数学模型重建出来。这类方法对应于超定线性方程的最小二乘求解方法。如果函数的数学模型的形式过于复杂，本节还将介绍基于最小二乘方法的曲线拟合技术。另外，本节还将通过例子演示基于最优化技术的微分方程边值问题与代数方程的求解方法。 其实，本节给出的几个例子是想演示一下，如何将看起来不相关的问题转换成最优化问题，通过最优化问题的建模与求解，得出原始问题的解。 3.5.1线性回归问题的求解 假设已知某函数的线性组合为 g(x)= c1f1(x)+ c2f2(x)+ c3f3(x)+ ··· + cnfn(x)（3-5-1） 式中， f1(x),f2(x), ··· ,fn(x)为已知函数； c1,c2, ··· ,cn为待定系数，这时假设已经 测出数据 (x1,y1), (x2,y2), ··· , (xm,ym)，则可以建立如下线性方程： Ac = y （3-5-2） ... 式中 f1(x1) f2(x1) ··· fn(x1) y1 ... A = f1(.. x2) f2(.. x2) · .· .. · fn(.. x2) , y = y.. 2 （3-5-3） ... . f1 (xm) f2 (xm) ··· fn (xm) ym ... 且 c =[c1,c2, ··· ,cn]T。故该方程的最小二乘解为 c = A\y。 　　例 3-27假设测出了一组 (xi,yi)由表 3-4给出，且已知函数原型如下，试用已知的 数据求出待定系数 ci的值。 2 y(x)= c1 + c2e.3x + c3 cos(.2x)e.4x + c4x 表 3-4例 3-27实测数据 ... . . xi 0 0.2 0.4 0.7 0.9 0.92 0.99 1.2 1.4 1.48 1.5 yi 2.88 2.2576 1.9683 1.9258 2.0862 2.109 2.1979 2.5409 2.9627 3.155 3.2052 　　解可以由表 3-4中给出的数据直接拟合出曲线方程中的 ci参数。这样，就可以依照各个子函数的表达式，由点运算的方式构造 A矩阵，再使用最小二乘命令求出各个待定系数。 >> x=[0,0.2,0.4,0.7,0.9,0.92,0.99,1.2,1.4,1.48,1.5]'; %样本点输入 y=[2.88;2.2576;1.9683;1.9258;2.0862;2.109;2.1979; 2.5409;2.9627;3.155;3.2052]; A=[ones(size(x)) exp(-3*x),cos(-2*x).*exp(-4*x) x.^2]; c=A\y; c1=c' %最小二乘 可以得出拟合参数 cT = [1.22, 2.3397, .0.6797, 0.87]，将更密集的 x向量代入该原型函数，拟合曲线和已知数据点如图 3-11所示，可见拟合效果是令人满意的。 >> x0=[0:0.01:1.5]'; B=[ones(size(x0)) exp(-3*x0) cos(-2*x0).*exp(-4*x0) x0.^2]; y1=B*c; plot(x0,y1,x,y,'x') %拟合曲线计算与绘制 3.5.2曲线的最小二乘拟合 前面介绍的线性回归方法中，首先假定 fi(x)为已知函数，所以待定系数的计算才能转换成线性代数方程的求解问题。在实际应用中， fi(x)这类函数本身还可以带有其他的待定系数，所以这里问题不能由线性代数方程求解出来，必须将其转换为最优化问题进行求解。 ·72·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 图 3-11原始数据与拟合曲线 　　定义 3-5假设有一组数据 xi，yi，i =1, 2, ··· ,m，且已知这组数据满足某一函数原型 y.(x)= f(a,x)，其中， a为待定系数向量，则最小二乘曲线拟合的目标就是求出这一组待定系数的值，可以将问题转换成下面的最优化问题： mm ∑[]2 ∑[]2 J = min yi . y.(xi)= min yi . f(a,xi)（3-5-4） aa i=1 i=1 该最优化问题可以由底层求解的方式求解，也可以调用 MATLAB的最优化工具箱中提供的 lsqcurvefit()函数，该函数可以解决最小二乘曲线拟合的问题。该函数的调用格式为 [a,Jm,cc,flag,out]=lsqcurvefit(Fun,a0,x,y,options)式中，Fun为原型函数的 MATLAB表示，可以是 M函数或匿名函数；a0为最优化的初值；x和 y为原始输入输出数据向量，对多元函数的拟合而言，这两个输入变元还可以是矩阵，后面将通过例子演示这种情况； options则为最优化工具箱通用的控制模板。调用该函数则将返回待定系数向量 a，以及在此待定系数下的目标函数的值 Jm。如果 flag的值为正，则说明求解成功。 　　例 3-28假设由下面的语句生成一组数据 x和 y >> x=0:0.4:10; %生成样本点数据 y=0.12*exp(-0.213*x)+0.54*exp(-0.17*x).*sin(1.23*x); 并已知该数据满足原型为 y(x)= a1e.a2x + a3e.a4x sin(a5x)，其中， ai为待定系数。采用最小二乘曲线拟合的目的就是获得这些待定系数，使得目标函数的值为最小。　　解显然，这类问题不能由前面介绍的线性回归方法直接求解，只能采用这里介绍的最小二乘曲线拟合方法来求解。 根据已知的函数原型，可以编写出如下的匿名函数。建立起函数的原型，则可以由下面的语句得出待定系数向量为 c = [0.12, 0.213, 0.54, 0.17, 1.23]，拟合残差为 1.7928×10.16。可以看出，得出的待定系数精度较高。下面语句还可以绘制出拟合曲线 与样本点，如图 3-12所示，可见拟合精度很高。 >> f=@(a,x)a(1)*exp(-a(2)*x)+a(3)*exp(-a(4)*x).*sin(a(5)*x); a0=[1,1,1,1,1]; [xx,res]=lsqcurvefit(f,a0,x,y); %最小二乘拟合 x1=0:0.01:10; y1=f(xx,x1); plot(x1,y1,x,y,'o') %拟合效果比较 图 3-12拟合效果比较 其实，由底层命令求解该最优化问题也不困难，只需依赖前面给出的匿名函数，定义一个新的匿名函数描述目标函数，则可以由下面的语句直接求解出所需的待定系数，这样得出的结果与前面的结果完全一致。 >> F=@(a)norm(f(a,x)-y); x1=fminunc(F,a0)　　例 3-29假设有一组实测数据由表 3-5给出，且已知该数据可能满足的原型函数为 y(x)= ax + bx2e.cx + d，试求出满足下面数据的最小二乘解 a、b、c和 d的值。 表 3-5例 3-29实测数据 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 yi 2.3201 2.6470 2.9707 3.2885 3.6008 3.9090 4.2147 4.5191 4.8232 5.1275 　　解表 3-5给出的样本点数据可以通过下面语句直接输入 MATLAB工作空间。 >> x=0.1:0.1:1; %输入样本点数据 y=[2.3201,2.6470,2.9707,3.2885,3.6008,3.9090,... 4.2147,4.5191,4.8232,5.1275]; 令 a1 = a，a2 = b，a3 = c，a4 = d，这样，原型函数可以写成 y(x)= a1x + a2x2e.a3x + a4，可以用匿名函数描述该原型函数。下面语句可以得出函数的待定参数 a = [3.1001, 1.5027, 4.0046, 2]T。注意，本例若不采用循环，则可能收敛不到真值。 >> f=@(a,x)a(1)*x+a(2)*x.^2.*exp(-a(3)*x)+a(4); a=[1;2;2;3]; %原型函数与初值 while (1) [a,f0,cc,flag]=lsqcurvefit(f,a,x,y); if flag>0, break; end, end %最小二乘 ·74·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 用下面的语句还可以计算出各个点处的值，可以将二者曲线绘制在同一坐标系下，如图 3-13所示。可见，二者还是很接近的，说明拟合效果较好。 >> y1=f(a,x); plot(x,y,x,y1,'o') %拟合效果比较 图 3-13拟合效果比较 如果某函数含有若干个自变量，且已知其原型函数 z = f(a,x1,x2, ··· ,xm)，则仍然可以使用 lsqcurvefit()函数来拟合参数 a，其中， a =[a1,a2, ··· ,an]。该函数仍需要用户编写一个匿名函数或 MATLAB函数描述原型函数，然后调用函数 lsqcurvefit()直接求解待定系数向量 a，下面将通过例子演示多变量函数最小二乘拟合的求解方法。 　　例 3-30假设某三元函数的原型函数为 a2xa4(x+y) a6(x+y+z) v = a1x + a3y + a5z 且已知一组输入输出数据，由文本文件 c3data1.dat给出，该文件的前三列为自变量 x、y和 z，第四列为返回向量，试采用拟合方法得出待定系数 ai。　　解解决这类问题第一步仍然需要引入向量型的自变量 x，如令 x1 = x，x2 = y， x3 = z，这样，原型函数可以重新表示为 a2x1 a4(x1+x2) a6(x1+x2+x3) v = a1x + a3x + a5x 12 3 因为给出的数据是纯文本文件，可以通过 load()函数将其读入 MATLAB工作空间，用子矩阵提取的方法将输入矩阵 X和输出向量 v提取出来，可以用下面语句拟合出待定系数的值 a = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6]，使得拟合误差的最小平方和最小，其值为 1.0904×10.7。注意，在匿名函数使用第 i个自变量时，一定要给出 X(:,i)命令来提取该自变量。 >> f=@(a,X)a(1)*X(:,1).^(a(2)*X(:,1))+... a(3)*X(:,2).^(a(4)*(X(:,1)+X(:,2)))+... a(5)*X(:,3).^(a(6)*(X(:,1)+X(:,2)+X(:,3))); %原型函数 XX=load('c3data1.dat'); X=XX(:,1:3); v=XX(:,4); %样本点数据读入 a0=[2 3 2 1 2 3]; [a,f,err,flag]=lsqcurvefit(f,a0,X,v) %最小二乘事实上，文件中给出的数据就是假设 a = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6]生成的，所以用这里给出的拟合方法可以很精确地得出待定系数。 3.5.3边值微分方程的打靶求解 本书卷 V将全面介绍微分方程的解析解与数值解，这里只给出一个简单的例子演示边值微分方程问题的打靶求解方法，并介绍如何将这类方法转换成最优化问题，得出问题的解。 对一般的微分方程而言，需要求解的是初值问题，即已知微分方程的数学描述 x ′ (t)= f(t, x),且已知 x(t0)（3-5-5） 这样，该方程可以调用 ode45()函数直接求出数值解。在实际应用中，如果已知 x(t0)中的几个分量，同时还已知 x(tn)的另外几个值，如何求这里微分方程的数值解呢？打靶方法是常用的求解方法。 打靶法的思路是这样的，先给 x(t0)的另外几个值赋初值，这样可以认为微分方程的初值为已知的，所以可以调用 ode45()求解微分方程，得出终点处微分方程的数值解 x.(tn)。这样，x.(tn)与事先给定 x(tn)中的几个已知的量相比，就存在误差量了。根据该误差量就可以调整 x(t0)的初值，重新求解微分方程。采用这样的方法就可以构造一个循环，对 x.(tn)与 x(tn)中的几个已知的量的误差作最小化，从而得出相容的初值，最终得出微分方程边值问题的数值解。 　　例 3-31假设已知微分方程 . x (t)= x2(t) 1′ . . x ′ 2(t)= .x1(t) . 3x2(t)+ e.5t x (t)= x4(t) 3′ . ′ . x4(t)=2x1(t) . 4x2(t) . 3x3(t) . 4x4(t) . sin t且已知边值 x1(0) = 1，x2(0) = 2，x3(10) = .0.021677，x4(10) = 0.15797，试求解该微分方程在 t → (0, 10)区间内的数值解。　　解可以将 x3(0)和 x4(0)作为决策变量，这样原始问题就可以转换成最优化问题 min ||x3(10) . x.3(10)|| + ||x4(10) . x.4(10)|| x3(0),x4(0) 令 y1 = x3(0)，y2 = x4(0)，则最优化问题的数学形式可以改写成 min ||y1 . x.3(10)|| + ||y2 . x.4(10)|| y 而 x.3(10)和 x.4(10)是以 [x1(0),x2(0),y1,y2]T为初值的微分方程解的后两项终值，所以 ·76·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 目标函数可以如下描述。注意，这里只能采用 M函数描述微分方程，不能采用匿名函数，因为该函数的内部是有中间变量的，不能采用匿名函数。另外， f、x0与 xn是附加参数，后面将演示不采用附加参数的函数格式。 function z=c3mode(y,f,x0,xn) x0=[x0(1:2); y]; [t,x]=ode45(f,[0,10],x0); z=norm([x(end,3:4)]'-xn); 现在可以用匿名函数表示微分方程模型，另外，将目标函数转换成不带有附加参数的匿名函数，可以由下面语句直接求解等效的 x3(0)与 x4(0)，得出 x3(0) = 0.101584621163685，x4(0) = .2.318606769109767。在等效初值下求解微分方程，则可以绘制出微分方程解的曲线，如图 3-14所示。另外可以看出，得出微分方程的终值与给定的初值是完全一致的，说明求解过程是成功的。 >> f=@(t,x)[x(2); -x(1)-3*x(2)+exp(-5*t); x(4); 2*x(1)-3*x(2)-3*x(3)-4*x(4)-sin(t)]; x0=[1; 2]; xn=[-0.021677; 0.15797]; g=@(x)c3mode(x,f,x0,xn); %将带有附加参数的函数转换成普通匿名函数 x2=rand(2,1); x3=fminunc(g,x2) %选择随机初值求解最优化问题 [t,x]=ode45(f,[0,10],[x0; x3]); plot(t,x), x(end,:) 图 3-14微分方程数值解曲线 事实上，这里反推初值的方法得出的解可能是不唯一的。例如，如果初值向量选择为 x3(0) = .1，x4(0) = 1，则得出的微分方程解的终值也满足给定数值的终值条件，如图 3-15所示，新初始条件下微分方程的解如图 3-15所示。在前面求解语句中，每次将 x2设置成不同的随机数向量，得出的相容初值都可能是不同的，由相容初值求解微分方程都得到满足给定的终值条件。 >> [t,x]=ode45(f,[0,10],[x0; -1; 1]); x(end,:) 图 3-15另一组微分方程数值解曲线 3.5.4方程求解问题转换为最优化问题 第 2章已经介绍了很多非线性方程求解的知识。其实，可以很容易地将方程求解问题转换为无约束最优化问题。假设想求解方程 f(x)= 0，那么方程的解 x会满足什么条件呢？当然，方程的解满足使得方程等于 fi(x)=0这样的函数条件，所以可以考虑将这些函数的平方和作为目标函数进行寻优计算，找到决策变量 x。这样，可以写出下面的最优化问题： n ∑ min f12(x)+ f22(x)+ ··· + fn2(x)= min fi 2(x)（3-5-6） x x i=1 　　例 3-32试求解下面给出的二元联立方程组 [4]。 { 32 4x1 +4x1x2 +2x2 . 42x1 = 14 32 4x2 +2x1 +4x1x2 . 26x2 = 22 　　解可以看出，这里的方程是例 2-45欠定方程的变换表示形式。很自然地，可以将方程求解的问题变成如下的最优化问题： 32 32 min (4x1 +4x1x2 +2x2 . 42x1 . 14)2 + (4x2 +2x1 +4x1x2 . 26x2 . 22)2 x 可以先用匿名函数分别定义出两个方程，可以绘制出两个方程的隐函数曲线，如图 3-16所示。从得出的曲线可见，原方程有 9个实根，当然用前面介绍的 more_sols()函数可以立即得出全部的 9个根。 >> f1=@(x1,x2)4*x1.^3+4*x1.*x2+2*x2.^2-42*x1-14; f2=@(x1,x2)4*x2.^3+2*x1.^2+4*x1.*x2-26*x2-22; fimplicit(f1), hold on, fimplicit(f2), hold off 现在考虑方程的最优化解法，借用前面定义的两个匿名函数构造出目标函数，然后选择初值求解，即可以得出方程的一个根。如果改变初值，则可能得出其他的解。 >> f=@(x)(f1(x(1),x(2)))^2+(f2(x(1),x(2)))^2; [x,f0]=fminunc(f,rand(2,1)); ·78·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 图 3-16二元联立方程解的图示 从求解效率和精度看，该函数低于专用于方程求解的 fsolve()函数，这里给出这样的例子主要演示如何将一个一般问题转换为最优化问题，并利用最优化技术求解问题的思想。 从这里给出的例子可以看出，学习了最优化的思想，就可以有意识地在某些问题的研究中应用该思想。在实际应用中比较关键的一步就是设计一个有物理意义的目标函数，并围绕该目标函数探讨最优化问题的求解方法。 本章习题 ∫1 3.1试求出使 (ex . cx)2 dx取极小值的 c值。 3. 2试求解下的无约束最优化问题。 面0 2 22 100(x2 . x1)2 + (1 . x1)2 + 90(x4 . x3) + (1 . x3)2 min [] x + 10.1(x2 . 1)2 +(x4 . 1)2+ 19.8(x2 . 1)(x4 . 1) 3.3试找出下面二元函数曲面的全局谷底。 √ sin 0.1+ (x1 . 4)2 +(x2 . 9)2f(x1,x2)= . 1+(x1 . 4)2 +(x2 . 9)2 3.4试求解 Griewangk(基准测试问题（ (n = 20）。 )) nn ∑ x2 ∏ xi i min 1+ . cos ∞ ,xi → [.600, 600] x 4000 i i=1 i=1 3.5试求解 Ackley基准测试问题 [12]。 (). pp ∑∑ min .20 + 10.20 exp..0.211 x 2 i. . exp 1 cos 2巾xi. x pp i=1 i=1 3.6试求解 Kursawe基准测试问题。 p ∑ J = min |xi|0.8 +5 sin3 xi +3.5828 x i=1 其中，可取 p =2或 p = 20。 3.7试求解扩展的 Freudenstein–Roth函数的最小值问题（n = 20）。 n/2 ∑()2 ()2f(x)= .13+x2i.1+((5.x2i)x2i.2)x2i+.29+x2i.1+((x2i+1)x2i.14)x2i i=1 初始点 x0 = [0.5, .2, ··· , 0.5, .2]T，解析解 xf = [5, 4, ··· , 5, 4]T，fopt =0。如果搜索范围变大，一般求解方法还能否得出问题的全局最优解？ 3.8试求解扩展三角函数的最小值问题（n = 20）。 . 2 nn ∑∑ f(x)= ..n . cos xj. + i(1 . cos xi) . sin xi. i=1 j=1 初始点 x0 = [1/n, 1/n, ··· , 1/n]T，解析解 xi =0，fopt =0。 3.9试求解扩展 Rosenbrock函数的最小值问题（n = 20）。 n/2 ∑()2 2 f(x) = 100x2i . x2i.1+ (1 . x2i.1)2i=1 初始点 x0 =[.1.2, 1, ··· , .1.2, 1]，解析解 xi =1，fopt =0。如果将 x2i . x22i.1替换成 x2i . x32i.1，则变成扩展 White–Holst问题，试求解该问题。 3.10试求解扩展 Beale函数的最小值问题（n = 20）。 n/2 ∑[]2 []2 []2 23 f(x)= 1.5.x2i.1(1.x2i)+2.25.x2i.1(1.x2i)2+2.625.x2i.1(1.x2i)2 i=1 初始值 x0 = [1, 0.8, ··· , 1, 0.8]T，解析解未知。 3.11试求解两个 Raydan函数的最小值问题（n = 20）。 nn ∑ i ()∑() exi f1(x)= exi . xi,f2(x)= . xi 10 i=1 i=1 n ∑ 初始值 x0 = [1, 1, ··· , 1]T，解析解 xi =1，f1opt = i/10，f2opt = n。 i=1 3.12试求解扩展 Miele–Cantrell函数的最小值问题（n = 20）。 n/4 ∑()2 8 f(x)= e4i.3 . x4i.2+ 100(x4i.2 . x4i.1)6 + tan4(x4i.1 . x4i)+ x 4i.3i=1 初始值 x0 =[1, 2, 2, 2, ···, 1, 2, 2, 2]T，解析解 xf =[0, 1, 1, 1, ···, 0, 1, 1, 1]T，fopt =0。 3. 13上面很多习题都假设 n = 20，如果 n = 200，试重新求解这些习题。 3. 14对下面给出的 Hartmann函数，试求解全局小值问题。 . 43 ∑∑ f(x)= . 向i exp .. aij(xj . pij )2. i=1j=1 ·80·薛定宇教授大讲堂（卷 IV）：MATLAB最优化计算 其中，α = [1, 1, 2, 3, 3, 2]T，且 . 3 10 30 3689 1170 2673 0.1 10 35 4699 4387 7470 A = , P = 10.4 × 3 10 30 1091 8732 5547 0.1 10 35 381.5 5743 8828 3.15试求解 Schwefel函数的全局最小化问题（n = 20）。 .. .. .. . .. n f(x) = 418.9829n . i=1 ∑ ..搜索区域 解析解 。.500500=1f=0xxopt，，ii √ ||xi xi sin 并绘制出 n =2时目标函数 的曲面。 3.16Eggholder试求函数的全局最小值。 √ sin.f()= (+47) x,y y x +(y + 47) . x sin 2 √ x . (y + 47) ) ( ) ( 并绘制出目标函数的曲面。 3.17试求出下面函数在 .2 . x . 11区域内的最小值，并用图形法检测结果。 52397179 1 54 f(x)= x 6 . x + x + x 3 . x 2 . x + 3. 18试分别求出函数 f(x)= 10, 20)内下面函数的最大值，并绘制图形验证其结果。 xsin2510巾x+802在(.110,2)和(.20区域10 3.19试求出 .10 . x, y . 10区域内下面函数的最小值，并用图形显示这样的解。 f(x, y)= sin2 3巾x +(x . 1)21+ sin2 3巾y+(y . 1)21+ sin2 2巾y 3.20例 3-20方程由许多波谷，如果仔细观察会发现谷底的函数值可能不同，试用图解法找出全局最优解（提示：考虑绘制侧视图与正视图，找出 x1和 x2的最优解）。 3.21如果 0 . x, y . 1，试求出下面的函数的最大值，并绘制图形验证结果。 f(x, y)= sin 19巾x +1x.7+ sin 19巾y +1y.7 +2 ) ( 3.22试求 .100 . x, y . 100区域内下面函数的最小值。 sin |x2 . y2|1+0.001(x2 + y2) cos2 . 0.5 [ f(x, y)=0.5+ 3.23假设某一化学过程中压力 P与温度 T之间的关系为 P ]2 =向e.T ，但向与 .是未知 的。现通过实验测得一些数据，在表 3-6中给出，试确定未知参数向与 . [13]。表 3-6习题 3.23中已知的实验数据 温度 T / .C 20253035405060 70压力 P / mmHg 15.45 19.23 26.54 34.52 48.32 68.11 98.34 120.45 3.24假设已知一组数据如表 3-7所示，且已知该数据满足原型函数。 1 e.(x.μ)2/2β2 y(x)= ∞ 2巾θ 试用最小二乘法求出 μ和 θ的值，并用得出的函数将函数曲线绘制出来，观察曲线拟合的效果。 表 3-7习题 3.24数据 xi .2 .1.7 .1.4 .1.1 .0.8 .0.5 .0.2 0.1 0.4 0.7 1 1.3 yi 0.1029 0.1174 0.1316 0.1448 0.1566 0.1662 0.1733 0.1775 0.1785 0.1764 0.1711 0.1630 xi 1.6 1.9 2.2 2.5 2.8 3.1 3.4 3.7 4 4.3 4.6 4.9 yi 0.1526 0.1402 0.1266 0.1122 0.0977 0.0835 0.0702 0.0577 0.0469 0.0373 0.0291 0.0224 3.25假设有一组数据在文件 data3ex5.dat中给出，其第 1列为散点的 x坐标，第 2列为 y坐标，第 3列为函数值，试由文件的数据计算下面原型函数的待定系数 a~e。 2 2 . bx)e.cx .dy2.exy f(x, y)=(ax 3. 26已知微分方程模型 x ′ 1(t)= x2(t)，x ′ 2(t)=2x1(t)x2(t)，且已知边值条件 x1(0) = .1，x1(巾/2) = 1，试求解相应的微分方程数值解。 3. 27已知某常微分方程模型如下，试求出向和 .，并求解本微分方程， x ′ 1 =4x1 .向x 1x2，x ′ = .2x2 + .x1x2，且已知 x1(0) = 2，x2(0) = 1，x1(3) = 4，x2(3) = 2。 3.28例 3-32给出了将方程求解问题转换为最优化问题的方法，试运行该函数 100次，并运行 fsolve()函数 100次解决同样的问题，比较方程求解的效率与精度。 中2 3.29试将习题 2.12、2.14中给出的方程求解问题转换成最优化问题并求解方程，试比较方程求解的精度与速度。