制御工学の計算で極を指定して特性方程式を求める際に必要となったのでここにメモとして残します。
極を[-1, -3, -2, -4, -4.5, -2]とした場合、MATLABでは
s = tf('s');
(s+1)*(s+3)*(s+2)*(s+4)*(s+4.5)*(s+2)
% 計算結果
>> s^6 + 16.5 s^5 + 109 s^4 + 367.5 s^3 + 664 s^2 + 606 s + 216のようにして計算できますが、これをシンボリック計算を使わずに計算するにはどうすれば良いでしょうか。規則性を探るために、\( (x+a) \)から\( (x+a)(x+b)… \)と増やしていって4次式になるまで展開してみます。
$$ (x + a)(x+ b) = x^2 + (a + b) x + ab $$
$$ (x + a)(x+ b)(x+c) = x^3 + (a + b + c)x^2 + \{ ab + (a + b)c \} x + abc $$
$$\begin{align} (x + a)(x+ b)(x+c)(x + d) &= x^4 + (a + b + c+ d) x^3 \\ &\quad+ [\{ ab+(a + b)c \} + (a + b + c)d] x^2 \\ &\quad+ [abc + \{ ab + (a + b)c \}d] + abcd \end{align}$$
ちょっと見づらいですが、上のように展開した式を見ると以下のような規則性があることがわかります。
\( x^n \)の係数 ＝ 一つ前の式の\( x^{n-1} \)の係数 ＋ 一つ前の式の\( x^n \)の係数 × 新しい定数項
例えば、\( (x+a)(x+b)(x+c) \)を展開した時の\(x\)の係数を計算する場合はそれぞれ
- 一つ前の式の\( x^0 \)の係数：\( ab \)
- 一つ前の式の\( x \)の係数：\( (a + b) \)
- 新しい定数項：\( c \)
となるので、\(x\)の係数は\( ab + (a + b)c \)となります。
規則性がわかったので、Pythonで実装してみました。
# (x + a)(x + b)(x + c)... のような式を展開した時のxの係数を求める
#
# x^n + k_(n-1) x^(n-1) + ... + k_1 x + k_0
import copy
def expand(nums):
# 計算結果を格納する配列（1で初期化しておく）
result = [1.0] * len(nums) # [k_0, k_1, ..., k_(n-1)]
result[0] = nums[0]
for i in range(1, len(result)):
pre_coef = copy.copy(result)
result[0] = pre_coef[0] * nums[i]
for j in range(1, i + 1):
result[j] = pre_coef[j-1] + pre_coef[j] * nums[i]
return result
if __name__ == '__main__':
nums = [1.0, 3.0, 2.0, 4.0, 4.5, 2.0]
print(expand(nums))これを実行すると
[216.0, 606.0, 664.0, 367.5, 109.0, 16.5]と表示され、最初にMATLABで計算した値と等しいことが確認できます（\(x^n\)の係数は必ず1なので計算結果には含まれません）。