向量化代码和预分配--北太天元学习6

向量化代码和预分配

在北太天元学习5中，我们学习了for循环。我们还可以使用for循环来生成
用递推式定义的数列的值。

例: 考虑序列x_{n+1} = x_{n} - 0.2*(x_{n}-25)，x_{0}=37。如果我们想找到 x_1，x_2，...
x_{10}，我们可以使用下面的代码


x = 37;
for n=1:10
    x= x - 0.2*( x - 25)；
end

在循环中的每一步，x的旧值都会被下一个值x - 0.2*(x-25) 所取代。
另一种策略将是对代码进行向量化，这意味着将数列的所有值作为向量进行跟踪。
此时，可以给出如下的代码

x=3；
for i=1:10
    x(i+1) = x(i) - 0.2*( x(i) - 25)；
end
x

这个代码把每一项x_{i+1}的值都保存到变量x中，而不是像上个代码那样简单地保存
最新计算的值。最后一行显示x的所有值。我们还可以用一个命令绘制向量x：

plot(x)

此时plot 仅仅用了一个输入参数，向量x的长度是11,
x的元素分别是 x_{0}, x_{1}, ..., x_{10}
我们使用plot(x), 绘制的横坐标的索引会是从1到11, 也就是 1:11。
由于我们实际上是绘制x_{0}到x_{10}, 希望把横坐标的索引改成0,1,..,10。
要解决此问题，可以使用0:10生成正确的11个x索引(或者叫x坐标）， 然后用

plot(0:10，x)

注意，在for循环的每次迭代中，我们将数列的下一个值附加到向量x上,
此时需要将x的大小(size)增加1。在循环的每一步上，北太天元都会花费时间分配更多
内存到变量x，以便在其中存储更多信息。对于目前这个小规模的例子，这不会带来麻烦。
例如，但如果我们正在执行数百万次迭代，则为每个迭代分配额外内存所需的时间就变得
不可忽略，甚至变得无法忍受。因此，通常优选提前设置向量的大小。
我们可以用下面的命令来执行提前设置向量大小的工作:

x = zeros(11,1）

这个赋值操作调用了内置函数zeros, 把x定义为大小为11x1矩阵(是一个列向量), 且每个元素
都赋值为0. 这里向量长度是11，因为我们记录x_{0}, ..., x_{10}, 这实际上要记录11个值。
通常，命令zeros(m,n) 用于创建大小为mxn的零矩阵 （一个矩阵称为零矩阵，如果它的所有的
元素都是0 ）。 以这种方式定义矩阵或向量的过程称为预分配，并且
如果知道向量或矩阵最终会有多大，则可以使用。通过预先分配变量，
北太天元将把所需的所有计算机内存空间分配好，而不是在程序的整个运行过程中，
不断为越来越大的矩阵分配内存。打个比方，如果我们预先知道在程序运行的整个过程中
需要装入11个鸡蛋，于是我们提前找好了一个能装11个鸡蛋的罐子（预分配内存),然后开始
一个鸡蛋一个鸡蛋撞到罐子里, 这种情况肯定比我们开始仅仅找了一个能装一个鸡蛋的罐子，
然后在运行的过程中，我们发现还需要装一个鸡蛋，于是我又找了一个能装两个鸡蛋的罐子，
此时需要把第一个鸡蛋从刚才能装一个鸡蛋的罐子里转移到刚刚找到的能装两个鸡蛋的罐子，
然后把新的鸡蛋装到了这个能装两个鸡蛋的罐子，但是在后续中发现又要换个更大的罐子，
这个不断换更大罐子的过程肯定是更加浪费时间的。
因此，我们建议如果能够预分配内存，尽量做预分配的工作，因为北太天元为变量分配内存
所花费的时间更少，代码运行速度也会更快。

使用预分配的代码如下所示：

x = zeros(11,1);
x(1) = 37;
for i = 1:10
    x(i+1) = x(i) - 0.2*( x(i) - 25)；
end

plot(0:10，x）

xlabel('n')

ylabel('x_{n}')

您不会注意到这个预先分配x的代码和之前的代码有什么区别，因为这两个版本的代码
看起来都是即时运行的。然而，如果问题的规模很大，预分配将会带来不同。

我们再举一个例子展示一下使用嵌套的for循环。我们来推广一个刚才的例子，
刚才的例子可以解释为一个物体的温度随着时间不断降低的过程，
x_{n+1} = x_{n} - k ( x_{n} - 25)
表示的是 n+1 时刻的温度比 n 时刻的温度下降的温度 和 n时刻的温度与环境温度25之间的
差成正比，这个比例系数k 称为冷却率。我们可以计算不同冷却率k的温度变换情况.
考虑序列x_{n+1} = x_{n} - k*(x_{n} - 25)，x_{0}=37， n=1，2，··，10，
其中k是0.05、0.01、0.15、0.02之一。 对于每一个k的值，我们要计算x_{n}, n=1,...,10.
我们把不同冷却率k的计算结果放在不同的列, 计算的结果存在 11x4的矩阵，
每一列对应一个k值.

x = zeros(11, 4);
x(1,:)=37;
k=[0.05 0.01 0.15 0.2];

for i = 1:4
  for j = 1:10
        x(j+1,i) = x(j,i) - k(i)*( x(j,i) - 25);
    end
end
    
plot(0:10, x(:,1), 0:10, x(:,2), 0:10, x(:,3), 0:10, x(:,4));
xlabel('n');
ylabel('x_{n}')
legend('k=0.05', 'k=0.1', 'k=0.15', 'k=0.2')
x(end,:)
我们可以从图中看出，改变冷却率k给出了x在n=10的值介于约
26到32之间，冷却率k越大，x_{10} 越小。
在代码的最后一行使用命令x(end,:)输出x的最后一行的精确值.