matlab/Matlab问题求解

函数和子函数

一个M文件中,可能会有多个函数,其中第一个称为主函数,后面的所有函数称为子函数

  • 脚本文件中,也可以直接在脚本的最后添加子函数,在当前文件夹内,如果有同名函数,按照子函数$\rightarrow$MATLAB内建函数$\rightarrow$其他$\rightarrow$M文件主函数的顺序访问。子函数最后的end不能省略

  • 一个M文件的主函数通常和M文件名相同(否则MATLAB仍以文件名主名作为识别标准),一个M包含多个函数时,每个函数最后的end或者都省略掉,或者都不省略。

  • 所有的子函数都可以被M文件内的脚本或主函数调用,但无法被其他M文件或命令行直接调用。因此,子函数是一种减少M文件数量,封装外部脚本不直接调用的函数的好方法。

子函数

编写一个内含子函数的M函数绘图文件

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
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的采样距离缩小,比如绘制正五边形采样点分布:

1
2
3
4
t=0:2*pi/5:2*pi;x=cos(t);y=sin(t);
HH('circle',x,y,t)
%利用m文件主函数返回的句柄可以间接调用到子函数。
%如果没有返回的句柄,则子函数cirline无法被外部调用

匿名函数

  • 适用于结构简单、无需创建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还可以作为参数代入一些更为复杂的函数体。

函数句柄

函数的句柄类似于指针,除plot绘图返回值,匿名函数定义外,还可对MATLAB内建函数或用户已定义函数创建句柄。

1
2
3
4
hm = @magic;
class(hm) % ans = 'function_handle'
functions(hm) % 查询句柄对应函数信息,包含function(函数名), type(函数类型sinple), file(文件位置)
hm(4)
1
2
3
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代表此函数成功找到了一个极值点。其余功能outputoption可以通过帮助系统了解用法,它们主要可以让MATLAB显示迭代过程中的各种运算指标。

  • [x,fval,exitflag]=fminsearch(fun,x0)可以利用无导数方法(如单纯形法),从$x_0$点出发,求多元函数fun在在多维空间中的一个极小值,fun建议使用多元匿名函数或多元函数句柄(输入值为向量)。返回值列表$x$为极小值点向量,fval为函数极小值,其余功能与fminbnd类似。

  • 注意到两个函数均为找一个极小值,有时也称为局部最小值,并不一定是定义域内的全局最小值。希望获取全局最小值,需要设立较好的区间或起始点,或反复遍历所有极值点。

划分区间求极值

例:求$y=e^{-0.1x}\text{sin}^2x-0.5(x+0.1)\text{sin}x$ 在 $-50\le x\le 50$的极小值

image-20201117170312325

1
2
3
4
5
6
7
8
9
10
11
12
13
14
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)

单峰函数求极值之黄金分割法

对于向内的试探法,从$[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})$,而且可以处理目标函数不可导的特殊情况。操作代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
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

牛顿迭代法

  • 牛顿迭代法起源于一阶泰勒展开近似

    此时, 设$\exist\bar{x},$ 满足 $f(\bar{x})=0$,则给定任意一点 $x,$有

    移项后即可得到

    对于线性函数 $f(x),$ 此法可快速求解 $f(x)=\mathbf{0}$

  • 几何上讲, 牛顿迭代法类似于寻找切线方向并移动至x轴交点, 迭代公式记为:

  • 对于求局部最小值问题, 因

    因此对应的牛顿迭代法公式只需改为导数方程问题,迭代公式:

例:求$y=e^{-0.1x}\text{sin}^2x-0.5(x+0.1)\text{sin}x$ 在 $-20\le x\le -18$的最小值,代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
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)将以$x_0$为初值,尝试寻找函数句柄或匿名函数fun的一个零点,fval为对应函数值

  • [x,fval]=fsolve(fun,x0)将以$x_0$为初值(向量),尝试寻找函数向量fun的一个零点,fval为对应函数向量值

例:求 $f(t)=sin^2t·e^{-0.1t}-0.5|t|$ 的零点,代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
% 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$的极小值点,代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
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)
%求解(前述四种初值点均可以成功找到准确的解)

梯度下降法求解多元最小值问题

将一元最优化问题推广到多元问题$\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$的极小值点,代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
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

牛顿法解多元最小值问题

多元问题$\text{argmin}_{x\in \mathbb{R}^n}\phi(x)$在使用最速梯度下降法求解时,由于相邻迭代的梯度方向正交,导致移动方向呈现锯齿型,收敛速度与效率不佳。且对于病态矩阵无健壮性。

定义Hesse矩阵 :

多元问题 $argmin_{x \in \mathbb{R}^{n}} \varphi(x)$ 的牛顿迭代公 式为 (注意矩阵相乘后方向可能改变) :

牛顿法的一个缺点即Hessian不能出现不可逆的点。但此方法既可以增加收敛的概率,也可以大大减少运行时间。

例:用牛顿法求$f(x,y)=100(y-x^2)^2+(1-x)^2$的极小值点,代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
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
  • 版权声明: 本博客所有文章除特别声明外,著作权归作者所有。转载请注明出处!

扫一扫,分享到微信

微信分享二维码
  • Copyrights © 2020-2021 chenk
  • 由 帅气的CK本尊 强力驱动
  • 访问人数: | 浏览次数:

请我喝杯咖啡吧~

支付宝
微信