向量化代碼和預分配

在北太天元學習5中，我們學習了for循環(huán)。我們還可以使用for循環(huán)來生成
用遞推式定義的數(shù)列的值。

例: 考慮序列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

在循環(huán)中的每一步，x的舊值都會被下一個值x - 0.2*(x-25) 所取代。
另一種策略將是對代碼進行向量化，這意味著將數(shù)列的所有值作為向量進行跟蹤。
此時，可以給出如下的代碼

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 僅僅用了一個輸入?yún)?shù)，向量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循環(huán)的每次迭代中，我們將數(shù)列的下一個值附加到向量x上,
此時需要將x的大小(size)增加1。在循環(huán)的每一步上，北太天元都會花費時間分配更多
內存到變量x，以便在其中存儲更多信息。對于目前這個小規(guī)模的例子，這不會帶來麻煩。
例如，但如果我們正在執(zhí)行數(shù)百萬次迭代，則為每個迭代分配額外內存所需的時間就變得
不可忽略，甚至變得無法忍受。因此，通常優(yōu)選提前設置向量的大小。
我們可以用下面的命令來執(zhí)行提前設置向量大小的工作:

x = zeros(11,1）

這個賦值操作調用了內置函數(shù)zeros, 把x定義為大小為11x1矩陣(是一個列向量), 且每個元素
都賦值為0. 這里向量長度是11，因為我們記錄x_{0}, ..., x_{10}, 這實際上要記錄11個值。
通常，命令zeros(m,n) 用于創(chuàng)建大小為mxn的零矩陣 （一個矩陣稱為零矩陣，如果它的所有的
元素都是0 ）。 以這種方式定義矩陣或向量的過程稱為預分配，并且
如果知道向量或矩陣最終會有多大，則可以使用。通過預先分配變量，
北太天元將把所需的所有計算機內存空間分配好，而不是在程序的整個運行過程中，
不斷為越來越大的矩陣分配內存。打個比方，如果我們預先知道在程序運行的整個過程中
需要裝入11個雞蛋，于是我們提前找好了一個能裝11個雞蛋的罐子（預分配內存),然后開始
一個雞蛋一個雞蛋撞到罐子里, 這種情況肯定比我們開始僅僅找了一個能裝一個雞蛋的罐子，
然后在運行的過程中，我們發(fā)現(xiàn)還需要裝一個雞蛋，于是我又找了一個能裝兩個雞蛋的罐子，
此時需要把第一個雞蛋從剛才能裝一個雞蛋的罐子里轉移到剛剛找到的能裝兩個雞蛋的罐子，
然后把新的雞蛋裝到了這個能裝兩個雞蛋的罐子，但是在后續(xù)中發(fā)現(xiàn)又要換個更大的罐子，
這個不斷換更大罐子的過程肯定是更加浪費時間的。
因此，我們建議如果能夠預分配內存，盡量做預分配的工作，因為北太天元為變量分配內存
所花費的時間更少，代碼運行速度也會更快。

使用預分配的代碼如下所示：

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的代碼和之前的代碼有什么區(qū)別，因為這兩個版本的代碼
看起來都是即時運行的。然而，如果問題的規(guī)模很大，預分配將會帶來不同。

我們再舉一個例子展示一下使用嵌套的for循環(huán)。我們來推廣一個剛才的例子，
剛才的例子可以解釋為一個物體的溫度隨著時間不斷降低的過程，
x_{n+1} = x_{n} - k ( x_{n} - 25)
表示的是 n+1 時刻的溫度比 n 時刻的溫度下降的溫度 和 n時刻的溫度與環(huán)境溫度25之間的
差成正比，這個比例系數(shù)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的最后一行的精確值.