【学习笔记】【算法】【智能优化】粒子群优化(PSO)

【学习笔记】【算法】【智能优化】粒子群优化(PSO)

文章目录

课程地址:https://www.bilibili.com/video/BV1uY41187rK?spm_id_from=333.999.0.0

1 算法背景

1.1 背景

粒子群算法,也称粒子群优化算法或鸟群觅食算法(Particle Swarm Optimization),缩写为PSO。

粒子群优化算法是一种进化计算技术(evolutionary computation),1995年由 Eberhart 博士和 kennedy 博士提出,源于对鸟群捕食的行为研究。

鸟群捕食时,会通过个体对信息的共享使群体找到最优的目的地-食物量最多的位置。我们可以设想一下鸟群捕食的场景:鸟群中的鸟随机散落在森林各处搜索食物,它们的目标是搜索到食物量最多的位置。但是所有的鸟都不知道食物量最多的位置具体在森林中的哪个位置,只能猜测食物大概在哪个方向。每只鸟沿着自己判定的方向进行搜索,并在搜索的过程中记录自己曾经找到过食物量最多的位置,同时所有的鸟都共享自己曾经找到过食物量最多的位置以及食物的量,这样鸟群就知道当前在哪个位置食物的量最多。在搜索的过程中每只鸟都会根据自己记忆中食物量最多的位置和当前鸟群记录的食物量最多的位置调整自己接下来搜索的方向。鸟群经过一段时间的搜索后就可以找到森林中哪个位置的食物量最多(全局最优解)。

该算法最初是受到飞鸟集群活动的规律性启发,进而利用群体智能建立的一个简化模型。粒子群算法在对动物集群活动行为观察基础上,利用群体中的个体对信息的共享使整个群体的运动在问题求解空间中产生从无序到有序的演化过程,从而获得最优解。

1.2 基础知识

问:已知A为全局最优,B和C如何移动才能到达A处?

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-K3b8LPyp-1662473600279)(pic1.png)]

这个过程如何用数学表达式描述?

  1. 某个粒子(点)的移动,是有大小,有方向的。

  2. 有大小,有方向的东西叫向量。

  3. 位置就是坐标。

( 1 , 1 ) = ( 2 , 3 ) + a → a = ( 1 , 1 ) − ( 2 , 3 ) = ( − 1 , − 2 ) (1,1)=(2,3)+a→ a=(1,1)-(2,3)=(-1,-2) (1,1)=(2,3)+aa=(1,1)(2,3)=(1,2)

( 1 , 1 ) = ( 2 , 3 ) + ( − 1 , − 2 ) (1,1)=(2,3)+(-1,-2) (1,1)=(2,3)+(1,2)

( x , y ) = ( 2 , 3 ) + r a n d ∗ ( − 1 , − 2 ) (x,y)=(2,3)+rand * (-1,-2) (x,y)=(2,3)+rand(1,2)

2 算法原理

2.1 基本原理

假设在一个D维的目标搜索空间中,有N个粒子组成一个群落,其中第 i 个粒子表示为一个 D 维的向量:
X i = ( x i 1 , x i 2 , ⋅ ⋅ ⋅ , x i D ) , i = 1 , 2 , ⋅ ⋅ ⋅ , N % MathType!MTEF!2!1!+- % feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn % hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr % 4rNCHbWexLMBbXgBd9gzLbvyNv2CaeHbl7mZLdGeaGqiVu0Je9sqqr % pepC0xbbL8F4rqqrFfpeea0xe9Lq-Jc9vqaqpepm0xbba9pwe9Q8fs % 0-yqaqpepae9pg0FirpepeKkFr0xfr-xfr-xb9adbaqaaeGaciGaai % aabeqaamaabaabauaakeaacaWGybWaaSbaaSqaaiaabMgaaeqaaOGa % eyypa0ZaaeWaaeaacaqG4bWaaSbaaSqaaiaabMgacaaIXaaabeaaki % aamYcacaqG4bWaaSbaaSqaaiaabMgacaaIYaaabeaakiaamYcacqGH % flY1cqGHflY1cqGHflY1caaJSaGaaeiEamaaBaaaleaacaqGPbGaam % iraaqabaaakiaawIcacaGLPaaacaaJSaGaaWyAaiaam2dacaaJXaGa % aWilaiaamkdacaaJSaGaeyyXICTaeyyXICTaeyyXICTaaWilaiaam6 % eaaaa!6335! {X_{\rm{i}}} = \left( {{{\rm{x}}_{{\rm{i}}1}},{{\rm{x}}_{{\rm{i}}2}}, \cdot \cdot \cdot ,{{\rm{x}}_{{\rm{i}}D}}} \right),i = 1,2, \cdot \cdot \cdot ,N Xi=(xi1,xi2,,xiD),i=1,2,,N
第 i 个粒子的“飞行”速度也是一个 D 维的向量,记为:
V i = ( v i 1 , v i 2 , ⋅ ⋅ ⋅ , v i D ) , i = 1 , 2 , ⋅ ⋅ ⋅ , N % MathType!MTEF!2!1!+- % feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn % hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr % 4rNCHbWexLMBbXgBd9gzLbvyNv2CaeHbl7mZLdGeaGqiVu0Je9sqqr % pepC0xbbL8F4rqqrFfpeea0xe9Lq-Jc9vqaqpepm0xbba9pwe9Q8fs % 0-yqaqpepae9pg0FirpepeKkFr0xfr-xfr-xb9adbaqaaeGaciGaai % aabeqaamaabaabauaakeaacaWGwbWaaSbaaSqaaiaabMgaaeqaaOGa % eyypa0ZaaeWaaeaacaqG2bWaaSbaaSqaaiaabMgacaaIXaaabeaaki % aamYcacaqG2bWaaSbaaSqaaiaabMgacaaIYaaabeaakiaamYcacqGH % flY1cqGHflY1cqGHflY1caaJSaGaaeODamaaBaaaleaacaqGPbGaam % iraaqabaaakiaawIcacaGLPaaacaaJSaGaaWyAaiaam2dacaaJXaGa % aWilaiaamkdacaaJSaGaeyyXICTaeyyXICTaeyyXICTaaWilaiaam6 % eaaaa!632D! {V_{\rm{i}}} = \left( {{{\rm{v}}_{{\rm{i}}1}},{{\rm{v}}_{{\rm{i}}2}}, \cdot \cdot \cdot ,{{\rm{v}}_{{\rm{i}}D}}} \right),i = 1,2, \cdot \cdot \cdot ,N Vi=(vi1,vi2,,viD),i=1,2,,N
在第 t 代的第 i 个粒子向第 t+1 代进化时,根据如下式子更新:
V i j ( t + 1 ) = w v i j ( t ) + c 1 r 1 ( t ) [ p i j ( t ) − x i j ( t ) ] + c 2 r 2 ( t ) [ p g j ( t ) − x i j ( t ) ] % MathType!MTEF!2!1!+- % feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn % hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr % 4rNCHbWexLMBbXgBd9gzLbvyNv2CaeHbl7mZLdGeaGqiVu0Je9sqqr % pepC0xbbL8F4rqqrFfpeea0xe9Lq-Jc9vqaqpepm0xbba9pwe9Q8fs % 0-yqaqpepae9pg0FirpepeKkFr0xfr-xfr-xb9adbaqaaeGaciGaai % aabeqaamaabaabauaakeaacaWGwbWaaSbaaSqaaiaabMgacaqGQbaa % beaakiaacIcacaWG0bGaey4kaSIaaGymaiaacMcacqGH9aqpcaWG3b % GaamODamaaBaaaleaacaWGPbGaamOAaaqabaGccaGGOaGaamiDaiaa % cMcacqGHRaWkcaWGJbWaaSbaaSqaaiaaigdaaeqaaOGaamOCamaaBa % aaleaacaaIXaaabeaakiaacIcacaWG0bGaaiykaiaacUfacaWGWbWa % aSbaaSqaaiaadMgacaWGQbaabeaakiaacIcacaWG0bGaaiykaiabgk % HiTiaadIhadaWgaaWcbaGaamyAaiaadQgaaeqaaOGaaiikaiaadsha % caGGPaGaaiyxaiabgUcaRiaadogadaWgaaWcbaGaaGOmaaqabaGcca % WGYbWaaSbaaSqaaiaaikdaaeqaaOGaaiikaiaadshacaGGPaGaai4w % aiaadchadaWgaaWcbaGaam4zaiaadQgaaeqaaOGaaiikaiaadshaca % GGPaGaeyOeI0IaamiEamaaBaaaleaacaWGPbGaamOAaaqabaGccaGG % OaGaamiDaiaacMcacaGGDbaaaa!7618! {V_{{\rm{ij}}}}(t + 1) = w{v_{ij}}(t) + {c_1}{r_1}(t)[{p_{ij}}(t) - {x_{ij}}(t)] + {c_2}{r_2}(t)[{p_{gj}}(t) - {x_{ij}}(t)] Vij(t+1)=wvij(t)+c1r1(t)[pij(t)xij(t)]+c2r2(t)[pgj(t)xij(t)]

x i j ( t + 1 ) = x i j ( t ) + v i j ( t + 1 ) % MathType!MTEF!2!1!+- % feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn % hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr % 4rNCHbWexLMBbXgBd9gzLbvyNv2CaeHbl7mZLdGeaGqiVu0Je9sqqr % pepC0xbbL8F4rqqrFfpeea0xe9Lq-Jc9vqaqpepm0xbba9pwe9Q8fs % 0-yqaqpepae9pg0FirpepeKkFr0xfr-xfr-xb9adbaqaaeGaciGaai % aabeqaamaabaabauaakeaacaWG4bWaaSbaaSqaaiaadMgacaWGQbaa % beaakiaacIcacaWG0bGaey4kaSIaaGymaiaacMcacqGH9aqpcaWG4b % WaaSbaaSqaaiaadMgacaWGQbaabeaakiaacIcacaWG0bGaaiykaiab % gUcaRiaadAhadaWgaaWcbaGaamyAaiaadQgaaeqaaOGaaiikaiaads % hacqGHRaWkcaaIXaGaaiykaaaa!5467! {x_{ij}}(t + 1) = {x_{ij}}(t) + {v_{ij}}(t + 1) xij(t+1)=xij(t)+vij(t+1)

其中
速度更新公式的第一项 w v i j ( t ) wv_{ij}(t) wvij(t) 是惯性部分,由惯性权重和粒子自身速度构成,表示粒子对先前自身运动状态的信任; 第二项 c 1 r 1 ( t ) [ p i j ( t ) − x i j ( t ) ] c_1r_1(t)[p_{ij}(t)−x_{ij}(t)] c1r1(t)[pij(t)xij(t)] 是认知部分,表示粒子本身的思考,即粒子自己经验的部分,可理解为粒子当前位置与自身历史最优位置之间的距离和方向;第三项 c 2 r 2 ( t ) [ p g j ( t ) − x i j ( t ) ] c_2r_2(t)[p_{gj}(t)−x_{ij}(t)] c2r2(t)[pgj(t)xij(t)] 是社会部分,表示粒子之间的信息共享与合作,即来源于群体中其他优秀粒子的经验,可理解为粒子当前位置与群体历史最优位置之间的距离和方向。

注意
理论上,我们可以使用行(row)或者列(column)来表示向量(vector)。而在线性代数中,目前的符号表示惯例是用列(column)来表示向量(vector),我猜测原因在于简化计算公式。所以,如果只是单纯展示向量,向量可以采用行或列的表示方式。如果涉及到矩阵计算(多项式计算)则建议使用列(column)来表示向量(vector)。

假设某粒子当前位置C,个体极值位置B,全局最优位置A,那么该粒子下一步的运动状态如下图所示。

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-j3Xu3WyP-1662473600280)(pic2.png)]

V i j ( t + 1 ) = w v i j ( t ) + c 1 r 1 ( t ) [ p i j ( t ) − x i j ( t ) ] + c 2 r 2 ( t ) [ p g j ( t ) − x i j ( t ) ] % MathType!MTEF!2!1!+- % feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn % hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr % 4rNCHbWexLMBbXgBd9gzLbvyNv2CaeHbl7mZLdGeaGqiVu0Je9sqqr % pepC0xbbL8F4rqqrFfpeea0xe9Lq-Jc9vqaqpepm0xbba9pwe9Q8fs % 0-yqaqpepae9pg0FirpepeKkFr0xfr-xfr-xb9adbaqaaeGaciGaai % aabeqaamaabaabauaakeaacaWGwbWaaSbaaSqaaiaabMgacaqGQbaa % beaakiaacIcacaWG0bGaey4kaSIaaGymaiaacMcacqGH9aqpcaWG3b % GaamODamaaBaaaleaacaWGPbGaamOAaaqabaGccaGGOaGaamiDaiaa % cMcacqGHRaWkcaWGJbWaaSbaaSqaaiaaigdaaeqaaOGaamOCamaaBa % aaleaacaaIXaaabeaakiaacIcacaWG0bGaaiykaiaacUfacaWGWbWa % aSbaaSqaaiaadMgacaWGQbaabeaakiaacIcacaWG0bGaaiykaiabgk % HiTiaadIhadaWgaaWcbaGaamyAaiaadQgaaeqaaOGaaiikaiaadsha % caGGPaGaaiyxaiabgUcaRiaadogadaWgaaWcbaGaaGOmaaqabaGcca % WGYbWaaSbaaSqaaiaaikdaaeqaaOGaaiikaiaadshacaGGPaGaai4w % aiaadchadaWgaaWcbaGaam4zaiaadQgaaeqaaOGaaiikaiaadshaca % GGPaGaeyOeI0IaamiEamaaBaaaleaacaWGPbGaamOAaaqabaGccaGG % OaGaamiDaiaacMcacaGGDbaaaa!7618! {V_{{\rm{ij}}}}(t + 1) = w{v_{ij}}(t) + {c_1}{r_1}(t)[{p_{ij}}(t) - {x_{ij}}(t)] + {c_2}{r_2}(t)[{p_{gj}}(t) - {x_{ij}}(t)] Vij(t+1)=wvij(t)+c1r1(t)[pij(t)xij(t)]+c2r2(t)[pgj(t)xij(t)]
B 对应个体极值位置 p i p_i pi,A 对应群体极值位置 p g p_g pg,黄色向量为当前速度方向,绿色向量为向个体极值飞行步长,红色为向全局最值飞行步长。

2.2 算法流程


输入:参数 w w w: 0.5 ~ 0.8 , c 1 c_1 c1 c 2 c_2 c2:0.1 ~ 2(建议取值范围0.1 ~ 1) , r 1 ( t ) r_1(t) r1(t) , r 2 ( t ) r_2(t) r2(t) :区间[0,1]内的随机数,增加搜索的随机性,使粒子能充分搜索求解空间, vmax v m a x v_{max} vmax x m a x x_{max} xmax:取决于优化函数


Step1:初始化种群 x x x ;

Step2:计算个体 适应度 * ;

Step3:更新粒子速度 -> 更新粒子位置 * ;

Step4:并计算新位置的适应度,若新位置适应度更高,则将该粒子的位置进行更新,否则不更新。

Step5:判断是否满足终止条件 * ,是则退出,否则返回Step2。


输出:输出最优值。


注:

*1. 一般的,我们优化目标是 最小化一个函数值。所以个体计算出的函数值越小,适应度越高。如果优化目标是 最大化一个函数值 可以转化为 最小化一个函数值,因为最大值加负号是最小值, max f = min -f

*2. 注意更新速度后,先进行速度边界检测,一般采用 v ( v > v m a x ) = v m a x v(v>v_{max}) = v_{max} v(v>vmax)=vmax ,位置同理。

*3. 常见终止条件为设定迭代进化次数、适应度n代不再变化等。

3 算法分析

优点:

  • 原理比较简单,实现容易,参数少。

缺点:

  • 易早熟收敛至局部最优、迭代后期收敛速度慢的。

解释:

  • 标准粒子群算法的参数( w w w c 1 c_1 c1 c 2 c_2 c2)是固定的。惯性权重 w w w 描述的是粒子的“惯性”,在进化前期 w w w 应该大一些,保证各个粒子独立飞行充分搜索空间,后期应该小一点,多向其他粒子学习。自我学习因子 c 1 c_1 c1,群体学习因子 c 2 c_2 c2 分别为向个体极值和全局极值的飞行步长权重。前期 c 1 c_1 c1, 应该大一些,后期 c 2 c_2 c2 ,应该大一些,这样就能平衡粒子的全局搜索能力和局部搜索能力。3个参数共同影响了粒子的飞行方向,导致即使其他粒子找到更好的,但是当前粒子惯性太大,不能很快的飞向更优的位置。

4 算法拓展

针对标准PSO的缺点,通常有如下的改进:

  1. 实现参数的自适应变化。
  2. 引入一些其他机制,比如随机的因素,速度、位置的边界变化—后期压缩最大速度等。
  3. 结合其他智能优化算法:遗传算法、免疫算法、模拟退火算法等等,帮助粒子跳出局部最优,改善收敛速度。

5 案例实操

5.1 二维可行解

求解 𝑓 ( 𝑥 ) = 𝑥 𝑠 𝑖 𝑛 ( 𝑥 ) 𝑐 𝑜 𝑠 ( 2 𝑥 ) − 2 𝑥 𝑠 𝑖 𝑛 ( 3 𝑥 ) + 3 𝑥 𝑠 𝑖 𝑛 ( 4 𝑥 ) 𝑓(𝑥)=𝑥𝑠𝑖𝑛(𝑥)𝑐𝑜𝑠(2𝑥)−2𝑥𝑠𝑖𝑛(3𝑥)+3𝑥𝑠𝑖𝑛(4𝑥) f(x)=xsin(x)cos(2x)2xsin(3x)+3xsin(4x) 在区间 [0,50] 上的最小值

%% 计算函数最小值

% 转载请注明出处,谢谢~ 
% QQ : 1366420642

close all % 关闭所有Figure图窗
clear       % 从当前工作空间中删除所有变量,并系统内存中释放它们。

f= @(x) x .* sin(x) .* cos(2 * x) - 2 * x .* sin(3 * x) +3 * x .* sin(4 * x); % 函数表达式 % 求这个函数的最小值     

%% 初始化种群
N = 20;                         % 初始化种群内个体个数 
d = 1;                          % 可行解维数(即 X 独立参数的数目,X 的列数) 
ger = 100;                      % 最大迭代次数(求解中止条件)      
limit = [0, 50];                % 设置位置参数限制  
vlimit = [-10, 10];             % 设置速度限制  
w = 0.8;                        % 惯性权重(可通过单一变量原则迭代调整)  
c1 = 0.5;                       % 自我学习因子(可通过单一变量原则迭代调整)  
c2 = 0.5;                       % 群体学习因子(可通过单一变量原则迭代调整)  

figure(1);						% 新建画布1
ezplot(f,[0,0.01,limit(2)]);    % 绘制曲线(自变量范围为0 ~ limit(2)@50,绘图步长为0.01) % 不推荐使用 ezplot。请改用 fplot。有关详细信息,请参阅兼容性考虑。

x = limit(1) + (limit(2)-limit(1)) .* rand(N, d);% 采用随机分布初始化种群的位置,注意此处的 x 是 1 维的,x 的行数代表种群个体数,x 的列数代表粒子的维度 

v = rand(N, d);                  % 采用随机分布初始化种群的速度(N X d 阶取值范围(0,1)矩阵) 
xm = x;                          % 初始化个体的历史最佳适应度位置  
gm = zeros(1, d);                % 初始化种群的历史最佳适应度位置(1 X d 阶double类零矩阵) 
fxm = ones(N, 1)*inf;            % 初始化个体的历史最佳适应度(N X 1 阶正无穷inf矩阵)  
fgm = inf;                       % 初始化种群历史最佳适应度(正无穷inf)   
hold on                          % 继续在前一画布(画布1)上绘制曲线
plot(xm, f(gm), 'ro');           % 绘制粒子群初始散落位置
title('初始状态图');               % 画布1标题

%% 群体更新 
iter = 1;  						  % 初始化粒子群进化代数
% record = zeros(ger, 1);         % 粒子群进化代数记录器  
while iter <= ger  				  % 粒子群进化

     fx = f(x) ;                  % 计算粒子群所有个体当前适应度  
     
     for i = 1:N                  % 遍历粒子群
        if fx(i)  <fxm(i)         % 判断计算的适应度是否是个体历史最佳适应度
            fxm(i) = fx(i);       % 更新个体历史最佳适应度  
            xm(i,:) = x(i,:);     % 更新个体历史最佳适应度位置(取值)  
        end   
     end  
     
    if  min(fxm)  < fgm 
        [fgm, nmin] = min(fxm);   % 更新群体历史最佳适应度  
        gm = xm(nmin, :);         % 更新群体历史最佳位置  
    end  
    
    v = v * w + c1 * rand * (xm - x) + c2 * rand * (repmat(gm, N, 1) - x);% 速度更新 
    v(v > vlimit(2)) = vlimit(2); % 边界速度处理(逻辑索引、标量扩展) 
    v(v < vlimit(1)) = vlimit(1); % 边界速度处理(逻辑索引、标量扩展) 
    
    x = x + v;					  % 位置更新  
    x(x > limit(2)) = limit(2);   % 边界位置处理(逻辑索引、标量扩展)  
    x(x < limit(1)) = limit(1);   % 边界位置处理(逻辑索引、标量扩展)  
    
    record(iter) = fgm;           % 最小值记录  
    
    figure(2)                         % 查找 Number 属性等于 n 的图窗,并将其作为当前图窗。如果不存在具有该属性值的图窗,MATLAB® 将创建一个新图窗并将其 Number 属性设置为 n。
    x0 = 0 : 0.01 : limit(2);  
    subplot(1,2,1);
    plot(x0, f(x0), 'b-', x, f(x), 'ro'); % 每一次调用plot绘制图形,会将之前画布上对象清除。
    title('状态位置变化');
    subplot(1,2,2);
    plot(record);
    title('最优适应度进化过程')  
    pause(0.01);
    iter = iter+1;  

end  

figure(3);
x0 = 0 : 0.01 : limit(2);  
plot(x0, f(x0), 'b-', x, f(x), 'ro');
title('最终状态位置')  
disp(['最大值:',num2str(fgm)]);  % 功能上类似于printf
disp(['变量取值:',num2str(gm)]);  

5.2 三维可行解

求解 20 + 𝑥 2 + 𝑦 2 − 10 c o s ( 2 π x ) − 10 c o s ( 2 π 𝑦 ) 20+𝑥^2+𝑦^2−10cos(2πx)−10cos(2π𝑦) 20+x2+y210cos(2πx)10cos(2πy) 在区间 [−5.12,5.12] 上的最小值

% 转载请注明出处,谢谢~ 
% QQ : 1366420642

close all % 关闭所有Figure图窗
clear  % 从当前工作空间中删除所有变量,并系统内存中释放它们。
clc    % 清除命令窗口的内容。

f = @(x,y) 20 +  x.^2 + y.^2 - 10*cos(2*pi.*x)  - 10*cos(2*pi.*y) ;% 函数表达式
 
x0 = [-5.12:0.05:5.12];
y0 = x0 ;
[X,Y] = meshgrid(x0,y0); % 返回网格坐标
Z = f(X,Y);              % 通过网格坐标保证函数f计算了所有的离散化不重复的(x,y)组合
figure(1); 
mesh(X,Y,Z);             % 绘制网格曲面图,该函数将矩阵 `Z` 中的值绘制为由 `X` 和 `Y` 定义的 x-y 平面中的网格上方的高度。
colormap(parula(5));     % 曲面图中网格线边着色,边颜色因 `Z` 指定的高度而异。

%% 初始化种群 
N = 50;                         % 初始化种群内个体个数  
d = 2;                          % 可行解维数  
ger = 100;                      % 最大迭代次数       
limit = [-5.12,5.12];           % 设置位置参数限制  
vlimit = [-.5, .5];             % 设置速度限制  
w = 0.8;                        % 惯性权重  
c1 = 0.5;                       % 自我学习因子  
c2 = 0.5;                       % 群体学习因子   

x = limit(1) + (limit(2)-limit(1)) .* rand(N, d);% 采用随机分布初始化种群的位置,注意此处的 x 是 2 维的,x 的行数代表种群个体数,x 的列数代表粒子的维度

v = rand(N, d);                  % 初始化种群的速度,注意此处的 v 是 2 维的
xm = x;                          % 初始化每个个体的历史最佳适应度位置  
gm = zeros(1, d);                % 初始化种群的历史最佳适应度位置  
fxm = ones(N, 1)*inf;            % 初始化每个个体的历史最佳适应度   
fgm = inf;                       % 种群历史最佳适应度  

hold on; 
% [X,Y] = meshgrid(x(:,1),x(:,2)); % 种群的初始位置三维网格图
% Z = f( X,Y ) ;
scatter3( x(:,1),x(:,2) ,f( x(:,1),x(:,2) ),'r*' );% 种群的初始位置三维散点图
hold off;

%% 群体更新  
record=[];
iter = 1;  
% record = zeros(ger, 1);       % 记录器  
while iter <= ger  
     fx = f( x(:,1),x(:,2) ) ;  % 个体当前适应度     
     for i = 1:N        
        if  fx(i)  <fxm(i) 
            fxm(i) = fx(i);     % 更新个体历史最佳适应度  
            xm(i,:) = x(i,:);   % 更新个体历史最佳位置(取值)  
        end   
     end  
    if   min(fxm)<  fgm
        [fgm, nmin] = min(fxm); % 更新群体历史最佳适应度  
        gm = xm(nmin, :);       % 更新群体历史最佳位置  
    end  
    v = v * w + c1 * rand * (xm - x) + c2 * rand * (repmat(gm, N, 1) - x);% 速度更新  
    %% 边界速度处理(逻辑索引、标量扩展)  
    v(v > vlimit(2)) = vlimit(2);  
    v(v < vlimit(1)) = vlimit(1);  
    x = x + v;% 位置更新  
    %% 边界位置处理(逻辑索引、标量扩展)  
    x(x > limit(2)) = limit(2);  
    x(x < limit(1)) = limit(1);  
    
    record(iter) = fgm;%最大值记录  
    
    figure(2)                     % 查找 Number 属性等于 n 的图窗,并将其作为当前图窗。如果不存在具有该属性值的图窗,MATLAB® 将创建一个新图窗并将其 Number 属性设置为 n。
    subplot(1,2,1)
    mesh(X,Y,Z)
    hold on 
    scatter3( x(:,1),x(:,2) ,f( x(:,1),x(:,2) ) ,'r*');
    title(['状态位置变化','-迭代次数:',num2str(iter)])
    
    subplot(1,2,2);
    plot(record);
    title('最优适应度进化过程')  
    
    pause(0.01)  
    iter = iter+1; 

end  

figure(3);
mesh(X,Y,Z); 
hold on 
scatter3( x(:,1),x(:,2) ,f( x(:,1),x(:,2) ) ,'r*');
title('最终状态位置')  
disp(['最优值:',num2str(fgm)]);  
disp(['变量取值:',num2str(gm)]);

6 查漏补缺

6.1 MATLAB 控制指令

注意:内容参考于于MathWorks 帮助中心-…

MATLAB 控制指令一般写在 MATLAB 程序的第一行。

  • clc:清除命令窗口的内容。

    % `clc` 清除命令行窗口显示的所有输入和输出,给你一个“干净的屏幕”。使用 `clc` 后,您无法使用滚动条查看函数历史记录,但您仍然可以使用向上箭头键 ↑ 从命令历史记录中调用语句。
    clc
    
  • clear:从工作区中删除项目,释放系统内存。

    % 从当前工作空间中删除所有变量,从系统内存中释放它们。
    clear
    % 从内存中删除变量、脚本、函数或 MEX 函数 name1 ... nameN。
    	% 如果名字是一个函数,那么 clear 会重新初始化函数中的任何持久变量。
    	% 如果名字是当前正在执行的脚本或函数,那么 clear 不会将其删除
    	% 如果名字是一个被 mlock 锁定的函数,那么 clear 不会将其删除
    	% 如果名字是一个全局变量,那么 clear 会将其从当前工作区中删除,但仍保留在全局工作区中
    clear name1 ... nameN
    % 删除与列出的任何正则表达式匹配的所有变量,此选项仅删除变量。
    clear -regexp expr1 ... exprN
    % 删除由 ItemType 指示的项目类型,例如 all、functions、global、variables 、classes、import、java 或 mex。调用 clear all、clear classes 和 clear 函数会降低代码性能,并且通常是不必要的。
    	% 要从当前工作区中清除一个或多个特定变量,请使用 clear name1 ... nameN。
    	% 要清除当前工作区中的所有变量,请使用 clear 或 clearvars。
    	% 要清除所有全局变量,请使用 clear global 或 clearvars –global。
    	% 要清除特定类,请使用 clear myClass。
    	% 要清除特定函数或脚本,请使用 clear functionName。
    	% 要清除所有 MEX 函数,请使用 clear mex。
    	% clear 函数可以删除您指定的变量。要删除除少数指定变量之外的所有变量,请改用 clearvars。
    	% 如果清除图形或图形对象的句柄,则不会删除对象本身。使用 delete 删除对象。另一方面,删除对象不会删除用于存储其句柄的变量(如果有)。
    	% clear 函数不会清除 Simulink® 模型。请改用 bdclose。
    	% 在 UNIX® 系统上,clear 不会影响分配给 MATLAB 进程的内存量。
    clear ItemType
    
  • clf:清除当前图形窗口中的图形对象

    % 从当前图形中删除所有未隐藏句柄的图形对象(即,它们的 HandleVisibility 属性设置为 on)。
    clf
    % 从当前图形中删除所有图形对象,无论其 HandleVisibility 属性的设置如何,并将除 Position、Units、PaperPosition 和 PaperUnits 之外的所有图形属性重置为其默认值。
    clf('reset')
    % 清除带有句柄 fig 的单个图窗。
    clf(fig)
    % 清除带有句柄 fig 的单个图窗。
    clf(fig,'reset')
    % 返回图形的句柄。 这在图形 IntegerHandle 属性关闭时很有用,因为使用重置选项时非整数句柄变得无效(即 IntegerHandle 重置为 on,这是默认设置)。clf 命令在命令行上发出时的行为与在回调例程中的行为相同——它不识别回调的 HandleVisibility 设置。 这意味着当从回调例程中发出时,clf 只删除那些 HandleVisibility 属性设置为 on 的对象。
    figure_handle = clf(...)
    
  • close:删除当前图窗或指定图窗。 它可以选择返回关闭操作的状态。

    % 删除当前图窗(相当于 close(gcf))。
    close
    % 删除由 h 标识的图形。 如果 h 是一个数组,close 删除所有由 h 标识的图形。 h 也可以是数字。
    close(h)
    % 删除具有指定名称的图窗。
    close name
    % 删除所有未隐藏句柄的图窗。
    close all
    % 删除所有图窗,包括带有隐藏句柄的图窗。
    close all hidden
    % 删除所有图形,包括已将 CloseRequestFcn 更改为不关闭窗口的 GUI。
    close all force
    % 如果指定的窗口已被删除,则返回 1,否则返回 0。
    status = close(...)
    

6.2 MATLAB 匿名函数

注意:内容参考于MathWorks 帮助中心-匿名函数

6.2.1 什么是匿名函数?

匿名函数是 存储在程序文件中、但与数据类型是 function_handle 的变量相关的函数。匿名函数可以接受多个输入并返回一个输出。它们可能只包含一个可执行语句。

匿名函数是 MATLAB7.0 以后引入的新概念,可以用其来替代以往各种的函数,尤其是复杂函数的表达形式。 例如,对于一个复杂函数 y=f (x),倘若f (x)的表达式非常复杂,在后续代码中,调用这个函数的时候,只需要用y (1)即可代替当x=1时的那一串以x为自变量的表达式。 而且,这还为用隐函数来进行非线性拟合、提高隐函数或复杂函数求解效率,以及进行多重函数定义等提供了解决方案。

匿名函数可以在任意地方定义(包含命令行窗口),且由于没有固定的名称,函数是可以像变量一样被传递的。与之相对的即为按照.m文件存储的命名函数。

例如,创建用于计算平方数的匿名函数的句柄:

% 变量 sqr 是一个函数句柄。@ 运算符创建句柄,@ 运算符后面的圆括号 () 包括函数的输入参数。该匿名函数接受单个输入 x,并显式返回单个输出,即大小与包含平方值的 x 相同的数组。
% 注意,其实常数可以视为1维向量
sqr = @(x) x.^2;

通过将特定值 (5) 传递到函数句柄来计算该值的平方,与您将输入参数传递到标准函数一样。

a = sqr(5)
a =
   25

许多 MATLAB® 函数接受将函数句柄用作输入,这样您可以在特定值范围内计算函数。您可以为匿名函数或程序文件中的函数创建句柄。使用匿名函数的好处是不必为仅需要简短定义的函数编辑和维护文件。

例如,通过将函数句柄传递到 integral 函数,计算 sqr 函数从 01 范围内的积分:

q = integral(sqr,0,1);

您无需在工作区中创建变量以存储匿名函数。可以在表达式内创建临时函数句柄,例如这次对 integral 函数的调用:

q = integral(@(x) x.^2,0,1);
6.2.2 表达式中的变量

函数句柄不仅可以存储表达式,还能存储表达式进行计算需要的变量。

例如,为需要系数 abc 的匿名函数创建句柄。

a = 1.3;
b = .2;
c = 30;
parabola = @(x) a*x.^2 + b*x + c;

由于 abc 在您创建 parabola 时可用,该函数句柄包含这些值。即使您清除变量,这些值仍持久保留在函数句柄内:

clear a b c
x = 1;
y = parabola(x)
y =
   31.5000

要为这些系数提供不同值,您必须创建新的函数句柄:

a = -3.9;
b = 52;
c = 0;
parabola = @(x) a*x.^2 + b*x + c;

x = 1;
y = parabola(1)
y =
   48.1000

可以将函数句柄及其相关值存储在 MAT 文件中,然后使用 saveload 函数在后续的 MATLAB 会话中加载它们,例如

save myfile.mat parabola

在构造匿名函数时仅使用显式变量。如果匿名函数访问未在参数列表或主体中显式引用的任何变量或嵌套函数,则 MATLAB 会在您调用该函数时引发错误。隐式变量和函数调用通常会在 evalevalinassigninload 等函数中遇到。请避免在匿名函数主体中使用这些函数。

6.2.3 多个匿名函数

匿名函数中的表达式可以包含其他匿名函数。这可用于将不同的参数传递到在某一值范围内计算的函数。例如,您可以针对不同的
g ( c ) = ∫ 0 1 ( x 2 + c x + 1 ) d x % MathType!MTEF!2!1!+- % feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn % hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr % 4rNCHbWexLMBbXgBd9gzLbvyNv2CaeHbl7mZLdGeaGqiVu0Je9sqqr % pepC0xbbL8F4rqqrFfpeea0xe9Lq-Jc9vqaqpepm0xbba9pwe9Q8fs % 0-yqaqpepae9pg0FirpepeKkFr0xfr-xfr-xb9adbaqaaeGaciGaai % aabeqaamaabaabauaakeaacaqGNbGaaiikaiaadogacaGGPaGaeyyp % a0Zaa8qmaeaacaGGOaGaamiEamaaCaaaleqabaGaaGOmaaaakiabgU % caRiaadogacaWG4bGaey4kaSIaaGymaiaacMcaaSqaaiaaicdaaeaa % caaIXaaaniabgUIiYdGccaWGKbGaamiEaaaa!508D! {\rm{g}}(c) = \int_0^1 {({x^2} + cx + 1)} dx g(c)=01(x2+cx+1)dx
c 值求解以下方程,方法是合并使用两个匿名函数:

g = @(c) (integral(@(x) (x.^2 + c*x + 1),0,1));

下面介绍得出该语句的步骤:

  1. 将被积函数编写为匿名函数,

    @(x) (x.^2 + c*x + 1)
    
  2. 通过将函数句柄传递到 integral 在从 0 到 1 的范围内计算函数,

    integral(@(x) (x.^2 + c*x + 1),0,1)
    
  3. 通过为整个方程构造匿名函数以提供 c 的值,

    g = @(c) (integral(@(x) (x.^2 + c*x + 1),0,1));
    

最终的函数可以针对任何 c 值来求解方程。例如:

g(2)
ans =
   2.3333
6.2.4 不带输入的函数

如果您的函数不需要任何输入,请在定义和调用匿名函数时输入空的圆括号。例如:

t = @() datestr(now);
d = t()
d =
26-Jan-2012 15:11:47

在赋值语句中省略圆括号会创建另一函数句柄,并且不执行函数:

d = t
d = 
    @() datestr(now)
6.2.5 带有多个输入或输出的函数

匿名函数需要您像对标准函数一样显式指定输入参数,用逗号隔开多个输入。例如,以下函数接受两个输入 xy

myfunction = @(x,y) (x^2 + y^2 + x*y);

x = 1;
y = 10;
z = myfunction(x,y)
z = 111

然而,匿名函数只返回一个输出。如果函数中的表达式返回多个输出,您可以在调用该函数句柄时请求它们。

例如,ndgrid 函数可以返回与输入向量数量一样多的输出。调用 ndgrid 的匿名函数只返回一个输出 (mygrid)。通过调用 mygrid 可访问 ndgrid 函数返回的输出。

c = 10;
mygrid = @(x,y) ndgrid((-x:x/c:x),(-y:y/c:y));
[x,y] = mygrid(pi,2*pi);

您可以使用来自 mygrid 的输出创建网格图或曲面图:

z = sin(x) + cos(y);
mesh(x,y,z)
6.2.6 匿名函数的数组

虽然大多数的 MATLAB 基本数据类型支持多维数组,但函数句柄必须是标量(单个元素)。但您可以使用元胞数组或结构体数组存储多个函数句柄。最常见的方式是使用元胞数组,例如

f = {@(x)x.^2;
     @(y)y+10;
     @(x,y)x.^2+y+10};

创建元胞数组时,记住 MATLAB 将空格解释为列分隔符。如上面的代码所示,省略表达式中的空格,或将表达式括在圆括号中,例如

f = {@(x) (x.^2);
     @(y) (y + 10);
     @(x,y) (x.^2 + y + 10)};

使用花括号访问元胞内容。例如,f{1} 返回第一个函数句柄。要执行该函数,请在花括号之后的圆括号中传递输入值:

x = 1;
y = 10;

f{1}(x)
f{2}(y)
f{3}(x,y)
ans =
     1

ans =
    20

ans =
    21

6.3 MATLAB 数组、向量、矩阵

注意:内容参考于 MathWorks 帮助中心-矩阵和数组

MATLAB 是“matrix laboratory”的缩写形式。MATLAB® 主要用于处理整个的矩阵和数组,而其他编程语言大多逐个处理数值。

所有 MATLAB 变量都是多维数组,与数据类型无关。矩阵 是指通常用来进行线性代数运算的二维数组。

6.3.1 数组创建

要创建每行包含四个元素的数组,请使用逗号 (,) 或空格分隔各元素。

a = [1 2 3 4]
a = 1×4

1     2     3     4

这种数组为行向量。

要创建每列包含四个元素的数组,请使用分号 ( ; ) 分隔各元素。

a = [1 ; 2 ; 3 ; 4]

a =

     1
     2
     3
     4

这种数组为列向量

线性代数中,目前的符号表示惯例是用列(column)来表示向量(vector)。向量维数是 ,因为向量的坐标只有一行,列数表示它的维数。普遍使用列向量,是为了方便进行矩阵与向量的乘法运算。

要创建包含多行的矩阵,请使用分号分隔各行。

a = [1 3 5; 2 4 6; 7 8 10]
a = 3×3
 
 1     3     5
 2     4     6
 7     8    10

创建矩阵的另一种方法是使用 ones、zeros 或 rand 等函数。例如,创建一个由零组成的 5×1 列向量。

矩阵的本质就是线性方程式,两者是一一对应关系

z = zeros(5,1)
z = 5×1

 0
 0
 0
 0
 0
6.3.2 矩阵和数组运算

MATLAB 允许您使用单一的算术运算符或函数来处理矩阵中的所有值。

a + 10
ans = 3×3

11    13    15
12    14    16
17    18    20
sin(a)
ans = 3×3

0.8415    0.1411   -0.9589
0.9093   -0.7568   -0.2794
0.6570    0.9894   -0.5440

要转置矩阵,请使用单引号 ('):

 a'
ans = 3×3
 
 1     2     7
 3     4     8
 5     6    10

您可以使用 * 运算符执行标准矩阵乘法,这将计算行与列之间的内积。例如,确认矩阵乘以其逆矩阵可返回单位矩阵:

p = a*inv(a)
p = 3×3

1.0000         0         0
     0    1.0000         0
     0   -0.0000    1.0000

请注意,p 不是整数值矩阵。MATLAB 将数字存储为浮点值,算术运算可以区分实际值与其浮点表示之间的细微差别。使用 format 命令可以显示更多小数位数:

format long
p = a*inv(a)
p = 3×3

1.0000         0         0
     0    1.0000         0
     0   -0.0000    1.0000

使用以下命令将显示内容重置为更短格式

format short

format 仅影响数字显示,而不影响 MATLAB 对数字的计算或保存方式。

要执行元素级乘法(而非矩阵乘法),请使用 .* 运算符:

p = a.*a
p = 3×3

 1     9    25
 4    16    36
49    64   100

乘法、除法和幂的矩阵运算符分别具有执行元素级运算的对应数组运算符。例如,计算 a 的各个元素的三次方:

a.^3
ans = 3×3

  1          27         125
  8          64         216
343         512        1000
6.3.3 串联

串联是连接数组以便形成更大数组的过程。实际上,第一个数组是通过将其各个元素串联起来而构成的。成对的方括号 [] 即为串联运算符。

A = [a,a]
A = 3×6

 1     3     5     1     3     5
 2     4     6     2     4     6
 7     8    10     7     8    10

使用逗号将彼此相邻的数组串联起来称为水平串联。每个数组必须具有相同的行数。同样,如果各数组具有相同的列数,则可以使用分号垂直串联。

A = [a; a]
A = 6×3

 1     3     5
 2     4     6
 7     8    10
 1     3     5
 2     4     6
 7     8    10
6.3.4 复数

复数包含实部和虚部,虚数单位是 -1 的平方根。

sqrt(-1)
ans = 0.0000 + 1.0000i

要表示复数的虚部,请使用 i 或 j。

c = [3+4i, 4+3j; -i, 10j]
c = 2×2 complex

   3.0000 + 4.0000i   4.0000 + 3.0000i
   0.0000 - 1.0000i   0.0000 +10.0000i

6.4 MATLAB 索引

注意:内容参考于MathWorks 帮助中心-索引

6.4.1 下标

A 的行 i 和列 j 中的元素通过 A(i,j) 表示。例如,A(4,2) 表示第四行和第二列中的数字。在幻方矩阵中,A(4,2) 为 15。

>> A=magic(4)

A =

    16     2     3    13
     5    11    10     8
     9     7     6    12
     4    14    15     1

因此,要计算 A 第四列中的元素的总和,请键入

A(1,4) + A(2,4) + A(3,4) + A(4,4)

此下标生成

ans =
     34

但这不是计算某列总和的最佳方法。

此外,还可以使用单一下标 A(k) 引用矩阵的元素。单一下标是引用行向量(N X 1 阶)或列向量(1 X N 阶)的常见方法。但是,也可以对满二维矩阵(2 X 2 阶)应用单一下标。在这种情况下,索引数组被视为一个由原始矩阵的所有列从左至右依次取出构成的长列向量。因此,在幻方矩阵中,A(12) 是另一种引用存储在 A(4,2) 中的值 15 的方法。

A( 1)=A(1,1)=16
A( 2)=A(2,1)=5
A( 3)=A(3,1)=9
A( 4)=A(4,1)=4
A( 5)=A(1,2)=2
A( 6)=A(2,2)=11
A( 7)=A(3,2)=7
A( 8)=A(4,2)=14
A( 9)=A(1,3)=3
A(10)=A(2,3)=10
A(11)=A(3,3)=6
A(12)=A(4,3)=15
A(13)=A(1,4)=13
A(14)=A(2,4)=8
A(15)=A(3,4)=12
A(16)=A(4,4)=1

如果尝试使用矩阵外部元素的值,则会生成错误:

t = A(4,5)

索引超出矩阵维度。
相反,如果将值存储在矩阵外部元素中,则会增大大小以便容纳新元素:

X = A;
X(4,5) = 17

X =
    16     3     2    13     0
     5    10    11     8     0
     9     6     7    12     0
     4    15    14     1    17
6.4.2 冒号运算符

冒号 : 是最重要的 MATLAB® 运算符之一。它以多种不同形式出现。表达式

1:10

是包含从 1 到 10 之间的整数的行向量:

1     2     3     4     5     6     7     8     9    10

要获取非单位间距,请指定增量。例如,

100:-7:50

100    93    86    79    72    65    58    51

0:pi/4:pi

0    0.7854    1.5708    2.3562    3.1416

包含冒号的下标表达式引用部分矩阵:

A(1:k,j)

表示 A 第 j 列中的前 k 个元素。因此,

sum(A(1:4,4))

计算第四列的总和。但是,执行此计算有一种更好的方法。冒号本身引用矩阵行或列中的所有元素,而关键字 end 引用最后一个行或列。因此,

sum(A(:,end))

计算 A 最后一列中的元素的总和:

ans =
     34

为什么 4×4 幻方矩阵的幻数和等于 34?如果将介于 1 到 16 之间的整数分为四个总和相等的组,该总和必须为

sum(1:16)/4

当然,也即

ans =
     34
6.4.3 串联

串联是连接小矩阵以便形成更大矩阵的过程。实际上,第一个矩阵是通过将其各个元素串联起来而构成的。成对的方括号 [] 即为串联运算符。例如,从 4×4 幻方矩阵 A 开始,组成

B = [A  A+32; A+48  A+16]

结果会生成一个 8×8 矩阵,这是通过连接四个子矩阵获得的:

B =

16     3     2    13    48    35    34    45
 5    10    11     8    37    42    43    40
 9     6     7    12    41    38    39    44
 4    15    14     1    36    47    46    33
64    51    50    61    32    19    18    29
53    58    59    56    21    26    27    24
57    54    55    60    25    22    23    28
52    63    62    49    20    31    30    17

此矩阵是一个接近于幻方矩阵的矩阵。此矩阵的元素是经过重新排列的整数 1:64。此矩阵的列总和符合 8×8 幻方矩阵的要求:

sum(B)

ans =
   260   260   260   260   260   260   260   260

但是其行总和 sum(B’)’ 并不完全相同。要使其成为有效的 8×8 幻方矩阵,需要进行进一步操作。

6.4.4 删除行和列

只需使用一对方括号即可从矩阵中删除行和列。首先

X = A;

然后,要删除 X 的第二列,请使用

X(:,2) = []

这会将 X 更改为

X =
    16     2    13
     5    11     8 
     9     7    12
     4    14     1

如果您删除矩阵中的单个元素,结果将不再是矩阵。因此,以下类似表达式

X(1,2) = []

将会导致错误。但是,使用单一下标可以删除一个元素或元素序列,并将其余元素重构为一个行向量。因此

X(2:2:10) = []

生成

X =
    16     9     2     7    13    12     1
6.4.5 标量扩展

可以采用多种不同方法将矩阵和标量合并在一起。例如,通过从每个元素中减去标量而将其从矩阵中减去。幻方矩阵的元素平均值为 8.5,因此

B = A - 8.5

形成一个列总和为零的矩阵:

B =
      7.5     -5.5     -6.5      4.5
     -3.5      1.5      2.5     -0.5
      0.5     -2.5     -1.5      3.5
     -4.5      6.5      5.5     -7.5

sum(B)

ans =
     0     0     0     0

通过标量扩展,MATLAB 会为范围中的所有索引分配一个指定标量。例如,

B(1:2,2:3) = 0

将 B 的某个部分清零:

B =
      7.5        0        0      4.5
     -3.5        0        0     -0.5
      0.5     -2.5     -1.5      3.5
     -4.5      6.5      5.5     -7.5
6.4.6 逻辑下标

根据逻辑和关系运算创建的逻辑向量可用于引用子数组。

假定 X 是一个普通矩阵,L 是一个由某个逻辑运算生成的与 X 同等大小的逻辑型元素的矩阵。那么 X(L) 指定 X 的元素,其中 L 的元素非逻辑零。注意:X(L) 的本质是一个由 X 所有列向量索引依次串联形成的一个新的索引列向量,L 的逻辑零元素会屏蔽对应位置 X 的元素索引,即该位置的元素索引不会出现在新的索引列向量中。

>> A=magic(4)

A =

    16     2     3    13
     5    11    10     8
     9     7     6    12
     4    14    15     1

>> L=A>8

L =

  4×4 logical array

   1   0   0   1
   0   1   1   0
   1   0   0   1
   0   1   1   0

% A(L)其实是A的元素索引组成的列向量,直接输出显示的是对应元素值
>> A(L)

ans =

    16
     9
    11
    14
    10
    15
    13
    12

>> B=A(L)

B =

    16
     9
    11
    14
    10
    15
    13
    12

% A(L)其实是A的元素索引组成的列向量,在该赋值语句中可以体现出来
>> A(L)=0

A =

     0     2     3     0
     5     0     0     8
     0     7     6     0
     4     0     0     1

通过将逻辑运算指定为下标表达式,可以在一个步骤中完成这种下标。假定您具有以下数据集:

x = [2.1 1.7 1.6 1.5 NaN 1.9 1.8 1.5 5.1 1.8 1.4 2.2 1.6 1.8];

NaN 是用于缺少的观测值的标记,例如,无法响应问卷中的某个项。要使用逻辑索引删除缺少的数据,请使用 isfinite(x),对于所有有限数值,该函数为 true;对于 NaN 和 Inf,该函数为 false:

x = x(isfinite(x))
x =
  2.1 1.7 1.6 1.5 1.9 1.8 1.5 5.1 1.8 1.4 2.2 1.6 1.8

现在,存在一个似乎与其他项很不一样的观测值,即 5.1。这是一个离群值。下面的语句可删除离群值,在本示例中,即比均值大三倍标准差的元素:

x = x(abs(x-mean(x)) <= 3*std(x))
x =
  2.1 1.7 1.6 1.5 1.9 1.8 1.5 1.8 1.4 2.2 1.6 1.8

标量扩展的另一示例:请使用逻辑索引和标量扩展将非质数设置为 0,以便高亮显示丢勒幻方矩阵中的质数的位置。(请参阅 magic 函数。)

A(~isprime(A)) = 0

A =
     0     3     2    13
     5     0    11     0
     0     0     7     0
     0     0     0     0
6.4.7 find 函数

find 函数可用于确定与指定逻辑条件相符的数组元素的索引。find 以最简单的形式返回索引的列向量。转置该向量以便获取索引的行向量。例如,再次从丢勒的幻方矩阵开始。(请参阅 magic 函数。)

k = find(isprime(A))'

使用一维索引选取幻方矩阵中的质数的位置:

k =
     2     5     9    10    11    13

使用以下命令按 k 确定的顺序将这些质数显示为行向量

A(k)

ans =
     5     3     2    11     7    13

将 k 用作赋值语句的左侧索引时,会保留矩阵结构:

A(k) = NaN

A =
    16   NaN   NaN   NaN
   NaN    10   NaN     8
     9     6   NaN    12
     4    15    14     1

6.5 MATLAB 二维和三维网格坐标

注意:内容参考于MathWorks 帮助中心-二维和三维网格

二维和三维网格坐标

meshgrid
6.5.1 语法
% 基于向量 x 和 y 中包含的坐标返回二维网格坐标。X 是一个矩阵,每一行是 x 的一个副本;Y 也是一个矩阵,每一列是 y 的一个副本。坐标 X 和 Y 表示的网格有 length(y) 个行和 length(x) 个列。
[X,Y] = meshgrid(x,y)

% 与 [X,Y] = meshgrid(x,x) 相同,并返回网格大小为 length(x)×length(x) 的方形网格坐标。
[X,Y] = meshgrid(x)

% 返回由向量 x、y 和 z 定义的三维网格坐标。X、Y 和 Z 表示的网格的大小为length(y) × length(x) × length(z)。
[X,Y,Z] = meshgrid(x,y,z)

% 与 [X,Y,Z] = meshgrid(x,x,x) 相同,并返回网格大小为 length(x)×length(x)×length(x) 的三维网格坐标。
[X,Y,Z] = meshgrid(x)
6.5.2 二维网格坐标

使用向量 x 定义的 x 坐标和向量 y 定义的 y 坐标创建二维网格坐标。二维网格中的点可以表示为: ( x i j , y i j ) (x_{ij},y_{ij}) xijyij

x = 1:3;
y = 1:5;
[X,Y] = meshgrid(x,y)

%% 对 x,y 网格化生成的 X,Y 矩阵形式,保证了对应的多项式取出来,所有离散的不重复的(x,y)元素组合
X = 5×3
 1     2     3
 1     2     3
 1     2     3
 1     2     3
 1     2     3
 
Y = 5×3
 1     1     1
 2     2     2
 3     3     3
 4     4     4
 5     5     5
 
 
[X,Y]

ans =

     1     2     3     1     1     1
     1     2     3     2     2     2
     1     2     3     3     3     3
     1     2     3     4     4     4
     1     2     3     5     5     5

在二维网格上计算表达式 x 2 + y 2 x^2+y^2 x2+y2

X.^2 + Y.^2
ans = 5×3

 2     5    10
 5     8    13
10    13    18
17    20    25
26    29    34
6.5.3 绘制曲面图

使用均匀分布的 x 坐标和 y 坐标在区间 [-2,2] 内创建二维网格。

x = -2:0.25:2;
y = x;
[X,Y] = meshgrid(x);

在二维网格上计算并绘制函数 f ( x , y ) = x e − x 2 − y 2 f(x,y)=xe^{−x^2−y^2} f(x,y)=xex2y2

F = X.*exp(-X.^2-Y.^2);
surf(X,Y,F)

从 R2016b 开始,操作网格之前并不总是需要先创建网格。例如,计算表达式 x e − x 2 − y 2 xe^{−x^2−y^2} xex2y2,将隐式扩展向量 x 和 y。有关隐式扩展的详细信息,请参阅数组与矩阵运算。

surf(x,y,x.*exp(-x.^2-(y').^2))
6.5.4 三维网格坐标

在区间 [0,6] 内使用定义的 x、y 和 z 坐标创建三维网格坐标,并计算表达式 x 2 + y 2 + z 2 x^2+y^2+z^2 x2+y2+z2

x = 0:2:6;
y = 0:1:6;
z = 0:3:6;
[X,Y,Z] = meshgrid(x,y,z);
F = X.^2 + Y.^2 + Z.^2;

确定网格的大小。三个坐标向量具有不同的长度,构成一个网格点矩形框。

gridsize = size(F)

gridsize = 1×3

7     4     3

使用单输入语法,基于 x 中定义的坐标生成均匀分布的三维网格。新网格构成一个网格点正方体。

[X,Y,Z] = meshgrid(x);
G = X.^2 + Y.^2 + Z.^2;
gridsize = size(G)

gridsize = 1×3
4     4     4

6.6 MATLAB 网格曲面图

注意:内容参考于MathWorks 帮助中心-网格曲面图

网格曲面图

mesh
6.6.1 语法
% 创建一个网格图,该网格图为三维曲面,有实色边颜色,无面颜色。该函数将矩阵 `Z` 中的值绘制为由 `X` 和 `Y` 定义的 x-y 平面中的网格上方的高度。边颜色因 `Z` 指定的高度而异。
mesh(X,Y,Z)
% 创建一个网格图,并将 `Z` 中元素的列索引和行索引用作 x 坐标和 y 坐标。
mesh(Z)
%  进一步指定边的颜色。
mesh(Z,C)
% 进一步指定边的颜色。
mesh(___,C)
% 将图形绘制到 `ax` 指定的坐标区中,而不是当前坐标区中。指定坐标区作为第一个输入参数。
mesh(ax,___)
% 使用一个或多个名称-值对组参数指定曲面属性。例如,`'FaceAlpha',0.5` 创建半透明网格图。
mesh(___,Name,Value)
% 将返回一个图曲面对象。在创建网格图后,使用 `s` 修改网格图。有关属性列表,请参阅 [Surface 属性](https://ww2.mathworks.cn/help/matlab/ref/matlab.graphics.chart.primitive.surface-properties.html)。
s = mesh(___)
6.6.2 创建网格图

创建三个相同大小的矩阵。然后将它们绘制为一个网格图。该绘图使用 Z 确定高度和颜色。

[X,Y] = meshgrid(-8:.5:8);
R = sqrt(X.^2 + Y.^2) + eps;
Z = sin(R)./R;
mesh(X,Y,Z)
6.6.3 为网格图指定颜色图颜色

通过包含第四个矩阵输入 C 来指定网格图的颜色。网格图使用 Z 确定高度,使用 C 确定颜色。使用颜色图指定颜色,该颜色图使用单个数字表示色谱上的颜色。使用颜色图时,CZ 大小相同。向图中添加颜色栏以显示 C 中的数据值如何对应于颜色图中的颜色。

[X,Y] = meshgrid(-8:.5:8);
R = sqrt(X.^2 + Y.^2) + eps;
Z = sin(R)./R;
C = X.*Y;
mesh(X,Y,Z,C)
colorbar
6.6.4 为网格图指定真彩色

通过包含第四个矩阵输入 CO 来指定网格图的颜色。网格图使用 Z 确定高度,使用 CO 确定颜色。使用真彩色指定颜色,真彩色使用三个数字(即三元组)表示所有可能的颜色。使用真彩色时,如果 Zm×n,则 COm×n×3。数组的第一页指示每种颜色的红色分量;第二页指示绿色分量;第三页指示蓝色分量。

[X,Y,Z] = peaks(25);
CO(:,:,1) = zeros(25); % red
CO(:,:,2) = ones(25).*linspace(0.5,0.6,25); % green
CO(:,:,3) = ones(25).*linspace(0,1,25); % blue
mesh(X,Y,Z,CO)
6.6.5 修改网格图外观

指定值为 0.5FaceAlpha 名称-值对组,以创建半透明网格曲面。要允许进一步修改,请将曲面对象赋给变量 s

[X,Y] = meshgrid(-5:.5:5);
Z = Y.*sin(X) - X.*cos(Y);
s = mesh(X,Y,Z,'FaceAlpha','0.5')
s = 
  Surface with properties:

       EdgeColor: 'flat'
       LineStyle: '-'
       FaceColor: [1 1 1]
    FaceLighting: 'none'
       FaceAlpha: 0.5000
           XData: [21x21 double]
           YData: [21x21 double]
           ZData: [21x21 double]
           CData: [21x21 double]

  Show all properties

使用 s 可在创建曲面后访问和修改网格图的属性。例如,通过设置 FaceColor 属性,为网格图的面添加颜色。

s.FaceColor = 'flat';

6.7 MATLAB 颜色图

注意:内容参考于MathWorks 帮助中心-颜色图

查看并设置当前颜色图

colormap
6.7.1 语法
% 将当前图窗的颜色图设置为预定义的颜色图之一。如果您为图窗设置了颜色图,图窗中的坐标区和图将使用相同的颜色图。新颜色图的长度(颜色数)与当前颜色图相同。当您使用此语法时,不能为颜色图指定自定义长度。有关颜色图的详细信息,请参阅什么是颜色图?。
colormap map
% 将当前图窗的颜色图设置为 map 指定的颜色图。
colormap(map)
% 为 target 指定的图窗、坐标区或图形设置颜色图,而不是为当前图窗设置颜色图。
colormap(target,map)
% 返回当前图窗的颜色图,形式为 RGB 三元组组成的三列矩阵。
cmap = colormap
% 返回 target 指定的图窗、坐标区或图的颜色图。
cmap = colormap(target)
6.7.2 更改图窗的颜色图
surf(peaks)
colormap winter  % 创建一个曲面图并将颜色图设置为 winter
colormap summer  % 将当前图窗的颜色图更改为 summer
colormap default % 将颜色图设置回您系统的默认值。如果您尚未指定不同默认值,则默认颜色图是 parula。
6.7.3 对图窗中的每个坐标区使用不同的颜色图

从 R2019b 开始,您可以使用 tiledlayoutnexttile 函数显示分块图。调用 tiledlayout 函数以创建一个 2×1 分块图布局。调用 nexttile 函数以创建坐标区对象 ax1ax2。通过将坐标区对象传递给 colormap 函数,为每个坐标区指定不同的颜色图。在上坐标区中,使用 spring 颜色图创建一个曲面图。在下坐标区中,使用 winter 颜色图创建一个曲面图。

tiledlayout(2,1)
ax1 = nexttile;
surf(peaks)
colormap(ax1,spring)

ax2 = nexttile; 
surf(peaks)
colormap(ax2,winter)
6.7.4 指定颜色图的颜色数

通过将整数作为内置颜色图的输入参数传递来指定用于颜色图的颜色数。使用 parula 颜色图中的五种颜色。

mesh(peaks)
colormap(parula(5))
6.7.5 创建自定义颜色图

通过定义一个由介于 0.0 和 1.0 之间的值组成的三列矩阵来创建一个自定义颜色图。每行定义一个三元素 RGB 三元组。第一列指定红色强度。第二列指定绿色强度。第三列指定蓝色强度。

通过将前两个列设置为零来使用蓝色值的颜色图。

map = [0 0 0.3
    0 0 0.4
    0 0 0.5
    0 0 0.6
    0 0 0.8
    0 0 1.0];

surf(peaks)
colormap(map)
6.7.6 返回用在绘图中的颜色图值

创建 peaks 函数的曲面图并指定颜色图。

mesh(peaks)
colormap(autumn(5))

返回定义用在绘图中使用的颜色的值的三列矩阵。每行是一个指定颜色图的一种颜色的 RGB 三元组颜色值。

cmap = colormap
cmap = 5×3

    1.0000         0         0
    1.0000    0.2500         0
    1.0000    0.5000         0
    1.0000    0.7500         0
    1.0000    1.0000         0
6.7.7 返回特定坐标区的颜色图值

通过将坐标区对象传递给 colormap 函数,返回特定坐标区的颜色图值。

使用 tiledlayoutnexttile 函数创建两个分块图,这两个函数是从 R2019b 开始推出的新函数。调用 tiledlayout 函数以创建一个 2×1 分块图布局。调用 nexttile 函数以创建坐标区对象 ax1ax2。然后显示两个以不同颜色图填充的等高线图。

tiledlayout(2,1)
ax1 = nexttile;
contourf(peaks)
colormap(ax1,hot(8))

ax2 = nexttile;
contourf(peaks)
colormap(ax2,pink)

通过将 ax1 传递给 colormap 函数,返回上部绘图中使用的颜色图值。每行是一个指定颜色图的一种颜色的 RGB 三元组颜色值。

cmap = colormap(ax1)
cmap = 8×3

    0.3333         0         0
    0.6667         0         0
    1.0000         0         0
    1.0000    0.3333         0
    1.0000    0.6667         0
    1.0000    1.0000         0
    1.0000    1.0000    0.5000
    1.0000    1.0000    1.0000
6.7.8 将图窗的颜色图更改为图像

加载 spine 数据集以返回 X 及其关联的颜色图 map。使用 image 函数显示 X 并将颜色图设置为 map

load spine
image(X)
colormap(map)

6.8 MATLAB 基本绘图函数

注意:内容参考于MathWorks 帮助中心-基本绘图函数

6.8.1 创建绘图

plot 函数具有不同的形式,具体取决于输入参数。

  • 如果 y 是向量,plot(y) 会生成 y 元素与 y 元素索引的分段线图。
  • 如果有两个向量被指定为参数,plot(x,y) 会生成 yx 的图形。

使用冒号运算符创建从 0 至 2πx 值向量,计算这些值的正弦,并绘制结果。

x = 0:pi/100:2*pi;
y = sin(x);
plot(x,y)

添加轴标签和标题。xlabel 函数中的字符 \pi 用于创建符号 πtitle 函数中的 FontSize 属性用于增大标题所用的文本大小。

xlabel('x = 0:2\pi')
ylabel('Sine of x')
title('Plot of the Sine Function','FontSize',12)
6.8.2 在一幅图形中绘制多个数据集

通过调用一次 plot,多个 x-y 对组参数会创建多幅图形。MATLAB® 对每条线使用不同的颜色。

例如,下列语句绘制 x 的三个相关函数:

x = 0:pi/100:2*pi;
y = sin(x);
y2 = sin(x-.25);
y3 = sin(x-.5);
plot(x,y,x,y2,x,y3)

legend 函数提供了一种标识各条线的简单方法:

legend('sin(x)','sin(x-.25)','sin(x-.5)')
6.8.3 指定线型和颜色

使用 plot 命令绘制数据时,可以指定颜色、线型和标记(例如加号或圆圈):

plot(x,y,'color_style_marker')

color_style_marker 包含一至四个字符(包括在单引号中),这些字符根据颜色、线型和标记类型构造而成。例如,

plot(x,y,'r:+')

使用红色点线绘制数据,并在每个数据点处放置一个 + 标记。

color_style_marker 由下列元素的组合形式构成。

类型含义
颜色'c' 'm' 'y' 'r' 'g' 'b' 'w' 'k'青蓝 品红 黄 红 绿 蓝 白 黑
线型'-' '--' ':' '-.' 无字符实线 虚线 点线 点划线 没有线条
标记类型'+' 'o' '*' 'x' 's' 'd' '^' 'v' '>' '<' 'p' 'h' 无字符加号 空心圆 星号 字母 x 空心正方形 空心菱形 空心上三角 空心下三角 空心右三角 空心左三角 空心五角形 空心六角形 无标记
6.8.5 绘制线条和标记

如果指定标记类型,但未指定线型,MATLAB® 仅使用标记创建图形,而不会创建线条。例如,

plot(x,y,'ks')

在每个数据点绘制黑色正方形,但不会使用线条连接标记。语句

plot(x,y,'r:+')v

绘制红色点线,并在每个数据点处放置加号标记。

6.8.6 在每十个数据点处放置标记

此示例展示如何使用比绘制线条所用的数据点更少的数据点来绘制标记。它使用点线图和标记图(分别采用不同数目的数据点)绘制两次数据图:

x1 = 0:pi/100:2*pi;
x2 = 0:pi/10:2*pi;
plot(x1,sin(x1),'r:',x2,sin(x2),'r+')
6.8.7 绘制虚数和复数数据

将多个复数值作为参数传递给 plot 时,MATLAB 会忽略虚部,但传递一个复数参数时除外。对于这一特殊情况,该命令是绘制虚部对实部的图的一种快捷方式。因此,

plot(Z)

其中 Z 是复数向量或矩阵,等效于

plot(real(Z),imag(Z))

下列语句将绘制一个具有 20 条边的多边形,并在各顶点处绘制一个小圆圈。

t = 0:pi/10:2*pi;
plot(exp(1i*t),'-o')
axis equal

axis equal 命令使 xy 轴上的各刻度线增量的长度相同,这会使此绘图看起来更加圆润。

6.8.8 将绘图添加到现有图形中

hold 命令用于将绘图添加到现有图形中。当键入

hold on

时,MATLAB 不会在您发出其他绘图命令时替换现有图形。MATLAB 而会将新图形与当前图形合并在一起。

例如,下列语句首先创建 peaks 函数的曲面图,然后叠加同一函数的等高线图:

[x,y,z] = peaks;
% Create surface plot
surf(x,y,z)
% Remove edge lines a smooth colors
shading interp
% Hold the current graph 
hold on
% Add the contour graph to the pcolor graph
contour3(x,y,z,20,'k')
% Return to default
hold off
6.8.9 图窗窗口

如果尚未创建图窗窗口,绘图函数会自动打开一个新的图窗窗口。如果打开了多个图窗窗口,MATLAB 将使用指定为“当前图窗”(通常为上次使用的图窗)的图窗窗口。即 MATLAB 绘制函数图形使用的是获得焦点的“当前图窗”。

注意:如果绘制图形比较复杂,在未绘制完成之前,使用鼠标点击了其他图窗,会导致 MATLAB 接下来的绘制工作在刚点击的图窗中进行,与预期不一致。

要将现有图窗窗口设置为当前的图窗,请将指针放置在该窗口中并点击鼠标,或者也可以键入

figure(n)

其中 n 是图窗标题栏中的编号。

要打开新的图窗窗口并将其作为当前图窗,请键入

figure
6.8.10 清空图窗以便创建新绘图

如果某图窗已存在,大多数绘图命令会清除轴并使用此图窗创建新绘图。但是,这些命令不会重置图窗属性,例如,背景色或颜色图。如果已在以前的绘图中设置图窗属性,您可以先使用带有 reset 选项的 clf 命令。

clf reset

然后创建新绘图,以便将此图窗的属性恢复为其默认值。

6.8.11 在一幅图窗中显示多个绘图

subplot 命令用于在同一窗口中显示多个绘图,或者在同一张纸上打印这些绘图。键入以下命令

subplot(m,n,p)

会将图窗窗口划分为由多个小子图组成的 m×n 矩阵,并选择第 p 个子图作为当前绘图。这些绘图沿图窗窗口的第一行进行编号,然后沿第二行进行编号,依此类推。例如,下列语句在图窗窗口的三个子区域中绘制数据:

x = 0:pi/20:2*pi;
subplot(3,1,1); plot(sin(x))
subplot(3,1,2); plot(cos(x))
subplot(3,1,3); plot(sin(x).*cos(x))
6.8.12 控制轴

axis 命令提供了许多用于设置图形的比例、方向和纵横比的选项。

6.8.13 自动改变坐标轴范围和刻度线

默认情况下,MATLAB 查找数据的最大值和最小值,并选择坐标轴范围来覆盖此范围。MATLAB 选择范围和轴刻度线值,以便生成可清楚地显示数据的图形。但是,您可以使用 axisxlimylimzlim 函数来设置您自己的范围。

注意

更改某根轴的极限会导致其他极限也发生更改,以便更好地表示数据。要禁用自动极限设置,请输入 axis manual 命令。

6.8.14 设置坐标轴范围

axis 命令用于指定您自己的极限:

axis([xmin xmax ymin ymax])

或者对于三维图形,

axis([xmin xmax ymin ymax zmin zmax])

请使用命令

axis auto

重新启用自动极限选择。

6.8.15 设置轴纵横比

axis 命令还可用于指定多种预定义模式。例如,

axis square

使 x 轴和 y 轴的长度相同。

axis equal

使 x 轴和 y 轴上的各个刻度线增量的长度相同。这意味着

plot(exp(1i*(0:pi/10:2*pi)))

(后跟 axis squareaxis equal)会将椭圆形转变为正圆:

axis auto normal

将轴比例恢复为其默认的自动模式。

6.8.16 设置轴可见性

使用 axis 命令可以显示或隐藏轴。

axis on

显示轴。这是默认设置。

axis off

隐藏轴。

6.8.17 设置网格线

grid 命令启用和禁用网格线。语句

grid on

启用网格线,而

grid off

再次禁用网格线。

6.8.18 添加轴标签和标题

此示例展示如何创建图形并增强其显示:

  • 定义 x 和 y 轴的范围 (axis)
  • 对 x 和 y 轴添加标签(xlabelylabel
  • 添加标题 (title)
  • 在图形中添加文本附注 (text)

使用 LaTeX 表示法生成数学符号。

t = -pi:pi/100:pi;
y = sin(t);
plot(t,y)

axis([-pi pi -1 1])
xlabel('-\pi \leq {\itt} \leq \pi')
ylabel('sin(t)')
title('Graph of the sine function')
text(0.5,-1/3,'{\itNote the odd symmetry.}')

如需关于在图形中放置箭头、方框和圆圈的信息,请参阅 annotation 函数。

6.8.19 保存图窗

通过从文件菜单中选择保存来保存图窗。这会将图窗写入到文件,包括属性数据、图窗菜单、uicontrol 和所有注释(即整个窗口)。如果这个图窗以前未被保存过,另存为对话框则会出现。此对话框提供用于将图窗另存为 .fig 文件或将其导出为图形格式的选项。

如果以前保存过这个图窗,再次使用保存会以“静默”方式保存图窗,而另存为对话框不会出现。

要使用标准图形格式(例如,TIFF 或 JPG)保存图窗以便用于其他应用程序,请从文件菜单中选择另存为(如果需要其他控件,则选择导出设置)。

注意

当指定保存图窗的格式时,下次保存该图窗或新图窗时,将再次使用该文件格式。如果您不希望按以前使用的格式保存,请使用另存为,并确保将保存类型下拉菜单设置为要写入的文件类型。

也可通过以下命令行进行保存:

  • 使用 savefig 函数将图窗及其包含的图形对象保存为 .fig 文件。
  • 使用包含任意选项的 saveas 命令,以各种格式保存图窗。
6.8.20 加载图窗

您可以使用以下函数将图窗加载到 MATLAB:

  • 使用 openfig 函数加载保存为 .fig 文件的图窗。
  • 使用 imread 函数将标准图形文件(包括保存图窗)读入到 MATLAB 中。
6.8.21 生成 MATLAB 代码以便再建图窗

通过从图窗文件菜单中选择生成代码,可以生成用于再建图窗及其所包含的图形的 MATLAB 代码。如果您已使用绘图工具创建图形,并且希望使用相同或不同数据创建类似图形,此选项尤其有用。

6.8.22 保存工作区数据

通过从图窗文件菜单中选择将工作区另存为,可以保存工作区中的变量。使用图窗文件菜单中的导入数据项可以重新加载保存的数据。MATLAB 支持多种数据文件格式,包括 MATLAB 数据文件,该数据文件的扩展名为 .mat

6.9 MATLAB 函数绘图函数

注意:内容参考于MathWorks 帮助中心-函数绘图函数

与 MATLAB 基本绘图函数相比,MATLAB 函数绘图函数可用来绘制表达式或者匿名函数的图形,基本绘图函数是用来对数据点进行画图。

fplot 可根据参数函数的变化特性自适应地设置采样间隔。

不推荐使用 ezplot。请改用 fplot。有关详细信息,请参阅兼容性考虑

绘制表达式或函数

fplot
6.9.1 语法
% 在默认区间 `[-5 5]`(对于 `x`)绘制由函数 `y = f(x)` 定义的曲线。
fplot(f)
% 将在指定区间绘图。将区间指定为 `[xmin xmax]` 形式的二元素向量。
fplot(f,xinterval)
% 在默认区间 `[-5 5]`(对于 `t`)绘制由 `x = funx(t)` 和 `y = funy(t)` 定义的曲线。
fplot(funx,funy)
% 将在指定区间绘图。将区间指定为 `[tmin tmax]` 形式的二元素向量。
fplot(funx,funy,tinterval)
% 指定线型、标记符号和线条颜色。例如,`'-r'` 绘制一根红色线条。在前面语法中的任何输入参数组合后使用此选项。
fplot(___,LineSpec)
% 使用一个或多个名称-值对组参数指定线条属性。例如,`'LineWidth',2` 指定 2 磅的线宽。
fplot(___,Name,Value)
% 将图形绘制到 `ax` 指定的坐标区中,而不是当前坐标区 (`gca`) 中。指定坐标区作为第一个输入参数。
fplot(ax,___)
% 返回 `FunctionLine` 对象或 `ParameterizedFunctionLine` 对象,具体情况取决于输入。使用 `fp` 查询和修改特定线条的属性。有关属性列表,请参阅 [FunctionLine 属性](https://ww2.mathworks.cn/help/matlab/ref/matlab.graphics.function.functionline-properties.html) 或 [ParameterizedFunctionLine 属性](https://ww2.mathworks.cn/help/matlab/ref/matlab.graphics.function.parameterizedfunctionline-properties.html)。
fp = fplot(___)
% 返回函数的纵坐标和横坐标,而不创建绘图。在以后的版本中将会删除该语法。请改用线条对象 [`fp`](https://ww2.mathworks.cn/help/matlab/ref/fplot.html#bu6xntj-1-fp) 的 `XData` 和 `YData` 属性。
[x,y] = fplot(___)

注意

fplot 不再支持用于指定误差容限或计算点数量的输入参数。要指定计算点数,请使用 MeshDensity 属性。

6.9.2 绘制表达式

在 x 的默认区间 [-5 5] 绘制 s i n ( x ) sin(x) sin(x)

fplot(@(x) sin(x))
6.9.3 绘制参数曲线

绘制参数化曲线 x = c o s ( 3 t ) x=cos(3t) x=cos(3t) y = s i n ( 2 t ) y=sin(2t) y=sin(2t)

xt = @(t) cos(3*t);
yt = @(t) sin(2*t);
fplot(xt,yt)
6.9.4 指定绘图区间并绘制分段函数

绘制分段函数

e x e^x ex − 3 < x < 0 −3< x <0 3<x<0

c o s ( x ) cos(x) cos(x) 0 < x < 3 0<x<3 0<x<3

使用 hold on 绘制多个线条。使用 fplot 的第二个输入参数指定绘图区间。使用 'b' 将绘制的线条颜色指定为蓝色。在相同坐标区中绘制多个线条时,坐标轴范围会调整以容纳所有数据。

fplot(@(x) exp(x),[-3 0],'b')
hold on
fplot(@(x) cos(x),[0 3],'b')
hold off
grid on
6.9.5 指定线条属性并显示标记

绘制具有不同相位的三个正弦波。对于第一个,使用 2 磅的线宽。对于第二个,指定带有圆形标记的红色虚线线型。对于第三个,指定带有星号标记的青蓝色点划线线型。

fplot(@(x) sin(x+pi/5),'Linewidth',2);
hold on
fplot(@(x) sin(x-pi/5),'--or');
fplot(@(x) sin(x),'-.*c')
hold off
6.9.6 创建后修改线条属性

绘制 sin(x) 并将函数行对象指定给变量。

fp = fplot(@(x) sin(x))
fp = 
  FunctionLine with properties:

     Function: @(x)sin(x)
        Color: [0 0.4470 0.7410]
    LineStyle: '-'
    LineWidth: 0.5000

  Show all properties

通过使用圆点表示法设置属性,将线条更改为红色点线。添加交叉标记,并将标记颜色设置为蓝色。

fp.LineStyle = ':';
fp.Color = 'r';
fp.Marker = 'x';
fp.MarkerEdgeColor = 'b';
6.9.7 在相同坐标区中绘制多个线条

使用 hold on 绘制两个线条。

fplot(@(x) sin(x))
hold on 
fplot(@(x) cos(x))
hold off
6.9.8 添加标题和轴标签以及格式化刻度

使用函数句柄从 −2π 到 2π 绘制 sin(x)。网格线的显示方式。然后添加一个标题,并为 x 轴和 y 轴添加标签。

fplot(@sin,[-2*pi 2*pi])
grid on
title('sin(x) from -2\pi to 2\pi')
xlabel('x');
ylabel('y');

使用 gca 访问当前坐标区对象。沿 x 轴以 π/2 为间隔显示刻度线。通过设置坐标区对象的 XTickXTickLabel 属性,格式化 x 轴刻度值。y 轴存在类似属性。

ax = gca;
ax.XTick = -2*pi:pi/2:2*pi;
ax.XTickLabel = {'-2\pi','-3\pi/2','-\pi','-\pi/2','0',...
    '\pi/2','\pi','3\pi/2','2\pi'};

6.10 MATLAB 三维散点图

参考 MATLAB 三维散点图

三维散点图

scatter3
6.10.1 语法
%% 向量和矩阵数据
scatter3(X,Y,Z) % 在向量 `X`、`Y` 和 `Z` 指定的位置显示圆圈。
scatter3(X,Y,Z,S)% 使用 `S` 指定的大小绘制每个圆圈。要绘制大小相等的圆圈,请将 `S` 指定为标量。要绘制具有特定大小的每个圆,请将 `S` 指定为向量。
scatter3(X,Y,Z,S,C)% 使用 `C` 指定的颜色绘制每个圆圈。
						% - 如果 `C` 是 RGB 三元组,或者是包含颜色名称的字符向量或字符串,则使用指定的颜色绘制所有圆圈。
						% - 如果 `C` 是一个三列矩阵,其中 `C` 中的行数等于 `X`、`Y` 和 `Z` 的长度,则 `C` 的每行指定相应圆圈的 RGB 颜色值。
						% - 如果 `C` 是长度与 `X`、`Y` 和 `Z` 的长度相同的向量,则 `C` 中的值线性映射到当前颜色图中的颜色。
scatter3(___,'filled')% 使用前面的语法中的任何输入参数组合填充这些圆。
scatter3(___,markertype)% 定标记类型。

%% 表数据
scatter3(tbl,xvar,yvar,zvar)% 绘制表 tbl 中的变量 xvar、yvar 和 zvar。要绘制一个数据集,请为 xvar、yvar 和 zvar 各指定一个变量。要绘制多个数据集,请为其中至少一个参数指定多个变量。对于指定多个变量的参数,指定的变量数目必须相同。(自 R2021b 开始提供)
scatter3(tbl,xvar,yvar,zvar,'filled')% 用实心圆绘制表中的指定变量。(自 R2021b 开始提供)

%% 其他选项
scatter3(ax,___)% 将图形绘制到 ax 指定的坐标区中,而不是当前坐标区 (gca) 中。选项 ax 可以位于前面的语法中的任何输入参数组合之前。
scatter3(___,Name,Value)% 通过使用一个或多个名称-值参数设置属性来修改散点图。例如:
							% 有关属性的完整列表,请参阅 Scatter 属性。
							% scatter3(x,y,z,'LineWidth',2) 创建一个标记轮廓为两磅的散点图。
							% scatter3(tbl,'MyX','MyY','MyZ','ColorVariable','MyColors') 根据表中的数据创建一个散点图,并使用表中的数据自定义标记颜色。 
h = scatter3(___) % 返回 Scatter 对象。在创建散点图后,可使用 h 修改其属性。
6.10.2 创建三维散点图

创建三维散点图。使用 sphere 定义向量 xyz

figure
[X,Y,Z] = sphere(16);
x = [0.5*X(:); 0.75*X(:); X(:)];
y = [0.5*Y(:); 0.75*Y(:); Y(:)];
z = [0.5*Z(:); 0.75*Z(:); Z(:)];
scatter3(x,y,z)
6.10.3 改变标记大小

使用 sphere 定义向量 xyz

[X,Y,Z] = sphere(16);
x = [0.5*X(:); 0.75*X(:); X(:)];
y = [0.5*Y(:); 0.75*Y(:); Y(:)];
z = [0.5*Z(:); 0.75*Z(:); Z(:)];

定义向量 s 可指定标记大小。

S = repmat([100,50,5],numel(X),1);
s = S(:);

创建一个三维散点图并使用 view 更改图窗中坐标区的角度。

figure
scatter3(x,y,z,s)
view(40,35)

xyzs 中的相应项确定每个标记的位置和大小。

6.10.4 改变标记颜色

使用 sphere 定义向量 xyz

[X,Y,Z] = sphere(16);
x = [0.5*X(:); 0.75*X(:); X(:)];
y = [0.5*Y(:); 0.75*Y(:); Y(:)];
z = [0.5*Z(:); 0.75*Z(:); Z(:)];

定义向量 sc 以指定每个标记的大小和颜色。

S = repmat([50,25,10],numel(X),1);
C = repmat([1,2,3],numel(X),1);
s = S(:);
c = C(:);

创建一个三维散点图并使用 view 更改图窗中坐标区的角度。

figure
scatter3(x,y,z,s,c)
view(40,35)

xyzc 中的相应项确定每个标记的位置和颜色。

6.10.5 填充标记

将向量 xy 创建为带随机干扰的余弦和正弦值。

z = linspace(0,4*pi,250);
x = 2*cos(z) + rand(1,250);
y = 2*sin(z) + rand(1,250);

创建一个三维散点图并填充标记。使用 view 可更改图窗中坐标区的角度。

scatter3(x,y,z,'filled')
view(-30,10)
6.10.6 设置标记类型

初始化随机数生成器以使 rand 的输出可重复。将向量 xy 定义为带随机干扰的余弦和正弦值。

rng default
z = linspace(0,4*pi,250);
x = 2*cos(z) + rand(1,250);
y = 2*sin(z) + rand(1,250);

创建一个三维散点图并设置标记类型。使用 view 可更改图窗中坐标区的角度。

figure
scatter3(x,y,z,'*')
view(-30,10)
6.10.7 设置标记属性

初始化随机数生成器以使 rand 的输出可重复。将向量 xy 定义为带随机干扰的余弦和正弦值。

rng default
z = linspace(0,4*pi,250);
x = 2*cos(z) + rand(1,250);
y = 2*sin(z) + rand(1,250);

创建一个三维散点图并设置标记边颜色和标记面颜色。使用 view 可更改图窗中坐标区的角度。

figure
scatter3(x,y,z,...
        'MarkerEdgeColor','k',...
        'MarkerFaceColor',[0 .75 .75])
view(-30,10)
6.10.8 绘制表中的数据

自 R2021b 开始提供

绘制表中数据的一种便捷方法是将表传递给 scatter3 函数,并指定要绘制的变量。例如,将 patients.xls 以表 tbl 形式读取。通过将 tbl 作为第一个参数传递给 scatter3 函数,后跟变量名称,绘制 SystolicDiastolicWeight 变量之间的关系。默认情况下,轴标签与变量名称匹配。

tbl = readtable('patients.xls');
scatter3(tbl,'Systolic','Diastolic','Weight');

您也可以同时绘制多个变量。例如,通过将 xvar 参数指定为元胞数组 {'Systolic','Diastolic'},在同一 x 轴上绘制两个血压变量。然后,添加一个图例。图例标签与变量名称匹配。

scatter3(tbl,{'Systolic','Diastolic'},'Age','Weight');
legend
6.10.9 使用自定义标记大小和颜色绘制表数据

自 R2021b 开始提供

绘制表中数据并自定义颜色和标记大小的一种方法是设置 ColorVariableSizeData 属性。您可以在调用 scatter3 函数时将这些属性设置为名称-值参数,也可以稍后在 Scatter 对象上设置它们。

例如,将 patients.xls 以表 tbl 形式读取。用填充标记绘制 SystolicDiastolicWeight 变量之间的关系。通过指定 ColorVariable 名称-值参数来更改标记颜色。将 Scatter 对象返回为 s,以便以后可以设置其他属性。

tbl = readtable('patients.xls');
s = scatter3(tbl,'Systolic','Diastolic','Weight','filled', ...
    'ColorVariable','Diastolic');

通过设置 SizeData 属性,将标记大小更改为 100 磅。然后添加颜色栏。

s.SizeData = 100;
colorbar
6.10.10 指定三维散点图的坐标区

从 R2019b 开始,您可以使用 tiledlayoutnexttile 函数显示分块图。

加载 seamount 数据集以获取向量 xyz。调用 tiledlayout 函数以创建一个 2×1 分块图布局。调用 nexttile 函数以创建坐标区对象 ax1ax2。然后通过将坐标区对象指定为 scatter3 的第一个参数,在坐标区中创建单独的散点图。

load seamount
tiledlayout(2,1)
ax1 = nexttile;
ax2 = nexttile;
scatter3(ax1,x,y,z,'MarkerFaceColor',[0 .75 .75])
scatter3(ax2,x,y,z,'*')
6.10.11 使用句柄设置散点序列属性

使用 sphere 函数创建向量 xyz

[X,Y,Z] = sphere(16);
x = [0.5*X(:); 0.75*X(:); X(:)];
y = [0.5*Y(:); 0.75*Y(:); Y(:)];
z = [0.5*Z(:); 0.75*Z(:); Z(:)];

创建向量 sc 以指定每个标记的大小和颜色。

S = repmat([70,50,20],numel(X),1);
C = repmat([1,2,3],numel(X),1);
s = S(:);
c = C(:);

创建一个三维散点图并返回散点序列对象。

h = scatter3(x,y,z,s,c);

使用 RGB 三元组颜色值设置标记面颜色。使用圆点表示法设置属性。

h.MarkerFaceColor = [0 0.5 0.5];
;

创建一个三维散点图并设置标记类型。使用 view 可更改图窗中坐标区的角度。

figure
scatter3(x,y,z,'*')
view(-30,10)
6.10.7 设置标记属性

初始化随机数生成器以使 rand 的输出可重复。将向量 xy 定义为带随机干扰的余弦和正弦值。

rng default
z = linspace(0,4*pi,250);
x = 2*cos(z) + rand(1,250);
y = 2*sin(z) + rand(1,250);

创建一个三维散点图并设置标记边颜色和标记面颜色。使用 view 可更改图窗中坐标区的角度。

figure
scatter3(x,y,z,...
        'MarkerEdgeColor','k',...
        'MarkerFaceColor',[0 .75 .75])
view(-30,10)
6.10.8 绘制表中的数据

自 R2021b 开始提供

绘制表中数据的一种便捷方法是将表传递给 scatter3 函数,并指定要绘制的变量。例如,将 patients.xls 以表 tbl 形式读取。通过将 tbl 作为第一个参数传递给 scatter3 函数,后跟变量名称,绘制 SystolicDiastolicWeight 变量之间的关系。默认情况下,轴标签与变量名称匹配。

tbl = readtable('patients.xls');
scatter3(tbl,'Systolic','Diastolic','Weight');

您也可以同时绘制多个变量。例如,通过将 xvar 参数指定为元胞数组 {'Systolic','Diastolic'},在同一 x 轴上绘制两个血压变量。然后,添加一个图例。图例标签与变量名称匹配。

scatter3(tbl,{'Systolic','Diastolic'},'Age','Weight');
legend
6.10.9 使用自定义标记大小和颜色绘制表数据

自 R2021b 开始提供

绘制表中数据并自定义颜色和标记大小的一种方法是设置 ColorVariableSizeData 属性。您可以在调用 scatter3 函数时将这些属性设置为名称-值参数,也可以稍后在 Scatter 对象上设置它们。

例如,将 patients.xls 以表 tbl 形式读取。用填充标记绘制 SystolicDiastolicWeight 变量之间的关系。通过指定 ColorVariable 名称-值参数来更改标记颜色。将 Scatter 对象返回为 s,以便以后可以设置其他属性。

tbl = readtable('patients.xls');
s = scatter3(tbl,'Systolic','Diastolic','Weight','filled', ...
    'ColorVariable','Diastolic');

通过设置 SizeData 属性,将标记大小更改为 100 磅。然后添加颜色栏。

s.SizeData = 100;
colorbar
6.10.10 指定三维散点图的坐标区

从 R2019b 开始,您可以使用 tiledlayoutnexttile 函数显示分块图。

加载 seamount 数据集以获取向量 xyz。调用 tiledlayout 函数以创建一个 2×1 分块图布局。调用 nexttile 函数以创建坐标区对象 ax1ax2。然后通过将坐标区对象指定为 scatter3 的第一个参数,在坐标区中创建单独的散点图。

load seamount
tiledlayout(2,1)
ax1 = nexttile;
ax2 = nexttile;
scatter3(ax1,x,y,z,'MarkerFaceColor',[0 .75 .75])
scatter3(ax2,x,y,z,'*')
6.10.11 使用句柄设置散点序列属性

使用 sphere 函数创建向量 xyz

[X,Y,Z] = sphere(16);
x = [0.5*X(:); 0.75*X(:); X(:)];
y = [0.5*Y(:); 0.75*Y(:); Y(:)];
z = [0.5*Z(:); 0.75*Z(:); Z(:)];

创建向量 sc 以指定每个标记的大小和颜色。

S = repmat([70,50,20],numel(X),1);
C = repmat([1,2,3],numel(X),1);
s = S(:);
c = C(:);

创建一个三维散点图并返回散点序列对象。

h = scatter3(x,y,z,s,c);

使用 RGB 三元组颜色值设置标记面颜色。使用圆点表示法设置属性。

h.MarkerFaceColor = [0 0.5 0.5];
  • 1
    点赞
  • 9
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值