matlab两个for循环嵌套加速,【MATLAB】并行计算初探！用parfor改写嵌套循环提高速度！...

咳咳……嗯，好，我们研究个问题——阻尼弹簧振子模型

我们先用牛二定律把无外力驱动情况下，阻尼弹簧振子的振动方程写出来！

设一个弹簧振子的质量为m，弹簧劲度系数为k，弹簧振子在有阻尼f=-bv情况下振动

根据牛顿第二定律，写出弹簧振子的运动定律：

移项即可得到

其中b称为阻尼系数，k则为劲度系数，工程上称为刚度系数

我们要研究这个二阶微分方程伴随k、b的动响应最大值是怎样的。也就是说，要研究在不同的k值、b值下，x的最大值是多少。物理图景角度讲，就是要找到振子偏离平衡位置所达到的最远距离伴随阻尼系数、刚度系数的变化。

可以预先设想，如果进行数据绘图，得到的结果是应该是三维曲面，曲面方程是这样的形式：

我们首先来想最简单的计算实现形式。大家可能会想到，对于一个双变量的函数最大值，可以使用for嵌套for的方式实现，其思想本源就是“控制变量”。内层for循环可以在k不变时遍历b的值，进而计算x的最大值，外层for循环可以遍历每一个k值(或者把关于k、b的循环反过来也行)，这样就可以遍历每一个k、b的值，计算出x的最大值了。至于微分方程的求解，可以直接利用ODE来求解，应该不是很困难。我们不妨按着这个思路来实现一下！

首先确定参数取值范围，在此我尝试了一些取值：

参数的生成可以使用linspace指令，写出来就是这样：

% MATLAB Code

ccc;%三清

m=5;%质量

b=linspace(0.1,5,50);%阻尼系数

k=linspace(1.5,5,70);%刚度系数

%——————

接下来需要构造ODE函数，来求这个二阶微分方程。按照MATLAB给出的odefun的一般构造方法，即降阶为一阶方程的方法来进行：

令

有

考虑到参数传递，为了避免麻烦，就采用了两层function嵌套的方式，外层fucntion用于传回ode求解结果和传参，内层function作为odefun：

% MATLAB Code

function [t,y]=sol_ode(m,b,k)

%传参嵌套方法

[t,y]=ode45(@odefun,[0,25],[0;1]);

function r=odefun(t,y)

%ODE函数

r=[y(2);(-b.*y(2)-k.*y(1))./m];

end

end

%——————

考虑到伴随循环的进行，竖坐标数据的增多，所以先行定义一个矩阵用于存储竖坐标的数据。一般来说，可以用nan指令或zeros指令定义非数矩阵或全零矩阵，注意其大小要与b、k的大小匹配：

% MATLAB Code

x_max=nan(length(b),length(k));%x最大值

%——————

下面开始构造循环，并将每次循环所得的x值填入x_max矩阵中，注意x的值是y矩阵的第一列，即x=y(:,1)

% MATLAB Code

for ik=1:length(k)%外层循环遍历k

for ib=1:length(b)%内层循环遍历b

[t,y]=sol_ode(m,b(ib),k(ik));%求解方程

x_max(ib,ik)=max(y(:,1));%填充数据

end

end

%——————

这样，就求得了不同b、k值下x的最大值，并将结果存入了x_max矩阵中

下面进行图形绘制，要知道，直接使用“surf(b,k,x_max)”是不能绘图的，因为b、k不是二维矩阵，而是向量数组，surf不支持用向量数组进行绘图。因此，我们用meshgrid进行数据网格构建后再绘图。

% MATLAB Code

[k,b]=meshgrid(k,b);%构建数据网格用于绘图

surf(b,k,x_max)%绘制曲面

shading interp%着色

%——————

可以看到，得到的图形是这样的：

完整代码：

% MATLAB Code

ccc;%三清

m=5;%质量

b=linspace(0.1,5,50);%阻尼系数

k=linspace(1.5,5,70);%刚度系数

x_max=nan(length(b),length(k));%x最大值

for ik=1:length(k)%外层循环遍历k

for ib=1:length(b)%内层循环遍历b

[t,y]=sol_ode(m,b(ib),k(ik));%求解方程

x_max(ib,ik)=max(y(:,1));%填充数据

end

end

[k,b]=meshgrid(k,b);%构建数据网格用于绘图

surf(b,k,x_max)%绘制曲面

shading interp%着色

function [t,y]=sol_ode(m,b,k)

%传参嵌套方法

[t,y]=ode45(@odefun,[0,25],[0;1]);

function r=odefun(t,y)

%ODE函数

r=[y(2);(-b.*y(2)-k.*y(1))./m];

end

end

%——————

其实，for-for嵌套的方法对于计算更多的数据时会显得十分吃力，一般来说，会采用矢量化编程的方法，也就是在程序最初确定b、k的值时就采用meshgrid生成数据网格，而后便可去除for循环。当然，odefun也需要做相应的改变，因为b、k不再是单独的一个数，而是向量或矩阵。大家可以尝试哦！

(我表示有点南，就不写了……)

但是，这不是今天的重点！今天的重点在这！

为大家介绍一个for-for循环的改进方法：

并行for循环！

什么是并行for循环呢？说的通俗一些，就是把一个for循环拆成若干部分，分别让不同的“工人(Worker)”去完成。相当于是：普通for循环是一个人按顺序做完一件事，并行for循环是多个人分块做完一件事。理论上讲，并行for循环会比普通for循环快一些。我们来做个尝试，比较两者的快慢。

并行for循环的关键字是“parfor”，需要注意的是，并行for循环有很多“讲究”，比如循环中不能出现循环变量是小数的情况，循环中不能出现迭代等。我们的程序恰好满足了这些条件，没有变量出现小数，没有循环的迭代，因此可以进行修改！

需要修改的位置是：

1.将外层的for改为parfor(注意MATLAB暂时不能parfor-parfor嵌套)

2.将内层循环的ib=1:length(b)修改为ib=1:bNum，bNum需要在parfor前声明为bNum=length(b);

修改部分代码：

% MATLAB Code

bNum=length(b);

parfor ik=1:length(k)%外层循环遍历k

for ib=1:bNum%内层循环遍历b

[t,y]=sol_ode(m,b(ib),k(ik));%求解方程

x_max(ib,ik)=max(y(:,1));%填充数据

end

end

%——————

这样就修改完成了。大家会看到，在第二个修改中，需要将内层for循环中的length函数部分替换为定值，这是因为“ib”是parfor循环中的嵌套循环变量。MATLAB在为每个Worker“分配”任务时，必须将任务中的变量确定，防止任务出错。因此内层循环中要求循环变量具有固定的范围。

当然，将内层for循环改为parfor循环就不需要上面第二步的修改，代码：

% MATLAB Code

for ik=1:length(k)%外层循环遍历k

parfor ib=1:length(b)%内层循环遍历b

[t,y]=sol_ode(m,b(ib),k(ik));%求解方程

x_max(ib,ik)=max(y(:,1));%填充数据

end

end

%——————

我们来对比一下for-for嵌套和parfor-for循环的时间，不妨将b、k划分的点数多一些，各改为500、700进行测试！在for循环上下分别加上tic、toc即可实现计时：

对比代码1(for-parfor循环)：

% MATLAB Code

tic

for ik=1:length(k)%外层循环遍历k

parfor ib=1:length(b)%内层循环遍历b

[t,y]=sol_ode(m,b(ib),k(ik));%求解方程

x_max(ib,ik)=max(y(:,1));%填充数据

end

end

toc

%——————

结果：历时 67.272047 秒。

对比代码2(for-for循环)：

% MATLAB Code

tic

for ik=1:length(k)%外层循环遍历k

for ib=1:length(b)%内层循环遍历b

[t,y]=sol_ode(m,b(ib),k(ik));%求解方程

x_max(ib,ik)=max(y(:,1));%填充数据

end

end

toc

%——————

结果：历时 169.156309 秒。

很明显，并行for循环会比for-for嵌套快一些！

但是，如果要处理更多的数据，运行还是很慢，要等很久……这时有同志会问：既然“包工头”能够让多个“工人”同时跟自己工作，那么能不能让“包工头”不干活、歇着呢？也就是说能不能让MATLAB后台计算一个程序，前台运行其他程序呢？就比如上面的程序，让MATLAB在后台运行它，释放前台并在前台运行其他程序，等到后台运行结束后回收数据。这个能不能实现？当然可以实现！

留个坑，下一期写！

先安装Toolbox去了……