一、函数和子函数
一个M文件中,可能会有多个函数,其中第一个称为主函数,后面的所有函数称为子函数
- 脚本文件中,也可以直接在脚本的最后添加子函数,在当前文件夹内,如果有同名函数,按照子函数
MATLAB内建函数其他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)
end
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”!')
end
shg
end
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
可使用字符串,匿名函数或函数句柄。返回值列表为极小值点,fval
为函数极小值,exitflag>0
代表此函数成功找到了一个极值点。其余功能output
,option
可以通过帮助系统了解用法,它们主要可以让MATLAB显示迭代过程中的各种运算指标。[x,fval,exitflag]=fminsearch(fun,x0)
可以利用无导数方法(如单纯形法),从点出发,求多元函数fun
在在多维空间中的一个极小值,fun
建议使用多元匿名函数或多元函数句柄(输入值为向量)。返回值列表为极小值点向量,fval
为函数极小值,其余功能与fminbnd
类似。- 注意到两个函数均为找一个极小值,有时也称为局部最小值,并不一定是定义域内的全局最小值。希望获取全局最小值,需要设立较好的区间或起始点,或反复遍历所有极值点。
1. 划分区间求极值
例:求
在
的极小值
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 %若有更小的极小值则更换最小值
end
fprintf('函数最小值%6.5f发生在x=%6.5f处',fc,xc)
2. 单峰函数求极值之黄金分割法
对于向内的试探法,从
开始,当
时,探索的效率最高,这样的试探法就成为黄金分割法(0.618法)
黄金分割法的算法复杂度为
,而且可以处理目标函数不可导的特殊情况。操作代码如下:
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);
end
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);
end
else
x2 = alpha2;
if(yx(alpha1)<minf)
minf = yx(alpha1);
end
end
if(alpha2-alpha1<1e-8)
break;
end
end
minf, alpha1
3. 牛顿迭代法
- 牛顿迭代法起源于一阶泰勒展开近似
此时, 设存在满足,则给定任意一点有移项后即可得到对于线性函数此法可快速求解
- 几何上讲, 牛顿迭代法类似于寻找切线方向并移动至x轴交点, 迭代公式记为:
- 对于求局部最小值问题, 因
因此对应的牛顿迭代法公式只需改为导数方程问题,迭代公式:
例:求
在
的最小值,代码如下:
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
end
x_old = x_new;
end
disp([iter, x_new, yx(x_new)]);
三、非线性方程组Matlab求解
求解非线性方程或方程组,除了单调函数二分法、光滑函数牛顿法外,还可以使用MATLAB函数进行求解,主要有以下两种形式:
[x,fval]=fzero(fun,x0)
将以为初值,尝试寻找函数句柄或匿名函数fun
的一个零点,fval
为对应函数值[x,fval]=fsolve(fun,x0)
将以为初值(向量),尝试寻找函数向量fun
的一个零点,fval
为对应函数向量值
例:求
的零点,代码如下:
% 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 %回到正常比例
四、多元函数最优化
例:求
的极小值点,代码如下:
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. 梯度下降法求解多元最小值问题
将一元最优化问题推广到多元问题
,设
满足一定的光滑性且其梯度除极值点外不为0向量。则
在
下降最快的方向为其负梯度方向:
。
因此对于凸问题,只需按
进行迭代即可无限接近理论最小值,其中
为待定步长。对于非凸问题,迭代有可能会收敛到某个
局部最小值。
对于多元情形:
- 多元问题同样存在最速下降法,不过有时由于最佳步长计算复杂,也可能用固定步长
代替。
- 多元问题是否为凸问题,等价于目标函数
的Hesse矩阵是否半正定的判定问题。
例:用梯度下降法求
的极小值点,代码如下:
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均保持稳定时结束迭代
end
x_old = x_new;y_old = y_new;
end
iter,x_new,y_new
err = ff(x_new,y_new)-0
本题使用的最速梯度下降法,虽然速度较慢,但在迭代的收敛性和准确性上,均优于基于单纯形法的MATLAB函数fminsearch
,速度上不如新函数fminunc
2. 牛顿法解多元最小值问题
多元问题
在使用最速梯度下降法求解时,由于相邻迭代的梯度方向正交,导致移动方向呈现锯齿型,收敛速度与效率不佳。且对于病态矩阵无健壮性。
定义Hesse矩阵 :
多元问题
的牛顿迭代公 式为 (注意矩阵相乘后方向可能改变) :
牛顿法的一个缺点即Hessian不能出现不可逆的点。但此方法既可以增加收敛的概率,也可以大大减少运行时间。
例:用牛顿法求
的极小值点,代码如下:
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;
p = hidff(x_old,y_old);
x_new=x_old-p(1); y_new=y_old-p(2);
if(abs(x_new-x_old)<1e-8 && abs(y_new-y_old)<1e-8) break; end
x_old = x_new;y_old = y_new;
end
iter,x_new,y_new