Z 変換を使用した差分方程式の求解

このワークフローでは、Symbolic Math Toolbox™ で Z 変換を使用して差分方程式を解きます。Z 変換の単純な例については、ztransおよびiztransを参照してください。

定義: Z 変換

関数 f(n) の Z 変換は、次のように定義されます。

$F (z) = \sum_{n = 0}^{\infty} \frac{f (n)}{z^{n}} .$

概念: シンボリックワークフローの使用

シンボリックワークフローは、計算を数値形式の代わりにもともとのシンボリック形式で保持します。この手法は、解のプロパティについて理解し、厳密なシンボリック値を使用するのに役立ちます。数値結果が必要な場合や、シンボリック的に続行できない場合に限り、シンボリック変数に数値を代入します。詳細は、数値演算またはシンボリック演算の選択を参照してください。通常、手順は以下になります。

方程式を宣言する。
方程式を解く。
値を代入する。
結果をプロットする。
結果を解析する。

ワークフロー: Z 変換を使用した "うさぎの増加" の問題の求解

方程式の宣言

Z 変換を使用すると、よく知られた "うさぎの増加" の問題 (ネズミ算) などの差分程式を解くことができます。1 組のうさぎが 1 年で成長して毎年 1 組のうさぎを産む場合、年 $n$ におけるうさぎの数 $p (n)$ は以下の差分方程式で表されます。

$p (n + 2) = p (n + 1) + p (n)$

この方程式を、右辺を 0 とする形式の式として表します。n は年を表すため、n は正の整数であると仮定します。この仮定により、結果が単純化されます。

syms p(n) z
assume(n>=0 & in(n,"integer"))
f = p(n+2) - p(n+1) - p(n)

f =  $p (n + 2) - p (n + 1) - p (n)$

方程式の求解

方程式の Z 変換を求めます。

fZT = ztrans(f,n,z)

fZT =  $z p (0) - z ztrans (p (n), n, z) - z p (1) + z^{2} ztrans (p (n), n, z) - z^{2} p (0) - ztrans (p (n), n, z)$

関数 solve は、シンボリック変数のみを解きます。したがって、solve を使用するには、最初に ztrans(p(n),n,z) を変数 pZT に置き換えます。

syms pZT
fZT = subs(fZT,ztrans(p(n),n,z),pZT)

fZT =  $z p (0) - pZT - z p (1) - pZT z - z^{2} p (0) + pZT z^{2}$

pZT の解を求めます。

pZT = solve(fZT,pZT)

pZT = 
 $- \frac{z p (1) - z p (0) + z^{2} p (0)}{- z^{2} + z + 1}$

pZT の逆 Z 変換を計算して、 $p (n)$ を計算します。結果を単純化します。

pSol = iztrans(pZT,z,n);
pSol = simplify(pSol)

pSol = 
 $\begin{array}{l} 2 {(- 1)}^{n / 2} p (1) \cos (n (\frac{π}{2} + asin (\frac{1}{2} i))) + \frac{2^{2 - n} \sqrt{5} {(\sqrt{5} + 1)}^{n - 1} σ_{1}}{5} - \frac{2 2^{1 - n} \sqrt{5} {(1 - \sqrt{5})}^{n - 1} σ_{1}}{5} \\ where \\ σ_{1} = \frac{p (0)}{2} - p (1) \end{array}$

値の代入

結果をプロットするには、最初に初期条件の値を置き換えます。p(0) と p(1) をそれぞれ 1 と 2 とします。

pSol = subs(pSol,[p(0) p(1)],[1 2])

pSol = 
 $4 {(- 1)}^{n / 2} \cos (n (\frac{π}{2} + asin (\frac{1}{2} i))) - \frac{3 2^{2 - n} \sqrt{5} {(\sqrt{5} + 1)}^{n - 1}}{10} + \frac{3 2^{1 - n} \sqrt{5} {(1 - \sqrt{5})}^{n - 1}}{5}$

結果のプロット

pSol プロットすることで、時間と共に増加するうさぎの数を示します。

nValues = 1:10;
pSolValues = subs(pSol,n,nValues);
pSolValues = double(pSolValues);
pSolValues = real(pSolValues);
stem(nValues,pSolValues)
title("Rabbit Population")
xlabel("Years (n)")
ylabel("Population p(n)")
grid on

結果の解析

このプロットは、解が指数関数的に増加する様子を示しています。しかし、解 pSol には多くの項が含まれているため、この動作を生み出す項を求めるには解析が必要です。

pSol のすべての関数は exp で表すことができるため、pSol を exp に書き換えます。追加の単純化ステップ 80 を指定した simplify を使用し、結果を単純化します。これで、pSol を解析できます。

pSol = rewrite(pSol,"exp");
pSol = simplify(pSol,"Steps",80)

pSol = 
 $\frac{2 2^{n}}{{(\sqrt{5} - 1)}^{n}} + \frac{2 {(1 - \sqrt{5})}^{n}}{2^{n}} - \frac{6 \sqrt{5} {(\sqrt{5} + 1)}^{n}}{5 2^{n} (\sqrt{5} + 1)} - \frac{6 \sqrt{5} {(1 - \sqrt{5})}^{n}}{5 2^{n} (\sqrt{5} - 1)}$

pSol を目視で検査します。pSol は項の和であることに注意してください。各項は、n の増加と共に増加または減少する可能性がある比です。各項について、複数の方法でこの仮定を確認できます。

limit を使用して n = Inf での極限が 0 になるのか、それとも Inf になるのかを確認します。
増加する n に対して項をプロットして動作を確認します。
大きい n の値で値を計算します。

簡単にするために、3 番目の方法を使用します。n = 100 で項を計算し、この方法を検証します。最初に children を使用して個々の項を求め、n に置き換えて、double に変換します。

pSolTerms = children(pSol);
pSolTerms = [pSolTerms{:}];
pSolTermsDbl = subs(pSolTerms,n,100);
pSolTermsDbl = double(pSolTermsDbl)

pSolTermsDbl = 1×4
10²¹ ×

    1.5841    0.0000   -0.6568   -0.0000

結果は、0 の項がある一方で大きく振れる項があることを示しています。大きい値の項によって指数関数的な動作が発生するとの仮説を立てます。これらの項で pSol を近似します。

idx = abs(pSolTermsDbl)>1; % use arbitrary cutoff
pApprox = pSolTerms(idx);
pApprox = sum(pApprox)

pApprox = 
 $\frac{2 2^{n}}{{(\sqrt{5} - 1)}^{n}} - \frac{6 \sqrt{5} {(\sqrt{5} + 1)}^{n}}{5 2^{n} (\sqrt{5} + 1)}$

pSol と pApprox 間の近似誤差をプロットして仮説を検証します。予想どおり、誤差は n の増加と共に 0 に近づきます。この結果は、シンボリック計算が問題の解析にどのように役立つかを示しています。

Perror = pSol - pApprox;
nValues = 1:30;
Perror = subs(Perror,n,nValues);
stem(nValues,Perror)
xlabel("Years (n)")
ylabel("Error (pSol - pApprox)")
title("Error in Approximation")