使用z -变换求解差分方程- MATLAB & Simulink - MathWorks澳大利亚金宝app

用z变换求解差分方程

通过使用符号数学工具箱™中的z变换来解决差分方程。有关z变换的简单示例，请参见ztrans而且iztrans．

定义:z变换

函数的z变换f（n)的定义为

$F （ z ）＝ \sum_{n ＝ 0}^{\infty} \frac{f （ n ）}{z^{n}} ．$

概念:使用符号工作流

符号工作流保持自然符号形式的计算，而不是数字形式。这种方法可以帮助您理解解决方案的属性，并使用精确的符号值。只有在需要数值结果或不能以符号形式继续时，才可以用数字代替符号变量。详细信息请参见选择数字或符号算术．通常，这些步骤是:

申报的方程。
解决方程。
替代的价值观。
阴谋的结果。
分析的结果。

工作流程:用z -变换解决“兔子生长”问题

申报的方程

你可以使用z变换来解差分方程，比如著名的“兔子生长”问题。如果一对兔子在一年内成熟，然后每年生产另一对兔子，兔子的数量p（n）在一年n由这个差分方程描述。

p（n+ 2) =p（n+ 1) +p（n）.

将方程声明为表达式，假设右侧为0．因为n代表年份，假设n为正整数。这个假设简化了结果。

信谊p (n) z假设(n > = 0 & (n,“整数”))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 = * p (0) - z * ztrans (p (n), n, z) - z * p (1) + z ^ 2 * ztrans (p (n), n, z) -…Z²*p(0) - ztrans(p(n) n, Z)

这个函数解决只求解符号变量。因此，使用解决，第一个代换ztrans (p (n), n, z)有了变量压电．

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

fZT = z * p(0) -压电陶瓷片的压电陶瓷- z * p (1) * z - z ^ 2 * p(0) +压电* z ^ 2

解出压电．

pZT = solve(fZT,pZT)

压电陶瓷= - * p (1) - (z z * p (0) + z ^ 2 * p (0)) / (z - z ^ 2 + + 1)

计算p（n）通过计算的z反变换压电．简化结果。

pSol = iztrans(pZT,z,n);pSol =简化(pSol)

pSol = 2 * (1) ^ (n / 2) * p (1) * cos (n *(π/ 2 + asin我/ 2)(1))+…(2 ^ (2 - n) * 5 ^ (1/2) * (5 ^ (1/2) + 1) ^ (n - 1) * (p (0) / 2 - p(1))) / 5 -…(2 * 2 ^ (1 - n) * 5 ^ (1/2) * (1 - 5 ^ (1/2)) ^ (n - 1) * (p (0) / 2 - p (1))) / 5

替代值

要绘制结果，首先替换初始条件的值。让p (0)而且(1页)是1而且2,分别。

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

pSol = 4 * (1) ^ (n / 2) * cos (n *(π/ 2 + asin我/ 2)(1))-…(3 * 2 ^ (2 - n) * 5 ^ (1/2) * (5 ^ (1/2) + 1) ^ (n - 1) / 10 +……3*2^(1 - n)*5^(1/2)*(1 - 5^(1/2))^(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可以用经验值,重写pSol来经验值．通过使用简化结果简化与80附加的简化步骤。现在，你可以分析了pSol．

pSol =重写(pSol，'exp');pSol = simplify(pSol，'Steps'，80)

pSol = (2 * 2 ^ n) / (5 ^ (1/2) - 1) ^ n + (2 * (1 - 5 ^ (1/2)) ^ n) / 2 ^ n -…(6 * 5 ^ (1/2) * (5 ^ (1/2) + 1) ^ n) / (5 * 2 ^ n *(5 ^(1/2) + 1)) -…(6*5^(1/2)*(1 - 5^(1/2))^n)/(5*2^ (5^(1/2) - 1))

视觉检查pSol．请注意,pSol是项的和。每一项都是一个可以增加或减少的比率n增加。对于每一项，你可以用几种方式来证实这个假设:

检查极限是否在n =无穷去0或正通过使用限制．
画出递增项n检查行为。
的较大值处计算值n．

为了简单起见，使用第三种方法。计算在N = 100，然后验证该方法。首先，通过使用找到各个术语孩子们，代替n，并转换为double。

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

pSolTermsDbl = 1.0e+21 * 1.5841 0.0000 -0.6568 -0.0000

结果表明，某些项是正确的0而其他项的幅度很大。假设大量值项产生指数行为。近似pSol用这些项。

idx = abs(pSolTermsDbl)>1;%使用任意截断pApprox = pSolTerms(idx);pApprox = sum(pApprox)

pApprox = (2*2^n)/(5^(1/2) - 1)^n -…(6*5^(1/2)*(5^(1/2) + 1)^n)/(5*2^ (5^(1/2) + 1))

通过绘制之间的近似误差来验证假设pSol而且pApprox．正如预期的那样，错误变为0作为n增加。这个结果演示了符号计算如何帮助您分析问题。

Perror = pSol - pApprox;nValues = 1:30;Perror = subs(Perror,n,nValues);stem(nValues,Perror) xlabel('Years (n)') ylabel('错误(pSol - pApprox)') title('错误在近似')

参考文献

安德鲁，l.c.， Shivamoggi, b.k.，工程师和应用数学家的积分变换，麦克米伦出版公司，纽约，1986年

克兰德尔，r.e.，科学计算项目，斯普林格出版社，纽约，1994年

G.斯特朗，应用数学概论，威尔斯利-剑桥出版社，马萨诸塞州威尔斯利，1986年