Matlab问题求解
一、函数和子函数
一个M文件中，可能会有多个函数，其中第一个称为主函数，后面的所有函数称为子函数
-
脚本文件中，也可以直接在脚本的最后添加子函数，在当前文件夹内，如果有同名函数，按照子函数
\rightarrow
MATLAB内建函数
\rightarrow
其他
\rightarrow
M文件主函数的顺序访问。子函数最后的
end不能省略
-
一个M文件的主函数通常和M文件名相同（否则MATLAB仍以文件名主名作为识别标准），一个M包含多个函数时，每个函数最后的
end或者都省略掉，或者都不省略。
-
所有的子函数都可以被M文件内的脚本或主函数调用，但无法被其他M文件或命令行直接调用。因此，子函数是一种减少M文件数量，封装外部脚本不直接调用的函数的好方法。
1. 子函数
编写一个内含子函数的M函数绘图文件
HC = Draw_d('circle');
HL = Draw_d('line');
function Hr = Draw_d(flag)
% exm060301.m Demo for handles of primary functions and subfunctions
% flag %允许使用字符串’line’或’circle’
% Hr %返回子函数cirline的句柄
t = (0:50)/50*2*pi; %0~2pi等分了50个区间
x = sin(t);
y = cos(t);
Hr = @cirline; %创建cirline的句柄（一种函数名的不同解读，类似于C++的指针）
feval(Hr, flag, x, y, t); %feval结合句柄调用，等价于cirline(flag,x,y,t)
function cirline(wd, x, y, t)
% wd %主函数传递来的flag，可能为’line’或’circle’
% t,x,y %分别为绘图参数、横坐标与纵坐标
switch wd
case 'line'
plot(t, x, 'b', t, y, 'r', 'LineWidth', 2);
case 'circle'
plot(x, y, '-g', 'LineWidth', 8);
axis square off;
otherwise
error('输入变量只能取“line”或“circle”！')
HC
的输出结果为：
HL
的输出结果为：
另外我们可以将
t
的采样距离缩小，比如绘制正五边形采样点分布：
t=0:2*pi/5:2*pi;x=cos(t);y=sin(t);
HH('circle',x,y,t)
%利用m文件主函数返回的句柄可以间接调用到子函数。
%如果没有返回的句柄，则子函数cirline无法被外部调用
2. 匿名函数
- 适用于结构简单、无需创建m文件或子函数体来定义的“一句话函数”
-
例如：
F=@(x) sin(x).*x,就定义了一个自变量为x，函数值F(x)=sin(x).*x的匿名函数。命名上，F称为函数句柄，括号内的x称为参数（参数可以由多个变量构成参数列表），sin(x).*x称为函数主体表达式。 -
匿名函数可以通过如
F(3)或F(x)直接调用，也可以通过feval(F,x)间接调用,有时，函数句柄F还可以作为参数代入一些更为复杂的函数体。
3. 函数句柄
函数的句柄类似于指针，除
plot
绘图返回值，匿名函数定义外，还可对MATLAB内建函数或用户已定义函数创建句柄。
hm = @magic;
class(hm) % ans = 'function_handle'
functions(hm) % 查询句柄对应函数信息，包含function（函数名）, type（函数类型sinple）, file（文件位置）
hm(4)
16 2 3 13
5 11 10 8
9 7 6 12
4 14 15 1
调用子函数时，函数句柄可以“完全的代替”本身子函数的函数名，格式与子函数直接调用相同。但利用函数句柄多次调用可以大大节省调用的时间（不必重复搜索路径）。
二、函数极值的数学方法
-
[x,fval,exitflag]=fminbnd(fun,x1,x2)可以求一元函数fun在[x1,x2]的一个极小值，fun可使用字符串，匿名函数或函数句柄。返回值列表 x 为极小值点，fval为函数极小值，exitflag>0代表此函数成功找到了一个极值点。其余功能output，option可以通过帮助系统了解用法，它们主要可以让MATLAB显示迭代过程中的各种运算指标。
-
[x,fval,exitflag]=fminsearch(fun,x0)可以利用无导数方法（如单纯形法），从 x_0 点出发，求多元函数fun在在多维空间中的一个极小值，fun建议使用多元匿名函数或多元函数句柄（输入值为向量）。返回值列表 x 为极小值点向量，fval为函数极小值，其余功能与fminbnd类似。
-
注意到两个函数均为找一个极小值，有时也称为局部最小值，并不一定是定义域内的全局最小值。希望获取全局最小值，需要设立较好的区间或起始点，或反复遍历所有极值点。
1. 划分区间求极值
例：求 y=e^{-0.1x}\text{sin}^2x-0.5(x+0.1)\text{sin}x 在 -50\le x\le 50 的极小值
x1=-50;x2=5;
yx=@(x)(sin(x).^2.*exp(-0.1*x)-0.5*sin(x).*(x+0.1)); %定义函数句柄
[xc0,fc0,exitflag]=fminbnd(yx,x1,x2) %找到的其中一个极小值
ezplot(yx,[-50,5]); %ezplot同样可用于函数句柄，只需指定x的范围即可
xlabel('x'),grid on
xx=[-23,-20,-18]; %观察最小值疑似存在的两个区段
fc=fc0;xc=xc0; %暂时设立最小值点和最小值为初始搜索的结果
for k=1:2
[xw,fw]=fminbnd(yx,xx(k),xx(k+1)); %分别计算两个区段的极小值
if fw<fc, xc=xw;fc=fw;end %若有更小的极小值则更换最小值
fprintf('函数最小值%6.5f发生在x=%6.5f处',fc,xc)
2. 单峰函数求极值之黄金分割法
对于向内的试探法，从 [a_k,b_k] 开始，当 α_1≈0.618a_k+0.382b_k,α_2≈0.382a_k+0.618b_k 时，探索的效率最高，这样的试探法就成为黄金分割法（0.618法）
黄金分割法的算法复杂度为 O(\text{ln}\frac{b-a}{\epsilon}) ，而且可以处理目标函数不可导的特殊情况。操作代码如下：
x1=-20;x2=-18;
yx=@(x)(sin(x).^2.*exp(-0.1*x)-0.5*sin(x).*(x+0.1));
%设立初始最小值
minf = yx(x1);
if(yx(x2)<minf)
minf = yx(x2);
while(1)
alpha1 = x1*0.618+x2*0.382;
alpha2 = x1*0.382+x2*0.618;
if(yx(alpha1)>yx(alpha2))
x1 = alpha1;
if(yx(alpha2)<minf)
minf = yx(alpha2);
x2 = alpha2;
if(yx(alpha1)<minf)
minf = yx(alpha1);
if(alpha2-alpha1<1e-8)
break;
minf, alpha1
3. 牛顿迭代法
-
牛顿迭代法起源于一阶泰勒展开近似
f(x) \approx f\left(x_{0}\right)+f^{\prime}\left(x_{0}\right)\left(x-x_{0}\right) \\ 此时, 设存在 \bar{x}, 满足 f(\bar{x})=0 ,则给定任意一点 x, 有
f(x) \approx f(\bar{x})+f^{\prime}(\bar{x})(x-\bar{x})=f^{\prime}(\bar{x})(x-\bar{x}) \\ 移项后即可得到
\bar{x} \approx x-\frac{f(x)}{f^{\prime}(x)} \\ 对于线性函数 f(x), 此法可快速求解 f(x)=\mathbf{0}
-
几何上讲, 牛顿迭代法类似于寻找切线方向并移动至x轴交点, 迭代公式记为：
x_{n+1} = x_n-\frac{f(x_n)}{f^`(x_n)} \\ -
对于求局部最小值问题, 因
\text{argmin }_x f(x) \Leftrightarrow f^{\prime}(x)=0 \\ 因此对应的牛顿迭代法公式只需改为导数方程问题，迭代公式：
x_{n+1}=x_{n}-\frac{f^{\prime}\left(x_{n}\right)}{f^{\prime \prime}\left(x_{n}\right)} \\
例：求 y=e^{-0.1x}\text{sin}^2x-0.5(x+0.1)\text{sin}x 在 -20\le x\le -18 的最小值，代码如下：
x1=-20;x2=-18;syms x;
yx=@(x)(sin(x).^2.*exp(-0.1*x)-0.5*sin(x).*(x+0.1));
yxp = diff(yx,x); %MATLAB默认使用sym类型储存导函数表达式
yxp2 = diff(yxp,x); %二阶导函数表达式
x_old = -19;iter = 0; %迭代初始值与迭代次数
while(1)
iter = iter + 1;
x_new = x_old - double(subs(yxp, x, x_old) / subs(yxp2, x, x_old));
if(abs(x_new-x_old)<1e-8)
break
x_old = x_new;
disp([iter, x_new, yx(x_new)]);
三、非线性方程组Matlab求解
求解非线性方程或方程组，除了单调函数二分法、光滑函数牛顿法外，还可以使用MATLAB函数进行求解，主要有以下两种形式：
-
[x,fval]=fzero(fun,x0)将以 x_0 为初值，尝试寻找函数句柄或匿名函数fun的一个零点，fval为对应函数值
-
[x,fval]=fsolve(fun,x0)将以 x_0 为初值(向量)，尝试寻找函数向量fun的一个零点，fval为对应函数向量值
例：求 f(t)=sin^2t·e^{-0.1t}-0.5|t| 的零点，代码如下：
% solve函数求解
syms t;
ft=sin(t)^2*exp(-0.1*t)-0.5*abs(t);
S = solve(ft, t);
ftS=subs(ft,t,S);
disp([S, ftS]);
% 作图观察
y_C=@(t) sin(t).^2.*exp(-0.1*t)-0.5*abs(t); %函数句柄形式的定义
t=-10:0.01:10; Y=y_C(t); %句柄很多时候可以跟函数名一样使用
clf,plot(t,Y,'r');hold on
plot(t,zeros(size(t)),'k'); %另画一条黑色0函数曲线
xlabel('t');ylabel('y(t)')
hold off
% 观察零点位置
zoom on %鼠标拉出方框，将方框选定区域放大
[tt,yy]=ginput(5); %鼠标点击五个位置（零点近似值），并记录坐标
zoom off %回到正常比例
四、多元函数最优化
例：求 f(x,y)=100(y-x^2)^2+(1-x)^2 的极小值点，代码如下：
ff=@(x)(100*(x(2)-x(1)^2)^2+(1-x(1))^2);
%函数句柄，x为输入的(行或列)向量,利用两个元素分别进行计算
syms x y,ezsurfc(ff([x,y]),[-2,2,-2,2]) %将横纵坐标x,y认定为ff二维自定义变量，即可进行surfc操作
format short g %五位有效数字，省略小数点后尾数的0
x0=[-5,-2,2,5;-5,-2,2,5]; %设立4种不同的搜索起点（每一种为列向量）
[sx,sfval,sexit,soutput]=fminsearch(ff,x0)
%收敛到了四种不同的解，但仅有第一个x=1,y=1是正确的
format short e %短科学计数法
disp([ff(sx(:,1)),ff(sx(:,2)),ff(sx(:,3)),ff(sx(:,4))])
%比较四个极小值点对应的函数值，显然第一个点为四个极小值点中取最小值的
clear,clc
ff=@(x) (100*(x(2)-x(1)^2)^2+(1-x(1))^2); %匿名函数需以向量为变量
x0=[5,5]; %仅接受一种初值计算，以5,5为例
options = optimoptions(@fminunc,'OptimalityTolerance',1e-8);
%设置迭代的终止条件，以确保误差小于1e-8（默认1e-6）
[x,fval] = fminunc(ff,x0,options)
%求解（前述四种初值点均可以成功找到准确的解）
1. 梯度下降法求解多元最小值问题
将一元最优化问题推广到多元问题 \text{argmin}_{x\in \mathbb{R}^n}\phi(x) ，设 \phi(x) 满足一定的光滑性且其梯度除极值点外不为0向量。则 \phi(x) 在 x=x^{(k)} 下降最快的方向为其负梯度方向： \boldsymbol{p}^{(k)}=-\nabla \varphi\left(x^{(k)}\right) 。
因此对于凸问题，只需按 x^{(k+1)}=x^{(k)}+\lambda p^{(k)} 进行迭代即可无限接近理论最小值，其中 \lambda>\mathbf{0} 为待定步长。对于非凸问题，迭代有可能会收敛到某个 局部最小值 。
对于多元情形：
-
多元问题同样存在最速下降法，不过有时由于最佳步长计算复杂，也可能用固定步长
λ
代替。
-
多元问题是否为凸问题，等价于目标函数
φ(x)
的Hesse矩阵是否半正定的判定问题。
例：用梯度下降法求 f(x,y)=100(y-x^2)^2+(1-x)^2 的极小值点，代码如下：
clear,clc
ff=@(x,y) (100*(y-x^2)^2+(1-x)^2);%函数及其导数
dff =@(x,y) [2*x - 400*x*(- x^2 + y) - 2;- 200*x^2 + 200*y];
x0 = -5; y0 = -2; %初始条件（可改变）
x_old = x0;y_old = y0;iter=0;
while(1)
iter = iter+1;
Grad = dff(x_old,y_old); %梯度方向的获得
lsf = @(lambda) ff(x_old-lambda*Grad(1),y_old-lambda*Grad(2)); %生成对应方向关于步长的一元函数
[lambda,~]=fminbnd(lsf,0,10);%搜索最佳步长
x_new=x_old-lambda*Grad(1);
y_new=y_old-lambda*Grad(2);
if(abs(x_new-x_old)<1e-8 && abs(y_new-y_old)<1e-8) break; %当x与y均保持稳定时结束迭代
x_old = x_new;y_old = y_new;
iter,x_new,y_new
err = ff(x_new,y_new)-0
本题使用的最速梯度下降法，虽然速度较慢，但在迭代的收敛性和准确性上，均优于基于单纯形法的MATLAB函数
fminsearch
，速度上不如新函数
fminunc
2. 牛顿法解多元最小值问题
多元问题 \text{argmin}_{x\in \mathbb{R}^n}\phi(x) 在使用最速梯度下降法求解时，由于相邻迭代的梯度方向正交，导致移动方向呈现锯齿型，收敛速度与效率不佳。且对于病态矩阵无健壮性。
定义Hesse矩阵 :
\nabla^{2} f(x, y)=\left[\begin{array}{ll} \frac{\partial^{2} f}{\partial x^{2}} & \frac{\partial^{2} f}{\partial x \partial y} \\ \frac{\partial^{2} f}{\partial x \partial y} & \frac{\partial^{2} f}{\partial y^{2}} \end{array}\right] \\
多元问题 argmin_{x \in \mathbb{R}^{n}} \varphi(x) 的牛顿迭代公 式为 (注意矩阵相乘后方向可能改变) :
x^{(k+1)}=x^{(k)}-\left.\left(\nabla^{2} f(x, y)\right)^{-1} \nabla f(x, y)\right|_{x^{(k)}} \\
牛顿法的一个缺点即Hessian不能出现不可逆的点。但此方法既可以增加收敛的概率，也可以大大减少运行时间。
例：用牛顿法求 f(x,y)=100(y-x^2)^2+(1-x)^2 的极小值点，代码如下：
clear,clc
ff=@(x,y) (100*(y-x^2)^2+(1-x)^2);
hidff =@(x,y) [ (x - 1)/(200*x^2 - 200*y + 1), -(200*x^4 - 400*x^2*y - x^2 + 2*x + 200*y^2 - y)/(200*x^2 - 200*y + 1)];
%提前计算的迭代步长，即inv(Hessian)*Gradient
x0 = -5; y0 = -5;x_old = x0;y_old = y0;iter=0;
while(1)
iter = iter+1;