当当触屏版

get_product_contenthtml 第5 章积分变换与复变函数问题的求解积分变换技术可以将某些难以分析的问题通过映射的方式映射到其他域内的表达式后再进行分析。例如，Laplace变换可以将时域函数映射成复域函数，从而可以将某时域函数的微分方程映射成复域的多项式代数方程，使得原微分方程在诸多方面，如稳定性、解析解等方面更便于分析，这样的变换方法构成了经典自动控制理论的基础。在实际应用中，Fourier变换、Mellin变换及Hankel变换都是有其应用领域的。如何利用计算机求解积分变换的解析解是本章主要介绍的问题之一。如果读者没有学过积分变换与复变函数课程，也可以利用类似于第3章介绍的方法，直接由计算机求解相关问题。5.1 节将首先介绍Laplace变换与反变换的定义及基本性质，然后介绍用MATLAB语言中的符号运算工具箱函数求取Laplace变换及反变换问题的解析解方法，还给出了复杂函数Laplace反变换的数值求解方法与应用实例。5.2节将介绍Fourier变换及反变换的定义、性质和变换问题的MATLAB解法，并介绍Fourier余弦变换、正弦变换、离散Fourier正余弦变换等问题的计算机求解方法，并介绍快速Fourier变换的求解与应用。5.3节将介绍Mellin变换、Hankel变换等问题的MATLAB语言的求解算法，可以得出函数的相应变换及反变换。z变换是另一类实用的变换方法，该变换方法也是离散控制理论的数学基础。5.4节将介绍z变换及其反变换的定义和性质，并介绍基于MATLAB语言符号运算工具箱的z变换问题的计算机辅助求解方法。本章的另一个主要问题是复变函数问题及其MATLAB语言求解，可以用5.5节中介绍的方法计算复变函数的奇点与留数，进行部分分式展开等运算，讨论了有理函数Laplace反变换的求解方法和化简方法，基于留数定理还探讨了封闭曲线积分的求解方法。5.7节还将介绍各种差分方程的求解方法。5.1 Laplace变换及其反变换法国数学家Pierre-SimonLaplace（1749–1827）引入的积分变换可以巧妙地将一般常系数微分方程映射成代数方程，奠定了很多领域，如电路分析、自动控制原理等的数学模型基础。本节将首先介绍Laplace变换及其反变换的定义与性质，然后介绍利用计算机数学语言MATLAB求解Laplace变换及其反变换的方法与应用，后给出复杂函数Laplace反变换的数值求解方法与实用函数。5.1.1 Laplace变换及反变换的定义与性质一个时域函数f(t)的Laplace变换可以定义为∫L[f(t)]=∞ f(t)e.stdt=F(s)（5-1-1）式中，L[f(t)]为Laplace变换的简单记号。0下面将不加证明地列出一些常用的Laplace变换性质。· 150·高等应用数学问题的MATLAB求解（第四版）（1）线性性质。若a与b均为标量，则L[af(t)± bg(t)]=aL[f(t)]± bL[g(t)]。
（2）时域平移性质。L[f(t. a)]=e.asF(s)。（3）s-域平移性质[。L[e.atf](t)]=F(s a)。
（4）微分性质。df(t)/dt=sF(s). f(0 )，一般地，n阶微分可以由下式求出

ndnL [f(t)/dt]L= s nF(s). s n.1f(0 ).s n.2f′(0 ) .···. f(n.1)(0 )（5-1-2）若假设函数f(t)及其各阶导数的初值均为0，则式（5-1-2）可以简化成n此性质事实上是微分方程映射成代数方程的关键式子。t （5）积分性质。若假设零初始条件，L0 f(τ)dτ=Fs (s)，一般地，函数f(t)的n重积分的Laplace变换可以由下式求出L [∫[t dnf(∫t[t )/∫dt] = ] s n] F (Fs() s) （5-1-3）Lf(τ)dτn=（5-1-4）··· sn 00 （6）初值性质。tlim 0 f(t)=limsF(s)。s→∞（7）终值性质。F(s)没有s.0的极点，则limf(t)=limsF(s)。如果→t→∞ s→0 （8）卷积性质。L[f(t). g(t)]=L[f(t)]L[g(t)]，式中，卷积算子. 的定义为∫t ∫t f(t). g(t)=f(τ)g(t. τ)dτ=f(t. τ)g(τ)dτ（5-1-5）00 （9）其他性质。[] ∫∫nL[tnf(t)]=(.1)n dndFsn (s), L ft(n t)= ∞ ··· ∞ F(s)ds（5-1-6）ss 如果已知函数的Laplace变换表达式F(s)，则可以通过下面的反变换公式反演求出其Laplace反变换f(t)=L.1[F(s)]=1 ∫λ j∞ F(s)estds（5-1-7）2.jλ.j∞ 其中，σ大于F(s)的奇点实部，奇点的概念将在后面给出。5.1.2 Laplace变换的计算机求解求解已知函数的Laplace变换的解析解必须借助于计算机数学语言，如本书介绍的MATLAB语言。MATLAB语言的符号运算工具箱可以轻松地求解Laplace变换问题。具体的变换及反变换问题的求解步骤为：（1）用syms命令声明符号变量t，这样就能描述时域表达式f了。
（2）直接调用laplace()函数，就可以得出所需的时域函数Laplace变换式子

还可以考虑采用simplify()等函数对其进行化简。（3）对复杂的问题来说，得出的结果形式通常需要调用simplify()函数化简，此外，有时变换的结果难以阅读，所以需要调用pretty()函数或latex()函数对结果进行进一步处理。可以在屏幕上或利用LATEX的强大功能将结果用可读性更强的形式显示出来。如果已知Laplace变换式子，则应该首先给出Laplace变换式子F，然后采用符号运算工具箱中的ilaplace()函数对其进行反变换。该函数的调用格式为
获得变化式f 后也可以对之进一步化简和改变显示格式。例5-1已知函数f(t)=t2e.2t sin(t .)，试求取该函数的Laplace变换。解分析原题，可以先声明t为符号变量，再用MATLAB语句表示给定的f(t)函数，然后就可以用下面的语句对该函数进行Laplace变换
2 直接得出原函数的Laplace变换为F(s)=(2 )2 . (2(2s 4))3 。22(s 2) 1(s 2) 1例5-2假设给出的函数为f(x)=x2e.2x sin(x .)，试求其Laplace变换，并对结果进行Laplace反变换，看是否能变换回原函数。解同样可以采用laplace()函数求解该问题
可见，这样的结果和上面的例子是完全一致的，但需要按要求做一下变量替换。使用Laplace反变换的函数ilaplace(F)得出原函数.t2e.2t sint，因为sin(t .)=. sint。例5-3试求出下面函数的Laplace反变换。.17x5 . 7x4 2x3 x2 . x 1 G(x)=x6 11x5 48x4 106x3 125x2 75x 17 解从原来给出的问题，似乎用下面的语句就能直接求解出所需结果
得出的结果可读性很差。事实上，这样的问题因为方程x6 11x5 48x4 106x3 125x2 75x 17=0不存在解析解，所以导致原问题不存在解析解。若利用MATLAB的变精度算法vpa(f)则可以得出高精度的数值解，这里因篇幅所限将其解表示为y(t)=.556.2565e.3.2617t 1.7589e.1.0778tcos0.6021t 10.9942e.1.0778tsin0.6021t 0.2126e.0.5209t 537.2850e.2.5309tcos0.3998t. 698.2462e.2.5309tsin0.3998t[]例5-4对例5-1给出的原函数f(t)，试得出Ld5f(t)/dt5和s5L[f(t)]之间的关系。解如果想求解这样的问题，可以利用符号运算工具箱中的diff()函数对函数f(t)求五阶导数，再进行Laplace变换，则
对f(t)进行Laplace变换，并将变换结果乘以s5 ，将得出的结果与前面直接得出的结果相减，可以得到差为6s. 48。
由于二者之差不为0，所以看起来和式（5-1-3）是不同的。这是因为f(t)函数有f(0)=f′(0)=0，但其高阶导数在t=0处的值不为0，故不满足式（5-1-3），而满足式（5-1-2）。由式（5-1-2）直接可见，考虑初始条件后，得出的结果和上述的差完全一致。· 152·高等应用数学问题的MATLAB求解（第四版）
例5-5试推导出L[d2f(t)/dt2]的微分公式。解MATLAB的符号运算工具箱还可以进行一些简单的Laplace变换公式推导。假设想导出f(t)的二阶导数的Laplace变换，首先应该先定义一下f(t)函数，这可以通过如下语句实现，并推导出二阶导数的Laplace变换公式
得出的结果为s.2*laplace(f(t),t,s)-s*f(0)-D(f)(0)，其数学表示为s2F(s). sf(0). f′(0)，展开即可发现，该结果与式（5-1-2）中的式子完全一致。当然，该功能可以进一步引申，求出函数八阶导数的Laplace变换
例5-6已知f(t)=e.5t cos(2t 1) 5，试求出L[d5f(t)/dt5]。解这个例子是上个例子的引申。若已知某具体函数f(t)，则可以将diff()函数与laplace()函数结合起来使用，这样用下面的MATLAB命令则可以得出所需的结果
其结果为F = (1475 cos 1s . 1189 cos 1 . 24360 sin 1 . 4282sin1s)/(s2 10s 29) 。对化简后的结果其实还可以采用其他化简方法微调。例如，想将分子多项式合并同类项，则可以给出如下语句
则将得出的显示结果为(1475 cos 1 . 4282 sin 1) s . 1189 cos 1 . 24360 sin 1 F1(s)=s2 10s 29 5.1.3 Laplace变换问题的数值求解前面给出了Laplace变换的求解函数laplace()，该函数可以推导出很多时域函数的Laplace变换的解析解，同时也有很多函数Laplace变换的解析解不存在或不适合用解析解方法求解，所以应该考虑数值方法求解Laplace变换问题。JurajValsa开发了基于数值方法的Laplace反变换的MATLAB函数INVLAP()[1, 2]，该函数的调用格式为其中，原函数由含有字符s的字符
串表示，(t0,tn)为感兴趣的区间且t0.=0，N为用户选择的计算点数，用户可以选择不同的N值来检验运算的结果。“其他参数”的选取可以参考原函数的联机帮助，不过这里建议：除非特别需要没有必要去人为修改这些默认参数。对源程序代码进行扩展，可以建立新的函数清单如下

该函数的调用格式如下
该函数支持多种调用格式，其中，G为Laplace变换表达式的字符串，如果同时提供了H，则G为前向通路的传递函数模型字符串，H为负反馈回路传递函数的字符串；如果需要描述输入信号，则u可以为输入信号的Laplace变换字符串，或输入时域信号的匿名函数句柄；输入信号还可以由采样点(tx,ux)表示；如果只考虑G模型的响应，可以将H设置为0。例5-7试用数值方法重新求解例5-3中给出函数的Laplace反变换的问题。解由前面给出的例题可知，虽然原问题的解析解是未知的，可以通过符号运算工具箱相关函数求出高精度的数值解。对同样的问题，可以利用MATLAB语句将其变换成关于s的字符串变量G，对其进行数值求解，则可以得出该函数Laplace反变换的数值解。和精确的数值解相比，这个例子的相对误差为1.83×10.5％，应该满足一般科学运算的要求。
从计算量看，如果将计算点数从100增加到5000，采用INVLAP_new()函数所需时间为0.61s，而采用ilaplace()与subs()函数则需时261s，可见对某些问题来说，采用数值方法更高效。在一般应用中，如果某复杂系统G(s)输入信号的Laplace变换可求且已知为R(s)，可以得出输出信号的Laplace变换表达式为Y(s)=G(s)U(s)，则可以通过前面给出的方法得出输出信号的数值解。(s0.4 0.4s0.2 0.5)例5-8考虑已知的Laplace变换表达式G(s)=√s(s0.2 0.02s0.1 0.6)0.4(s0.3 0.5)0.6。试用数值方法绘制出t ∈ (0,1)区间内的Laplace反变换时域函数曲线。解该函数不能用ilaplace()函数求取Laplace反变换解析解，只能采用数值方法求解此问题。选择计算点数N=1000，则可以由下面的语句对G(s)进行Laplace反变换，得出的时域函数曲线如· 154·高等应用数学问题的MATLAB求解（第四版）图5-1所示。增加计算点数到N=5000仍然能得出一致的结果，说明这样得出的结果是正确的。
图5-1数值Laplace反变换时域函数曲线值得指出的是，分数阶多项式pΓ (x)本身的精确解应该对应于无穷级数，其Laplace反变换的解析解是不可能求出的，必须借助于数值方法来求解原始问题。另外，该数值算法速度足够快，实际求解时间只需0.3s。如果给出函数u(t)的Laplace变换的解析解是不存在的，则可以考虑对原始的问题用数值方法求解输入信号的Laplace变换。分析INVLAP()源程序可见，该函数主要采用循环结构，每一个循环步骤内生成一个s向量，根据该向量可以采用数值积分的方法求出输入信号的Laplace变换L[u(t)]=∫∞ u(t)e.stdt=U(s)（5-1-8）0 其中，s为向量。由于积分中含有e.st 项，在实际应用中只要求解的时间区间(0,tn)足够大，即可以用该有限时间内的积分代替无穷积分，得出近似的数值Laplace变换。如果输入信号只是已知一些样本点向量x0，u0，则实际u(t)信号在t时刻的值可以通过插值求出，插值方法在第8章中将详细叙述。例5-9假设前面给出的G(s)为某系统的分数阶传递函数模型，试求出该系统在u(t)=e.0.3tsin t2 激励下的响应曲线。解分数阶传递函数的输入方法与前面给出的完全一致。将输入信号用匿名函数描述出来，则可以由下面的语句计算出输出信号的曲线，如图5-2所示。该函数内部采用了数值积分运算，耗时比前面介绍的Laplace反变换数值算法长得多，大约需要9s。若想解决耗时问题需要在算法上有所突破。
现在重新考虑输入信号，假设已知在t ∈ (0,15)周期内已知的采样点可以直接如下求出，基于这些样本点可以求出输出信号的数值解，该解与前面得出的完全一致，但由于该函数内部采用了插值运算，所以耗时极长，计算点数缩减一半以后本例仍需要大约90s的时间。
图5-2复杂分数阶系统的输出曲线例5-10假设复杂开环无理模型如下[3]试绘制单位负反馈统的阶跃响应曲线。[]2s，inh(0.1√s)1 G(s)=0.1√s√ssinh(√s) 解开环无理传递函数可以由字符串表示，这样由下面语句直接绘制出系统的闭环阶跃响应曲线，如图5-3所示。
图5-3闭环系统的阶跃响应曲线 5.2 Fourier 变换及其反变换 5.2.1 Fourier 变换及反变换定义与性质 Fourier变换的一般定义为F[f(t)]=∫∞.∞f(t)e.jσtdt=F(ω)如果已知Fourier变换式F(ω)，则可以由Fourier反变换公式反演出f(t)函数为f(t)=F.1[F(ω)]=12.∫∞.∞F(ω)ejσtdω这里仍将不加证明地列出一些Fourier变换的性质。（1）线性性质。若a与b均为标量，则F[af(t)±bg(t)]=aF[f(t)]±bF[g(t)]。（5-2-1）（5-2-2）
· 156·高等应用数学问题的MATLAB求解（第四版）（2）平移性质。F[f(t± a)]=e±jaσF(ω)。
（3）复域平移性质。F[e±jatf(t)]=F(ω. a)。
（4）微分性质。Fdf(t)/dt=jωF(ω)，一般地，n阶微分可以由下式求出

dn F dtn f(t)=(jω)nF[f(t)]（5-2-3）t （5）积分性质。F[[∫.∞ f(τ] )dτ[] = F j(ωω] ) ，一般地，函数f(t)的n重积分的Fourier变换可以由下式求出F [∫t (··· )∫t f(τ)dτn] = F(j[ωf)(n t)]（5-2-4）.∞ .∞ （6）尺度变换。F[f(at)]=Fω/a/a。
（7）卷积性质。F[f(t). g(t)]=F[f(t)]F[g(t)]，其中，卷积定义仍为式（5-1-5）。

5.2.2 Fourier 变换的计算机求解和Laplace变换一样，应该先声明符号变量，并定义出原函数为f，这样就可以按如下的格式调用Fourier变换求解函数fourier()，得出该函数的Fourier变换式
给出了Fourier变换表达式，则可以通过ifourier()函数求解该函数的Fourier反变换问题。该函数的具体调用格式为
从上面的语句可以看出，如果定义了已知函数，则可以用一个语句求解出其Fourier变换或反变换式，变换函数的调用和Laplace变换一样简单。如果已知的是MATLAB得出的Fourier变换式，则可以直接使用ifourier()函数进行变换。例5-11考虑f(t)=1/(t2 a2)，a>0，试写出该函数的Fourier变换式。解可以用下面的语句得出原函数的Fourier变换，得出F=.e.a|σ|/a。
例5-12假设时域函数为f(t)=sin2(at)/t，a>0，试求出其Fourier变换。解由给定的式子，可以用下面的语句获得原函数的Fourier变换
得出的结果为..jheaviside(.2a . ω)/2. .jheaviside(2a. ω)/2 .jheaviside(.ω)，其中，heav-iside(x)函数为x的阶跃函数，又称为Heaviside函数，当x>0时，该函数的值为1，x=0取0.5，否则为0。当ω>2a，则三个heaviside()函数的值均为1，故F(ω)=0。若ω..2a，则三个函数的值均为0，故F(ω)=0。若0<ω<2a，则第二和第三个heaviside()的值为1，故F(ω)=.j./2。当0>ω>.2a，则F(ω)=j./2。综上所述，原函数的Fourier变换可以手工化简成{F[f(t)]=0, |ω| > 2a .j.sign(ω)/2,|ω| < 2a √ 例5-13假设函数为f(t)=e.a|t|/ |t|，试求解Fourier变换问题。解先考虑用现成的函数fourier()对给出的原函数进行变换
得出的结果为√2./(2 √ |w . ja|) √2./(2 √ |w ja|)。5.2.3 Fourier 正弦和余弦变换Fourier正弦变换的一般定义为∫∞ FS[f(t)]=f(t)sin(ωt)dt=Fs(ω)（5-2-5）0 Fourier余弦变换的一般定义为∫FC[f(t)]=∞ f(t)cos(ωt)dt=Fc(ω)（5-2-6）0 MATLAB语言的符号运算工具箱中并未直接提供余弦Fourier变换的函数，所以可以考虑采用符号积分的方法直接求取余弦Fourier变换。下面将提供具体例子演示正弦和余弦Fourier变换的推导方法。例5-14试求出f(t)=tne.at ，a>0，n=1,2,··· ,8的余弦Fourier变换。解解决这样的问题可以采用循环结构，对不同的i值，可以用直接积分的算法求取Fourier余弦变换，并将得出的结果进行化简，如表5-1所示。
其实，按照数学手册[4] 中给出的公式，对整数n，可以得出ne.at( 1 )n 1( a )n 1[n/2](.1)mC2m 1(ω[] ∑)2m 1FC t= n! Re a jω = n! a2 ω2 n 1a （5-2-7）m=0表5-1n取不同值时Fourier余弦变换结果n值FC[f(t)]a2 . ω2 , (.a2 3ω2)a, 6(.a2 2aω ω2)(.a2 . 2aω ω2) , 24 a(a4 . 10a2ω2 5ω4)1～4 (a2 ω2)2 . 2(a2 ω2)3 (a2 ω2)4 (a2 ω2)5 5,6 .120(.a ω)(a ω)(a2 (a. 2 4 ωa ω2 )6 ω2)(a2 4ωa ω2) , . 720(.a6 21a4(aω22 . ω352ω)74a2 7ω6)a 7 5040(a4 4a3ω . 6a2ω2 . 4aω3 ω4)(a4 . 4a3ω . 6a2ω2 4aω3 ω4)(a2 ω2)8
8 40320 a (.a2 3ω2)(.a6 33a4ω2 . 27a2ω4 3ω6)(a2 ω2)9
{例5-15试求取分段函数f(t)=cos t, 0 <x<a 的Fourier余弦变换。0,其他解研究Fourier余弦变换定义可见，式（5-2-6）的被积函数在t∈ (a,∞)区间的值为0，这样，其积分亦为0，故整个积分问题就变成t∈ (0,a)区间的积分问题了，故可以用下面的语句求出该函数的Fourier 余弦变换>> syms t w; syms a positive; f=cos(t); F=simplify(int(f*cos(w*t),t,0,a))
· 158·高等应用数学问题的MATLAB求解（第四版）其结果是分段函数. .a/2 sin2a/4,F(ω)=sin[a(ω. 1)]/[2(ω . 1)] sin[a(ω 1)]/[2(ω 1)],ω ∈ {.1, 1}
. ω .∈ {.1, 1}0,a=./2&ω=3如果将f用分段函数直接表示，也可以得到一致的结果。
5.2.4 离散Fourier正弦、余弦变换离散Fourier正弦、余弦变换又称为有限Fourier正弦、余弦变换，和前面介绍的Fourier正弦、余弦变换相比，其积分区间从t∈ (0, ∞) 变成了t ∈ (0,a)，故其定义为∫ ∫Fs(k)=a f(t)sink.tdt,Fc(k)=a f(t)cosk.tdt（5-2-8）0a 0 a 相应地，可以定义出有限Fourier正弦、余弦反变换为∑ ∑ f(t)=2 ∞Fs(k)sink.t,f(t)=1 Fc(0) 2 ∞Fc(k)cosk.t（5-2-9）a aaa a k=1k=1和前面定义的不同，反变换不再是积分式子，而是无穷级数的求和。下面通过例子介绍如何用MATLAB及其符号运算工具箱求取给定函数的有限变换。{例5-16考虑分段函数f(t)=t,t.a/2其中，a>0，试求其离散Fourier正弦变换。a . t,t>a/2，解函数的离散Fourier正弦变换可以由下面的语句直接求出
( )得出的结果为(.1)k/2a2((.1)k 1/2. j)/(k2.2)。如果采用分段函数的方式描述f函数，也可以得出完全一致的结果。
5.2.5 快速Fourier变换从前面的叙述看只有一部分简单的函数可以通过各种Fourier变换的公式得出其相应的变换表达式，限制了这样解析变换方式在实际中的应用。在实际应用中，更多地采用数值形式直接对信号的采样值进行离散Fourier变换。离散数据xi,i=1,2,··· ,N的Fourier变换是数字信号处理的基础。离散Fourier变换的数学表示为N∑ X(k)=xie.2.j(k.1)(i.1)/N ,其中，1.k.N（5-2-10）i=1N其逆变换定义为x(k)=N 1 ∑ X(i)e2.j(k.1)(i.1)/N ,其中，1.k.N（5-2-11）i=1快速Fourier变换（fastFouriertransform，FFT）技术是求解离散Fourier变换的实用、也是通用的方法。MATLAB提供了内核函数fft()，其调用格式很简单，为
可以直
· 160·高等应用数学问题的MATLAB求解（第四版）相应地，Mellin反变换可以定义为1 ∫c j∞ f(x)=M.1[M(z)]=2.jc.j∞ M(z)x.zdz（5-3-2）MATLAB符号运算工具箱并没有直接可用的Mellin正反变换的函数，但仍可以通过积分的形式计算Mellin变换。下面将通过例子解决这类问题。例5-18考虑时域函数f(t)=lnt/(t a)，其中，a>0，试求其Mellin变换。解根据定义可以用下面语句求取该函数的Mellin变换（注意：新版本无法得出所需的解）。
经过化简，可以立即得出结果M[f(t)]=az.1.(lnasin.z. .cos.z)csc2 .z。例5-19假设给定函数f(t)=1/(t a)n ,(a>0)，对若干个n值求取Mellin变换，并总结出对一般n值的规律。解下面的语句将给出n=1,2,··· ,8时Mellin变换的结果
上面的循环结构可以得出如下结果az.1. csc .z .a.2 z.(z. 1) csc .z 1/2a.3 z.(.2 z)(z. 1) csc .z .1/6az.4.(z. 1)(.2 z)(.3 z) csc .z 1/24a.5 z.(z. 4)(.3 z)(.2 z)(z. 1) csc .z .1/120a.6 z.(.5 z)(z. 4)(.3 z)(.2 z)(z. 1) csc .z 1/720a.7 z.(.6 z)(.5 z)(z. 4)(.3 z)(.2 z)(z. 1) csc .z .1/5040a.8 z.(.7 z)(.6 z)(.5 z)(z. 4)(.3 z)(.2 z)(z. 1) csc .z 所以，可以总结出一般的Mellin变换规律为n.1M [ 1 ] =(.1)k.1. a z.n ∏ (z . i) csc .z (t a)n (n . 1)! i=1和很多其他积分变换类似，绝大多数函数是不存在Mellin变换解析解的，所以可以考虑通过数值积分引入数值Mellin变换的概念。可以编写mellin_trans()来计算函数的数值Mellin变换，调用格式为其中，“属性对”与integral()函数要求
的完全一致。
例5-20试求出函数f(x)=sin(3x0.8)/(x 2)1.5的数值Mellin变换。解无论用什么工具都不能得出该函数Mellin变换的解析解，所以数值Mellin变换是的选择。可以用匿名函数先描述原函数，然后调用数值Mellin变换函数直接求解，结果如图5-5所示。
图5-5数值Mellin变换5.3.2 Hankel 变换及求解Hankel变换是另一类常用的数学变换，ν阶Hankel变换的数学表达式为∫H[f(t)]=∞ tf(t)Jγ(ωt)dt=Hγ(ω)（5-3-3）0 其中，Jγ(z)为ν阶Bessel函数，在MATLAB下可以用J还可以定义ν阶Hankel反变换公式为H .1[H(ω)]=∫∞ ωHγ (ω)Jγ(ωt)dω（5-3-4）0 由定义可见Henkel变换实际上是无穷积分问题，所以可以通过int()函数直接求解。22例5-21求取f(t)=e.at/2函数的零阶Hankel变换。解利用MATLAB的积分还是可以直接求取原函数的Hankel变换
可以得出变换结果为e.σ2/(2a 2)/a2 ，由反变换语句可以还原出f1(t)=f(t)。用积分方法只能求出极少函数的Hankel变换，并只能求取低阶Bessel变换，例如ν=0，而对绝大多数的Hankel变换无能为力，所以应该考虑数值Hankel变换方法。可以编写一个MATLAB函数来计算数值Hankel变换，
其中，“属性对”与integral()函数描述的一致。
例5-22考虑前面的例子，若取a=2，试绘制不同阶次的Hankel变换曲线。解下面的语句可以直接绘制出各阶Hankel变换的数值曲线，如图5-6所示，前面已经得出了ν=0阶变换的理论解，可以直接对其取值，也叠印在新图形上。可见，零阶Hankel变换的结果与理论值完全一致。此外，由得出的结果看，高阶Hankel变换的衰减很慢，不能由得出的结果直接求解数值· 162·高等应用数学问题的MATLAB求解（第四版）Hankel反变换，必须考虑足够大的ω范围才行。如果函数F(ω)已知，则可以由类似的MATLAB语句计算数值Hankel反变换。
ν=0，与理论值一致ν =1 ν =2 图5-6不同阶次的数值Hankel变换曲线5.4 z 变换及其反变换严格说来，z变换并不属于积分变换，但由于其定义、性质和求解方法也类似于Laplace变换，并且该方法可以在描述序列信号中起重要作用，所以本节将介绍z变换及其求解方法，并将给出通过长除法求取z反变换的数值方法的MATLAB实现。5.4.1 z 变换及反变换定义与性质离散序列信号f(k)，k=1,2,··· 的z 变换可以定义为∑∞Z[f(k)]=f(k)z.k =F(z)（5-4-1）k=0类似于前面介绍的Laplace变换和Fourier变换，z变换也有很多类似的性质，这里仍不加证明地列出一些性质如下：（1）线性性质。若a与b均为标量，则Z[af(k)± bg(k)]=aZ[f(k)]± bZ[g(k)]。
（2）后向平移性质。Z[f(k. n)]=z.nF(z)。
（3）前向平移性质。对非零初值的问题，前向平移的z变换可以由下式计算

n.1∑ Z[f(k n)]=znF(z). z n.if(i)（5-4-2）i=0特别地，对零初值问题有Z[f(k n)]=znF(z)。（4）s-域比例性质。Z[r.kf(k)]=F(rz)。（5）频域微分性质。Z[kf(k)]=.zdF(z)/dz。
（6）频域积分性质。Z[f(k)/k]=∞ Fω (ω)dω。
（7）初值性质。klim 0 f(k)=limF(z)。
→z→∞
（8）终值性质。如果F(z)无单位圆外的极点，则limf(k)=lim(z. 1)F(z)。

∫z k→∞ z→1（9）卷积性质。Z[f(k). g(k)]=Z[f(k)]Z[g(k)]，式中，离散信号的卷积算子. 定义为∑∞f(k). g(k)=f(k)g(k. l)（5-4-3）l=0给定z变换式子F(z)，则其z反变换的数学表示为∮ f(k)=Z.1[f(k)]=21 .j F(z)zk.1dz（5-4-4）5.4.2 z 变换的计算机求解利用MATLAB的符号运算工具箱，则z变换及其反变换可以很容易地求取出来，掌握这样的工具可以免除复杂问题的手工推导，既节省时间也能避免底层的低级错误。利用符号运算工具箱中提供的ztrans()和iztrans()函数可以得出给定函数的正反z变换。这两个函数的调用格式为
若原函数只有一个变量，则调用时无须给出k和z。例5-23求解f(kT)=akT. 2 (akT 2)e.akT 函数的z变换问题。解原函数的z变换可以用下面的语句来完成
该结果可以表示为aTz2zaTze.aT ( z ).1 Z[f(kT)]= 2zeaT (z . 1)2 . z . 1(z . e.aT )2 e.aT . 1例5-24考虑F(z)=q/(z.1 . p)m 函数的z反变换问题，这里可以对不同的m值进行反变换，并总结出一般规律。解根据要求，可以用符号运算工具箱求出m=1,2,··· , 8 的z 反变换
对不同的i 值循环可以得出如下结果.q/p(1/p)n ，q/p2(1 n)(1/p)n ， .1/2q(1/p)n(1 n)(2 n)/p3 1/6q(1/p)n(3 n)(2 n)(1 n)/p4 ， .1/24q(1/p)n(4 n)(3 n)(2 n)(1 n)/p5 1/120q(1/p)n(5 n)(4 n)(3 n)(2 n)(1 n)/p6 .1/720q(1/p)n(6 n)(5 n)(4 n)(3 n)(2 n)(1 n)/p71/5040q(1/p)n(7 n)(6 n)(5 n)(4 n)(3 n)(2 n)(1 n)/p8总结上述结果的规律，可以写出一般的z 反变换结果为
[]Z .1 =(.1)mq m∏.1 (n i)(z.1 .qp)m (m . 1)!pn mi=1· 164·高等应用数学问题的MATLAB求解（第四版）5.4.3 双边z 变换前面叙述的z变换由于描述的是n.0时序列信号的性质，所以又称为单边z变换，如果将n的范围扩展到整个整数集合，则可以定义出双边z变换∑∞Z[f(k)]=f(k)z.k =F(z)（5-4-5）k=.∞ MATLAB中并未提供双边z变换的求解函数，但可以从定义直接求出给定函数f的双边z变换表达式FMAT
LAB的symsum()函数有时不支持(.∞, ∞) 区间的求和。例5-25试求出分段函数f(n)的双边z变换[5]，其中，f(n)={ 2n ,n . 0 .3n ,n<0。解采用前面介绍的方法，整个求和可以分成两个区间单独计算，可以由下面的语句直接得出函数的双边z变换。下面语句将得出F=z/(z. 2) z/(z. 3)。
5.4.4 有理函数z 反变换的数值求解很多函数z反变换的解析解是不能由iztrans()求出的，即使对某些有理函数也不能解析地求出其z反变换。假设一般有理函数z变换表达式可以写成()b0 b1z.1 b2z.2 bm.1z.(m.1) bmz.m Fz.1=z.d ··· （5-4-6）a0 a1z.1 a2z.2 an.1z.(n.1) anz.n ··· 将该函数作关于z.k 的幂级数展开则可以写出∑() ∞Fz.1=f0 f1z.1 f2z.2 =fkz.k （5-4-7）··· k=0

)上式恰巧是z变换的定义式。可以通过长除法展开Fz.1，从而得出Fz.1函数z反变换的数值解。可以编写长除法函数，其调用格式为其中，d为纯延迟的步数，系数向量num=[b0,b1,b2,··· ,bm]，den=[a0,a1,a2,··· ,an]，默认的计算点数为N=10。应用数值方法求出的序列信号由y向量返回。
例5-26试用数值方法求G(z)=z2 0.4的z 反变换。z5 .4.1z4 6.71z3 .5.481z2 2.2356z.0.3645解原函数分子和分母同时乘以z.5 则可以得出下面所需的表达式F (z.1) = z.3 1 0.4z.2 1 . 4.1z.1 6.71z.2 . 5.481z.3 2.2356z.4 . 0.3645z.5 这样由下面的语句可以直接得出原函数的z反变换，得出的序列如图5-7所示。
图5-7z反变换的数值解5.5 复变函数问题的计算机求解5.5.1 复数矩阵及其变换前面介绍过，MATLAB可以用来直接表示复数矩阵。假设已知一个复数矩阵Z，则可以使用简单函数对该矩阵进行如下变换：（1）共轭复数矩阵
（2）实部、虚部提取
（3）幅值、相位表示

其中，相位的单位为弧度。
例5-27重新考虑例4-41中的Jordan标准型问题。由于原矩阵存在复数特征值，所以有时需要用手工的方法修改变换矩阵。这里采用下面语句修改变换矩阵，也能起到同样的作用。
5.5.2 复变函数的映射若某函数f(z)的自变量z为复数，则该函数称为复变函数。由于复数是MATLAB的基本的数据结构，在大多数算法中均未特别区分实数和复数，所以前面介绍的绝大多数内容均可以直接用于复变函数的分析。例如，前面的微积分运算的解析解和数值解均可以直接用于复变函数的微积分运算。例5-28已知某复变函数为f(z)=z2 3z 4 ，其中，z为复数变量，试求出f(3)(.j√5)的值。(z . 1)5 解由下面语句可以立即得出所需的结果为d3=0.8150. j0.6646。
在复变函数中一种很重要的数学变换形式是映射，即函数可以从一个变量z变换成另一个变量w的函数，其中，z=g(w)为给定的函数。经常使用的映射是平移映射z=w γ，反演映射z=1/w和双线性映射z=(aw b)/(cw d)，这里，γ为给定复数，a，b，c，d为给定实数。其中，平移映射将函数的原点平移到γ点；反演映射可以将单位圆内外的点相互映射，而双线性变· 166·高等应用数学问题的MATLAB求解（第四版）换实现直线与圆的相互映射。求解函数映射的直接的MATLAB函数是subs()，下面通过例子演示函数的映射。例5-29考虑例5-28中的复变函数f(z)=z2 3z 4 ，试得出z = s . 1 下的映射函数F(s)。(z . 1)5 s 1 解映射问题由subs()函数直接得出，得出的结果需要化简
该语句可以直接得出映射函数为F(s)=.(s 1)3(4s 2 3s 1)/16。例5-30还可以绘制出s=(z. 1)/(z 1)的映射效果，z平面单位圆内如图5-8（a）所示的点可以通过映射，映射成s左半平面的点，如图5-8（b）所示。
（a）z域散点与单位圆（b）s域映射效果图5-8z域到s域的映射示意图5.5.3 Riemann面绘制复变函数映射的三维图形表示和实函数是不一致的，复变函数映射应该首先采用函数cplxgrid()生成极坐标网格，用户可以根据给定的单值复变函数公式计算出变量f，然后由cplxmap()函数绘制该复变函数的映射曲面，这类曲面又称为Riemann面。这些函数具体调用格式为
例5-31试绘制出复变函数f(z)=z3 sin z2 的映射曲面。解由下面的语句可以立即绘制出相应的复变函数的映射曲面，如图5-9所示。
对于复数变量z，多值复变函数的Riemann面可能有多个分支，这些分支又称为Riemannn叶（Riemannsheet），比如f(z)=√z 就有n个分支。MATLAB提供了绘制其各次方根Rie-nmann曲面的函数
可以直接绘制出√z 的Riemann曲面。34例5-32试绘制√z 和√z 的Riemann曲面。34解由cplxroot()函数可以直接绘制出√z 和√z 的Riemann面，如图5-10（a）、（b）所示。
从上面给出的cplxroot()函数可见，其局限性在于只能绘制出方根函数的Riemann面，图5-9复变函数的Riemann面（a）三次方根（b）四次方根n图5-10.z 方根函数的Riemann面对其他类型的多值复变函数无能为力，所以可以考虑扩展cplxmap()函数：将该函数另存为cplxmap1()函数，再删除掉其中的mesh()函数和hold语句，这样就可以考虑多值复变函数Riemann面的绘制了。3例5-33利用上述思路重新绘制√z 的Riemann面。33解首先考虑重新绘制√z 函数的Riemann面。如果一个函数f1(z)是f(z)=√z 的一个分支，则另两3个分支可以由f1(z)e.2j./3和f1(z)e.4j./3求出。这样可以由下面语句直接绘制√z 的Riemann面，与图5-10（a）给出的完全一致。
5.6 复变函数问题的求解本节将介绍复变函数中一些常见问题的求解方法，首先将给出奇点行留数的概念与计算，然后计算有理函数的部分分式展开与Laurent级数展开问题，后介绍复变函数的封闭曲线积分问题的求解。5.6.1 留数的概念与计算在介绍留数概念之前，应该先介绍一下复变函数解析的概念。若函数f(z)在复平面的区域内各点处均为单值，且其导数为有限值，则称f(z)在复平面内为解析的，这样就将单值函数上不解析的点称为奇点。使得f(z)分母多项式等于零的奇点又称为极点。· 168·高等应用数学问题的MATLAB求解（第四版）假设z=a为f(z)函数上的奇点，若存在一个小整数m使得乘积(z. a)mf(z)在z=a点处解析，则称z=a为m重奇点。特别地，可以由[p,m]=poles(f)函数计算出函数的极点，如果极点个数大于1，则函数的极点由列向量p返回，极点的重数由向量m返回。如果想得出某个区域(a,b)内的极点，则可以调用下面的语句
若z=a为f(z)函数的单奇点，则函数在该奇点处的留数（residue）可以定义为Res[f(z),a] =lim(z. a)f(z)（5-6-1）za→若z=a为函数f(z)的m重奇点，则该点的留数定义为Res[f(z),a]=zlim a (m . 1 1)! ddzmm..11 [f(z)(z. a)m] （5-6-2）→利用MATLAB的符号运算工具箱求留数的方法很简单。假设已知奇点a和重数m，则用下面的MATLAB语句自然可以求出相应的留数。
()例5-34试求出函数f(z)=z3(z 1 . 1) sin z . 3 e.2z 的留数。解对原函数的分析可见，z=0是三重奇点，z=1是单奇点，故可以直接使用下面的MATLAB语句将这两个奇点处的留数分别求出。可以得出z=0处的留数为F1=.√3/4 1/2，z=1处的留数为F2=e.2 sin 1/2 √3 e.2 cos1/2。
基于上述考虑，可以编写出一个同时求解极点、重数与留数的函数residuesym()，其调用格式为其中，描述极点感兴趣区间的a 和b 可以略去。

例5-35试求函数f(z)=(sinz. z)/z6 的留数。解乍看该函数很容易认定z=0为6重奇点，所以用下面的语句很容易就可以求出该点处的留数值，其值为1/120，而极点的重数为m=3。
[]不严格地说，从k=1开始尝试，能够使得zlim a dk.1(z . a)k f(z)/dzk.1 < ∞ 成立的小的k →就是奇点的重数，记作m，可以考虑k=2，可见导数F1为无穷大，所以再试验更大的k值。对此例子来说，k的小值为m=3，留数的值F2仍然为1/120。试凑更大的k值，如k=20，也不会改变求出的留数值，F3=1/120。
可见，若选择的n值大于或等于奇点的实际重数，则可以正确得到该函数的留数。在一般应用时可选择一个较大的n值来求取留数。其实考虑对sinz进行Taylor幂级数展开，则可以看出，z=0极点的实际重数是3而不是6。f(z)=(z.z3/6 z5/120.z7/5040 ··· ).z = .1/6 z3/120.z5/5040 ··· z6 z3 ()例5-36试求出函数f(z)=1/zsinz的留数。解分析该函数，因为sinz在z=0点的收敛速度和z是一样的，显然，z=0点为f(z)的二重奇点，这时，相应的留数可以用下面语句求出c0=0。
进一步分析给定函数f(z)，可以发现该函数在z=± k.处均不解析，其中，k为正整数，且这些点是原函数的单奇点，由于MATLAB的符号运算工具箱并未给出整数的定义，所以这里只能对一些k值进行试探，求出它们的留数，后将结果归纳成所需的公式。
对向量k=[.4, 4, .3, 3, .2, 2, .1,1]，可以得出留数为c=[.1/(4.),1/(4.),1/(3.)，.1/(3.)，.1/(2.)，1/(2.)，1/.，.1/.]。综上，可以归纳出Res[f(z),±k.]=±(.1)k/(k.)。事实上，如果利用MATLAB新版的整型变量定义，可以直接得出上述的结果。
值得指出的是，这个问题不适合使用residuesym()函数，因为其分母不是多项式，不能得出所有的极点，所以必须由留数的定义直接求解。5.6.2 有理函数的部分分式展开考虑有理函数G(x)=B(x)= b1xm b2xm.1 ··· bmx bm 1（5-6-3）A(x)xn a1xn.1 a2xn.2 an.1x an··· 其中，ai和bi均为常数。有理函数的互质概念是一个非常重要的概念。所谓互质，就是指多项式A(x)和B(x)没有公约数。对一般给定的两个多项式来说，用手工方式判定多项式互质还是比较困难的，但利用MATLAB符号运算工具箱中的gcd()函数可以直接求出两个多项式的公约数。该函数的调用方法为
其中，A和B分别表示两个多项式，该函数将得出这两个多项式的公约数C，若得出的C为多项式，则两个多项式为非互质的多项式，这时两个多项式可以约简为A/C和B/C。例5-37给出两个多项式A(x)=x4 7x3 13x2 19x 20,B(x)=x7 16x6 103x5 346x4 655x3 700x2 393x 90，试判定它们是否互质。解求解这样的问题可以采用MATLAB语言提供的gcd()函数完成· 170·高等应用数学问题的MATLAB求解（第四版）
可见，两个多项式具有公约数d=x 5，故两个多项式不是互质的，这两个多项式可以进一步简化为(x3 2x2 3x 4)/[(x 2)(x 3)2(x 1)3]。
若互质多项式A(x)=0的根均为相异的值.pi,i=1,2,··· ,n，则可以将G(x)函数写成下面的部分分式展开形式G(x)=x r1 p1 x r2 p2 ··· x rn pn （5-6-4）其中，ri为留数，ri=Res[G(x),.pi]，其值可以由下面的极限式求出ri=Res[G(x),.pi]=lim(x pi)G(s)（5-6-5）x→.pi 如果分母多项式中含有(x pi)k 项，亦即.pi为k重根，则相对这部分特征值的部分分式展开项可以写成ri ri 1 ri k.1 （5-6-6）x pi (x pi)2 ···(x pi)k 这时，ri j.1 可以用下面的公式直接求出，其中一次项系数为留数，其他为普通系数。ri j.1 =(j . 1 1)! lim ddxjj..11 [(x pi)k G(x)],j =1, 2, ··· ,k（5-6-7）x→.pi MATLAB语言中给出了现成的数值函数residue()求取有理函数G(x)的部分分式展开表示，其中，a=[1,a1,a2,··· ,an]，
b=[b1,b2,··· ,bm]，返回的r和p向量为式（5-6-4）的ri，pi系数，若有重根则应该相应地由式（5-6-6）中给出的系数取代。k为余项，对m<n的函数来说该项为空矩阵。该函数并未给出.pi 是否为重根的自动判定功能，所以部分分式展开的结果需要手动写出。值得指出的是，数值运算对含有重奇点的问题经常会导致不精确的结果。例5-38试求下面有理函数的部分分式展开s3 2s2 3s 4 G(s)=s6 11s5 48s4 106s3 125s2 75s 18 解用下面的语句可以求出该函数的部分分式展开
其中，p为奇点向量，n，d1为每个p值对应系数的分子和分母数值。由数值方法直接求出的分母多项式的根为小数，有一些误差。事实上，该分母多项式的特征值为：.3 为二重奇点，.2 为单奇点，.1 为三重奇点。分析奇点的情况，可以写出部分分式展开为1772 11 G(s)=.8(s 3) . 4(s 3)2 s 2 8(s 11) . 2(s 1)2 2(s 1)3 例5-39写出下面式子的部分分式展开。2s7 2s3 8 G(s)=s8 30s7 386s6 2772s5 12093s4 32598s3 52520s2 45600s 16000 解采用MATLAB自带的residue()函数，只能求解得出数值解，对本例给出的问题，可以用下面的语句直接求出有理函数的部分分式展开式
从得出的结果较难判定重根情况，故难以准确写出部分分式展开式。由得出的数据，假定p1=.5为三重实根，p2=.4为三重实根，p3=.2，p4=.1 为单个实根，可以写出部分分式展开的表达式为49995.903093068628488.583258044113040.999976250750473.1527861460G1(s)=(s 5) (s 5)2 (s 5)3 .(s 4)21449.55550223475481.33332013621.22222222240.0023148148 (s 4)2 . (s 4)3 (s 2) (s 1) 应该指出，上面给出的展开方式在分母上作了近似，实际数值方法得出的分母不精确。另外，MATLAB数值运算在处理重根问题中经常会导致不精确的结果。所以对这样的问题来说，可以考虑编写更好的解析算法。新版的MATLAB符号运算工具箱提供了部分分式展开的函数partfrac()，其调用格式为
例5-40考虑例5-38中给出的函数f(s)，试用解析方式求出其部分分式展开。解用下面的语句可立即得出该函数的部分分式展开式，该结果与原例中数值结果完全一致。
得出的结果为1772 11 G1(s)=.8(s 3) . 4(s 3)2 s 2 8(s 11) . 2(s 1)2 2(s 1)3 例5-41仍考虑例5-39中给出的有理函数G(s)，试用解析方法写出其部分分式展开。解原函数可以由residue()函数进行部分分式展开。若将得出的部分分式展开减去原函数并化简，则结果为0，表示得出的结果是正确的。
这样，可以得出原分式的部分分式展开如下，经检验该展开式与原式完全一致。13041 341863 7198933 16444 193046 1349779 11 1 (s 5)3 12(s 5)2 144(s 5) . 3(s 4)3 9(s 4)2 . 27(s 4) 9(s 2) 432(s 1) 例5-42考虑例5-37中的非互质多项式构成的有理函数G(x)=A(x)/B(x)，试用数值方法和解析方法写出其部分分式展开。解用部分分式展开函数可以直接得出所需的展开，因为得出的公约数对应的展开项系数为0，会被直接忽略，所以无须事先求出公分母。
这样得出的解为A(x)71721 1 1 F(x)=B(x)= .4(x 3)2 . 8(x 3) (x 2) 2(x 1)3 . 2(x 1)2 8(x 1) 用前面介绍的解析解方法可以得出和x = .5奇点完全无关的解析解，亦即partfrac()函数同样适合于非互质有理函数的部分分式展开。· 172·高等应用数学问题的MATLAB求解（第四版）例5-43求有理函数的部分分式展开。.17x5 . 7x4 2x3 x2 . x 1 G(x)=x6 11x5 48x4 106x3 125x2 75x 17 解可以试图由partfrac()函数直接对原函数进行部分分式展开，然而因为原有理函数的分母不能进行因式分解，所以不能得出精确的部分分式展开表达式，该函数的调用也没有任何结果。后面将介绍近似的部分分数展开方法。
如果原有理函数分母D(x)的某无理数根为x0，通过任何的求根算法根本不能在有限位内得出x0的精确值，只能得到其近似值x.0，所以代入式（5-6-2）的留数公式则有1 dm.1 Res[f(x),x.0]=xlim x.0 (m . 1)!dxm.1 [(x . x.0)m f(.x0)]（5-6-8）→假设用vpa()函数能求出多项式方程所有的根xi,i=1,2,··· ,n，可以由这些根重组多项式的分母而取代原来的分母，则能得出新的函数f1(x)来取代式（5-6-8）中的f(x)，这样就能确保f1(x)和(x. x.0)m 在x.0处能真正相消，得出原函数的留数。在新版本下实际求取近似的部分分数展开系数时采用变量替换的方法而不能采用极限方法。该函数的清单为
该函数的调用格式为f其中，F 为有理函数的解析表达式，s 为自变量。返回的结果f 是部分分式展开的表达式。例5-44重新考虑例5-43中有理函数的部分分式展开问题。解由于原有理函数不能进行严格的因式分解，所以前面介绍的部分分数展开函数partfrac()对此问题无能为力，这里可以采用partfrac1()函数进行处理
得出的部分分式展开为（为排版需要这里只保留了有限几位有效数字）0.21255680.879464926 5.497076258j268.642522. 349.1231095jF(s)=x 0.5208596 x 1.07775887 0.602106591j x 2.53094582 0.399763j556.256530690.879464926. 5.497076258j268.642522 349.1231095j x 3.261731 x 1.07775887. 0.602106591j x 2.53094582. 0.399763j5.6.3 基于部分分式展开的Laplace反变换符号运算工具箱中提供的Laplace反变换函数ilaplace()可以较好地解决一般函数的Laplace反变换问题。但带有复特征值的有理函数的Laplace反变换问题不适合由该函数直接求解，这已经在例5-3中说明了，直接用该函数求解出的解可读性很差。观察例5-43的结果可以立即发现，若某项部分分式展开式为(a jb)/(s c jd)，则一定会含有其共轭项(a. jb)/(s c. jd)，这样就需要对下面的式子进行化简(a jb)e(c jd)t (a . jb)e(c.jd)t=αect sin(dt .)（5-6-9）其中，α = .2√a2 b2，且. = .arctan(b/a)。有了这样的算法，则可以用MATLAB语言编写出实现上述算法的数值函数pfrac()。该函数基于MATLAB的原始residue()函数，其内容为其中，p和K的定义和residue()函数

一致，r稍有不同，若相应的pi项为实数，则ri和residue()一致，若某个pi的值为复数，则ri，ri 1项分别为相应的α和.值。例5-45重新考虑例5-43中的问题，可以用MATLAB数值算法得出如下结果
这样，由得出的结果可以写出有理函数的可读性更好的Laplace反变换为L .1[F(s)]=.556.2565e.3.2617t. 881.0352e.2.5309tsin(0.3998t. 0.6559) 0.2126e.0.5209t. 11.1340e.1.0778tsin(0.6021t. 2.9829)比较此结果与例5-3中的结果，显而易见这个结果更简洁。5.6.4 Laurent级数展开第3章曾介绍过Taylor幂级数展开，可以将一个给定函数f(x)展开成(x. x0)的无穷级数形式，Laurent级数是Taylor级数的一个直接拓展。如果函数f(z)在圆环区域D：R1<<R2是解析的，且0.R1<R2< ∞，则可以将Laurent级数写成f(z)=|z . ∞z0| ck(z . z0)k （5-6-10）∑ k=.∞ 其中，系数可以如下直接计算ck =1 ∫ f(ζ)dζ,（5-6-11）2.j|z.z0|=θ (ζ . z0)k 1· 174·高等应用数学问题的MATLAB求解（第四版）式中，|z . z0| <ρ可以是满足R1<ρ<R2的任何圆周，如果函数f(z)在D区域解析，则Laurent级数是的。在实际应用中由式（5-6-11）得到级数系数是相当困难的，可以采用变通的方法来构造级数。假设原函数f(z)可以分解成两个子函数的积，记作f(z)=f1(z)f2(z)，其中，f2(z)适合于一般的Taylor级数展开，而另一个部分f1(z)可以展开成(z. z0)k 的级数，且k为负整数，则可以先作变量替换x=1/(z. z0)，将原函数替换成x的函数，即替换为z=(1 xz0)/x，则可以得出关于x的Taylor展开，记作F1(x)。这样，可以作变量替换x=1/(z. z0)将结果替换成z的函数。下面将通过例子演示这样的展开方法。例5-46试写出函数f(z)=z2e1/z 的Laurent级数展开。解由于函数中含有e1/z 项，所以在z=0处可能存在奇怪现象，不过这时z=0并不是极点，而是本征奇点，由原函数可以看出，可以将其分解为两个函数，其中，f1(z)=e1/z ，f2(z)=z2。由变量替换z=1/x，可以得出子函数f1(x)的Taylor级数展开，这样，可以利用x=1/z将变量替换回z的函数，再将f2(z)=z2 乘回f1(z)的Taylor级数展开，则可以得出Laurent级数。
得出的Laurent级数为1 z.1 z.2 z.3 z.4 z.5 F(z)=z2 z 2 6 24 120 720 5040 ··· 在一个相对较大的区间内，例如选z ∈ (.20,20)，可以绘制出原函数f(z)与有限项Laurent级数展开F(z)的曲线，如图5-11所示。可以看出，除在本征奇点z=0附近的一个很小的邻域内外，两条曲线几乎完全一致。
图5-11原函数与Laurent级数的曲线对比对给定的Laurent级数F(z)而言，因为本征奇点的存在，分母中存在无穷项z.k ，其留数定义为Laurent级数中z.1 项的系数c.1。例5-47现在考虑复数有理函数f(z)=1/(z. 1) 1/(z . 2j)，试写出其Laurent级数。解显而易见，f(z)函数在极点z=1与z=2j处不解析。这样，可以将复数平面分成三个区域：圆盘区域| z | <1；圆环区域1<| z | < 2 和圆环区域∞ > | z | >2。所以应该分三种情形考虑函数的Laurent级数展开问题：（1）若| z | <1，则意味着满足| z | < 2 或| z/2| <1，这时Taylor级数展开就足够了。可以利用Taylor级数的展开公式直接得到原函数的Taylor幂级数展开∑ 1 . 1 u = ∞u k ，若|u| <1，则展开式收敛k=0这样，原函数的Laurent级数（即Taylor级数）可以写成∑.1 .1/2j ∞(1 ) kF1(z)=1 . z 1 . z/(2j)= . 1 (2j)k 1zk=0（2）若1<| z | <2，函数f(z)解析，由于项1/(z. 1)不满足收敛条件，可以将其改写为(1/z)/(1. 1/z)，这样可以对1/z项作Taylor幂级数展开，而对第二项作Taylor幂级数展开，这样，Laurent级数可以写成
得出的Laurent级数为F2(z)=.1z k . ∞(2j)1 k 1z k ∑∑ k=.∞ k=0
（3）若z>2，两项的Taylor级数都不收敛，所以应该对1/z作Taylor级数展开这样得出的Laurent级数为F3(z)=∑(1 (2j)k 1) z k.11 k=.∞ 从这个例子可以看出，为使得原函数f(z)解析，应该将z平面分成三个不同的区域，分别得出Laurent级数展开，所以实际的Laurent级数应该是分段函数。仿照上述思路，可以对部分分式展开自动作区域分割，然后写出每个区域的Laurent级数展开，就此可以编写laurent_series()函数，以分段函数的形式得出各区域的Laurent级数。
· 176·高等应用数学问题的MATLAB求解（第四版）
例5-48再考虑例5-39中的有理函数，试关于z=0点求出其Laurent级数。解利用上面给出的函数laurent_series()可以直接得出关于z=0的Laurent级数
这样，得出的Laurent级数可以描述成如下的分段函数5 432 （1）若z<1，F1(z)=. 22818679z 221063z4981z 121z57z 1 || 6400000000 64000000 . 1600000 50000 . 40000 2000 （2）若1<| z | < 2 5 432216104333z1968701z34487z71z961zF2(z)=.172800000000 1728000000 . 43200000 675000 1080000 4911 1 1 1 . 27000 432z . 432z2 432z3 . 432z4 432z5 （3）若2<| z | < 4 3083895667z5 64031299z4 3265513z3 51527z2 330961z F3(z)= 172800000000 . 1728000000 43200000 . 337500 1080000 16549 529 1057 2113 4225 8449 . 27000 432z . 432z2 432z3 . 432z4 432z5 （4）若4<| z | < 5 342479989z5 62140157z4 56159897z3 252776089z2 226630597z F4(z)=. 56250000 2250000 . 450000 450000 . 90000 202363009 31883669 7198645 140679637 618454229 2709385813 ;18000144z2 .144z.144z3 144z4 .144z5 （5）若z>5，F5(z)=2 60 1028 13224 142048 1346208 11631876 || z . z2 z3 . z4 z5 . z6 z7 5.6.5 封闭曲线积分问题计算留数定理：考虑如下曲线积分∮ f(z)dz（5-6-12）其中，Γ为二维平面内的逆时针走行的任意形状的封闭曲线，且包围m个奇点pi(i=1,2,··· ,m)，则封闭曲线积分的值等于2.j乘以这些点的留数之和，即∮ m∑ f(z)dz=2.jRes[f(z),pi]（5-6-13）i=1如果积分线为顺时针走行的，则应该将被积函数乘以.1。例5-49试求函数f(z)在|z| =6逆时针曲线上的曲线积分。其中，函数f(z)为2z7 2z3 8 f(z)=z8 30z7 386z6 2772z5 12093z4 32598z3 52520z2 45600z 16000 解可以用例5-41中给出的方法求出其部分分式展开为13041 341863 7198933 16444 193046 1349779 11 1 (s 5)3 12(s 5)2 144(s 5) . 3(s 4)3 9(s 4)2 . 27(s 4) 9(s 2) 432(s 1) 可见，该函数的奇点为：p1=.1为单奇点，p2=.2亦为单奇点，p3=.4，p4=.5均为3重奇点。由上面的部分分式展开式可知，各个奇点的留数为部分分式展开的一次项的系数，这样可以得出所需的封闭曲线积分值为∮ []f(z)dz=2.j1798933 1349779 11 1 =4.j|z|=6 144.279432由第3章介绍的曲线积分方法，可以直接求出该积分。假设积分路径Γ由圆|z| =6 给出，可以表示为z=6cost j6sint，t[0,2.]，用下面语句直接积分也可以同样得出I=j4.如果要求解的问题是|z| =3的逆时针曲线积分，则在前面的求和式子中去除p3，p4两个奇点的留数（因为这两个奇点在封闭曲线外），可以终得出曲线积分的值为I=529j./216。用下面的语句直接求曲线积分也可以得出同样的结果
∮ 例5-50试求出下面的曲线积分[6] |z|=2 (z j)10(z 1 . 1)(z . 3) dz。解经过简单观察原函数，可以发现该函数在z=1，z=3处有单个奇点，在z=.j处有一个10重奇点，又因为给定的封闭曲线Γ为|z| =2正向圆周，所以z=1，z=.j在该圆周包围的范围内，z=3在该圆外，不必计算该留数，这样，原曲线积分的值可以用下面的语句直接求出，其值为(237/312500000 j779/78125000).。
根据文献[6]中给出的方法，手工求解时建议先计算z=3的留数，故得出整个环路积分值为..j/(3 j)10。二者之差为0，说明两个结果是完全一致的。由直接曲线积分方法也可以得出同样的结果
若曲线Γ的方程为|z| =4，则该曲线将z=1,3,.j三个奇点均包围在内，这时曲线积分应该和这三个留数的和有关，故可以用下面的语句求出曲线积分的值为0。
既然用符号运算的方法计算留数特别简单，所以没有必要再用间接的方法计算了，可以通过式（5-6-13）中给出的算法直接求出。该结果与直接曲线积分方法完全一致。
∫ 例5-51试计算封闭曲线积分I=| z |=1 z 2e1/zdz。解例5-46中曾经得出被积函数的Laurent级数，并指出z=0是原函数的本征奇点，这样该点处的· 178·高等应用数学问题的MATLAB求解（第四版）留数为c.1=1/6，所以根据留数定理，可以得出封闭曲线积分为I=2.jc.1=.j/3。还可以通过直接曲线积分验证得出的结果。
利用变量替换的方法可以将开区间积分变换成封闭曲线积分。例如，有的书上用变换的∫方法求解F(x)=∞ sin xdx这类无穷积分问题。而事实上，这样的积分直接用MATLAB中0 x 的int()函数更容易求解，所以本书不介绍这类方法。
5.7 差分方程的求解常系数线性差分方程的一般形式为y[(k n)T] a1y[(k n. 1)T] a2y[(k n. 2)T ] ··· any(kT)（5-7-1）=b1u[(k.d)T] b2u[(k.d.1)T] bmu[(k.d.m 1)T]··· 其中，T为采样周期。和微分方程描述的连续系统类似，这里的系数ai和bi也是常数，所以这类系统称为线性时不变离散系统。另外，对应系统的输入信号和输出信号也可以由u(kT)和y(kT)表示。u(kT)为第k个采样周期的输入信号，y(kT)为该时刻的输出信号。为方便起见，简记y(t)=y(kT)，且记y[(k i)T]为y(t i)，则前面的差分方程可以简记为y(t n) a1y(t n. 1) a2y(t n. 2) ··· any(t)（5-7-2）=b1u(t m. d) b2u(t m. d . 1) bm 1u(t. d)··· 本节将介绍线性差分方程的解析求解方法与数值求解方法，还将探讨基于仿真技术的非线性差分方程的数值求解方法。5.7.1 一般差分方程的解析求解方法前面给出了线性常系数差分方程的一般形式，若信号的初值y(0),y(1),··· ,y(n. 1)含有非零元素，则对式（5-7-2）两边进行z变换可以得出n∑.1n.2∑ z nY(z). z n.i y(i) a1zn.1Y(z). a1 z n.i y(i) anY(z)··· i=0i=0[]（5-7-3）m.1∑ =z.d b1zmU(z). b1 z n.i u(i) bm 1U(z)··· i=0由此得出Y(z)=(b1zm b2zm.1 ··· bm 1)z.dU(z) E(z)（5-7-4）zn a1zn.1 a2zn.2 an··· 其中，E(z)为输入、输出信号初值按式（5-4-2）中的变换计算出来的表达式n.1n.2n.3∑∑∑ E(z)=zn.i y(i). a1 z n.i y(i). a2 z n.i y(i).···. an.zy(0) .u(n)（5-7-5）i=0i=0i=0其中u.(n)=.b1 m.1z n.i u(i).···. bmzu(0)（5-7-6）∑ i=0对Y(z)进行z反变换则可以得出差分方程的解析解y(t)。根据前面的算法，可以编写出一般差分方程的通用求解函数
其调用语句为其中，A，B向量分别表示差分方程左侧和右侧的
系数向量，U为输入信号的z变换表达式，y0给出输出信号的初值向量，d为延迟步数，其默认值为0。调用该函数可以直接获得差分方程的解析解。该函数也可用于非首一化的差分方程。下面将给出具体实例来演示一般差分方程的求解方法。例5-52试求解差分方程48y(n 4). 76y(n 3) 44y(n 2). 11y(n 1) y(n)=2u(n 2) 3u(n 1) u(n)，其中，y(0)=1，y(1)=2，y(2)=0，y(3)=.1，且输入u(n)=(1/5)n ，试求出差分方程的解析解。解由给出的问题可以直接提取出A，B向量，将初始输出向量和输入信号送给计算机，再调用diff_eq()直接求解给出的差分方程
可以得出差分方程的解为()()()()()nn n nn4321 261 7521 1751 421 y(n)=5 3 . 52 . 54 35. 52 (n . 1) 其图形显示如图5-12所示，给出的几个初始点均在得出的解中。图5-12差分方程解的曲线表示5.7.2 线性时变差分方程的数值解法线性时变差分状态方程一般可以写成{ x(k 1)=F(k)x(k) G(k)u(k)x(0)=x0（5-7-7）y(k)=C(k)x(k) D(k)u(k),· 180·高等应用数学问题的MATLAB求解（第四版）可见，采用递推方法，则x(1)=F(0)x0 G(0)u(0)x(2)=F(1)x(1) G(1)u(1)=F(1)F(0)x0 F(1)G(0)u(0) G(1)u(1)终可以直接得出x(k)=F(k.1)F(k.2)F(0)x0 G(k.1)u(k.1)··· F(k. 1)G(k.2)u(k.2) ··· F(k.1) ··· F(0)G(0)u(0)（5-7-8）[]k.1k.1k.1∏ ∑∏ =F(j)x0 F(j)G(i)u(i)j=0i=0j=i 1若已知F(i),G(i)，则可以通过上面的递推算法直接求出离散状态方程的解。从数值求解的角度看，还可以用迭代方法求解该差分方程，即从已知的x(0)根据方程式（5-7-7）推出x(1)，再由x(1)计算x(2)，··· ，这样就可以递推地得出系统在各个时刻的状态。可见，迭代法更适合计算机实现。例5-53试求解离散线性时变差分方程[7] [ ][ ][ ][ ]x1(k 1)=01x1(k) sin(k./2)u(k)x2(k 1)1cos(k.)x2(k)1其中，x(0)=[1,1]T ，且u(k)=(.1)k ，k =0, 1, 2, 3, ··· 。解采用迭代方法，可以用下面的循环结构立即得出状态变量在各个时刻的值，如图5-13所示。
图5-13离散时变系统的响应曲线5.7.3 线性时不变差分方程的解法线性时不变差分方程有F(k)=··· =F(0)=F，G(k)=··· =G(0)=G，由式（5-7-8）可以立即得出k.1∑ x(k)=Fk x0 F k.i.1Gu(i)（5-7-9）i=0由于计算机数学语言并不能直接求出k是变量形式时Fk 的解析表达式，所以用上述表达式无法求出状态变量的解析解，必须考虑其他方法。再重新考虑式（5-7-7），其时不变形式可以写成{
x(k 1)=Fx(k) Gu(k)x(0)=x0（5-7-10）y(k)=Cx(k) Du(k),两端同时求z变换，由式（5-4-2）的性质可以得出X(z)=(zI. F).1这样可以推导出离散状态方程的解析解为[
zx0 GU(z). Gzu0]
（5-7-11）(zI. F).1 进一步观察上面的式子还可以发现，[F k = Z .1][]} 常数方阵F 的k 次方也可以通过z 反变换计算x0 Z .1{(zI. F).1（5-7-12）x(k)=Z.1GU(z). Gzu0z
][例5-54已知某离散系统的状态方程如下，试求出各个状态阶跃响应的解析解（5-7-13）z(zI. F).1.... 11/6 .5/4 3/4 .1/3 10 0 4 .. 0..x(k) .. 0 0 ..x(k 1)=u(k),x0=00 1/20 0 0 01/40 0 解直接套用下面的公式，则可以求解出状态方程的解析解为
从而得出各个状态的解析解为.. 48(1/3)k . 48(1/2)kk . 72(1/2)k . 24(1/2)kC2 k.1 48 144(1/3)k . 48(1/2)kk . 144(1/2)k . 48(1/2)kC2 k.1 48 kkx(k)=.... ....216(1/3). 192(1/2). 48(1/2)kCk2 .1 24 24(1/2)kk . 24(1/2)kCk2 .1 . 144(1/2)k 162(1/3)k 6 []。其实，事实上，得出结果中的nchoosek(n,k)是组合符号，其数学表示为Ckn =n!/(n . k)!k!C2 k.1 . 还可以进一步简化为(k. 1)(k . 2)/2，这样可以将得出的结果手工简化为. .12(8 k k2)(1/2)k 48(1/3)k 48 24(.8 k . k2)(1/2)k 144(1/3)k 48 24(.10 3k . k2)(1/2)k 216(1/3)kx(k)=... ... 12(.14 5k . k2)(1/2)k 162(1/3)k 6另外，因为原结果中只有Ck2 .1 项需要进一步简化，并对(1/2)k 合并同类项，还可以由下面的语句自动化简，这也将得出与手工化简一致的结果。
例5-55考虑例4-74中的Ak 计算问题，试用z反变换重新计算Ak。解由式（5-7-13）可见，矩阵乘方Ak 可以由下面的语句直接计算，其结果与例4-74中的也完全一致。· 182·高等应用数学问题的MATLAB求解（第四版）
5.7.4 一般非线性差分方程的数值求解方法假设已知差分方程的显式形式，即y(t)=f(t,y(t. 1),,y(t. n),u(t),,u(t. m))（5-7-14）······则可以通过递推的方法直接求解该方程，得出方程的数值解。例5-56考虑非线性差分方程y(t)=y(t. 1)2 1.1y(t. 2) 0.1u(t)，若输入信1 y(t. 1)2 0.2y(t. 2) 0.4y(t. 3)号为正弦函数u(t)=sint，采样周期为T=0.05s，试求解该方程的数值解。解引入一个存储向量y0，其三个分量y0(1)，y0(2)和y0(3)分别表示y(t. 3)，y(t. 2)和y(t. 1)，在每一步递推后更新一次y0向量。这样，用下面的循环结构就可以求解该方程，并绘制出输入信号和输出信号的曲线，如图5-14所示。可见，在正弦信号激励下，非线性系统的输出会产生畸变，这与线性系统响应是不同的。
输出信号y(t)输入信号u(t)图5-14非线性离散差分方程数值解曲线5.8 习题（1）对下列的函数f(t)进行Laplace变换：.fa(t)=sinαt/t，.fb(t)=t5 sinαt，.fc(t)=t8 cosαt，.fd(t)=t6e.t ，.fe(t)=5e.at t4e.at 8e.2t ，.ff(t)=e.t sin(αt θ)，.fg(t)=e.12t 6e9t。
（2）对上面的结果作Laplace反变换，看看能不能还原给定的函数。
（3）下面两个公式也是Laplace变换的性质，试选择不同的常数n，验证这两条性质。

. L [tnf(t)] =(.1)n dnL ds[fn (t)]， . L [tn.1/2] = √.(22n n . 1)! s.n.1/2。（4）对下面的F(s)式进行Laplace反变换：.Fa(s)=1 ，.Fb(s)=√s. a .√s . b，.Fc(s)=lns . a ，√s(s2 . a2)(s b) s . b .Fd(s)=√s(s1 a) ，.Fe(s)=s3 3 a2 a3 ，.Ff(s)=(s . s7 1)8 ， 22 .Fg(s)=lns a，.Fh(s)=s2 ∏ 3s 8 ，.Fi(s)=1 s α 。s2 b28(s i)2 s . α i=1（5）证明，对s的非整数次方，下面的Laplace变换的公式成立。.对不同的γ取值，试证明L[tΓ]= Γ(γ 1) ， sΓ 1[ ] (√).对任意a>0，试证明L1= . es/aerfc s/a。√t (1 at) a （6）Laplace变换的一个重要应用是求解常系数线性微分方程，可以利用当函数和各阶导数的零初始值下L[dnf(t)/dtn]= snL[f(t)]这一性质对微分方程进行Laplace变换的方法去求解微分方程，非零初值问题也可以利用相应方法求解。试使用这样的方法求解下面的微分方程。.y′′(t) 3y′(t) 2y(t)=e.t ,y(0)=y′(0)=0，
.y′′. y=4sint 5cos2t,y(0)=.1,y′(0)=.2，

. .x ′′ . x y z =0 . . .x y′′. y z=0x(0)=1,x′(0)=y(0)=y′(0)=z(0)=z′(0)=0。. x y z ′′ . z =0, （7）假设某分数阶系统是由两个子模型G1(s)和G2(s)并联而成，则系统的总模型可以由G(s)=G1(s) G2(s)计算出来。试对下面的两个子模型并联的总系统绘制出阶跃响应曲线。(s0.4 2)0.8s0.4 0.6s 3G1(s)=√s(s2 3s0.9 4)0.3,G2(s)=(s0.5 3s0.4 5)0.7（8）系统模型G1(s)，G2(s)串联连接构造的总系统可以由G(s)=G2(s)G1(s)表示，试求出上例两个子传递函数串联后的阶跃响应曲线。
（9）试求出下面函数的Fourier变换，对得出的结果再进行Fourier反变换，观察能否得出原函数。

.f(x)=x2(3. . 2|x|),0.x.2.，.f(t)=t2(t . 2.)2 ,0.t.2.，2 .f(t)=e.t, .l.t.l，.f(t)=te.|t|, ...t..。（10）试求出下面函数的Fourier正弦和余弦变换，并用Fourier正弦、余弦反变换对得出的结果进行处理，观察是否能还原出原函数。.f(t)=e.t lnt，.f(x)=cos x2 ，.f(x)=ln√1 ，x 1 x2
.对任意的a>0，f(x)=x(a2 . x2)，.f(x)=coskx。（11）试求函数.f(x)=ekx ，.f(x)=x3 的离散Fourier正弦、余弦变换。· 184·高等应用数学问题的MATLAB求解（第四版）sin(alnx),x.1（12）试对分段函数f(x)={ 0, 其他进行Mellin变换。
（13）请将下述时域序列函数f(kT)进行z变换，并对结果进行反变换检验。

.fa(kT)=cos(kaT)，.fb(kT)=(kT)2e.akT ，.fc(kT)=1(akT. 1 e.akT )， a.fd(kT)=e.akT . e.bkT ，.fe(kT)=sin(αkT)，.ff(kT)=1. e.akT (1 akT)。（14）已知下述各个z变换表达式F(z)，试对它们分别进行z反变换。10zz2 z .Fa(z)=，.Fb(z)=，.Fc(z)=，(z . 1)(z . 2) (z . 0.8)(z. 0.1)(z. a)(z. 1)2 z.1(1 . e.aT ) Az[zcosβ. cos(αT. β)].Fd(z)=(1 . z.1)(1 . z.1e.aT ) ，.Fe(z)=z2 . 2zcos(αT) 1。（15）对下面的Laplace变换式求出相应的z变换，并对结果进行检验。.G(s)=s2(sb a) ，.G(s)=s2(sb a)2 1 . s e.2s 。（16）已知函数G(s)=1 ，若用s = 2(z . 1) 替换G(s)，则可以得出函数H(z)。这样的变换(s 1)3 T (z 1) 1 Ts/2又称为双线性变换。若T=1/2，试求出H(z)。若对结果进行新变换z=1 . Ts/2，则得出双线性反变换的结果。试观察这样的反变换是否能恢复原函数。（17）试用计算机证明Z{1 . e.akT [cos(bkT) a sin(bkT)]} = z(Az B)b (z . 1)(z2 . 2e.aT cos(bT)z e.2aT) 式中，A =1 . e.aT cos(bT). ab e.aT sin(bT)，B=e.2aT ab e.aT sin(bT). e.aT cos(bT)。（18）试判定下面的多项式组是否互质。如果非互质，试化简B(s)/A(s)。
.B(x)=.3x4 x5 . 11x3 51x2 . 62x 24，A(x)=x7 . 12x6 26x5 140x4 . 471x3 . 248x2 1284x . 720，

.B(x)=3x6 . 36x5 120x4 90x3 . 1203x2 2106x . 1080，A(x)=x9 15x8 79x7 127x6 . 359x5 . 1955x4 . 3699x3 . 3587x2 . 1782x . 360。

（19）试绘制如下复变函数的Riemann曲面。.f(z)=zcosz2 ，.f(z)=ze.z 2 (cosz. sinz)。
（20）记z=x jy，其中，x和y满足x2 (y . 1)2 =1。很显然，z是一个圆。试利用w=1/z替换将该圆映射到w平面上，这个映射出来的图形是什么？
（21）考虑函数f(z)=z2 4z 3 e.5z ，试找出函数全部的极点，找到各个极点

z5 4z4 3z3 2z2 5z 2 的重数，并计算它们的留数。（22）试求出下面有理函数的部分分式展开。3x4 . 21x3 45x2 . 39x 12 .f(x)=，
.f(s)=，
s8 21s7 181s6 839s5 2330s4 4108s3 4620s2 3100s 1000 3x6 . 36x5 120x4 90x3 . 1203x2 2106x . 1080
.f(x)=，

x7 15x6 96x5 340x4 720x3 912x2 640x 192s 5 x7 13x6 52x5 10x4 . 431x3 . 1103x2 . 1062x . 360(x2 4x 3)e.5x .f(x)=x5 7x4 . 2x3 . 100x2 . 232x . 160 。（23）试求出下列函数奇点处的留数。.f(z)=1 . sinze.2z (z 4 10z 3 35z 2 50z 24)， z7 sin(z. ./3)2 .f(z)=(z . 3)4 (sinz. e.3z)，.f(z)=(1 . cos2z)(1. e.z ) 。z4 5z3 9z2 7z 2 z3 sin z （24）试写出下面函数的Laurent级数并找出其留数。2 .f(z)=ze.1/z [sin(1/z). cos(1/z)]，.f(z)=z5 cos(1/z2)。（25）试写出下面函数的Laurent级数。3 1 1 155 f(z)= （26）求下面的封闭环路积分。z . 1(z . 1)2 z . 2(z . 2)2 z i z . i ∮ 15 . (z2 1)z2(z4 2)3 dz，其中，Γ为|z| =3 正向圆周，Γ ∮ 3z. 1 z e1/zdz，其中，Γ为|z| =2 正向圆周，∮Γ 2 . ∫Γ cosz(1. e.z )sin(3z 2)dz，其中，Γ为z=1的正向圆周，z sin z || . z . 2 dz。| z |=2 z3(z . 1)(z . 3) （27）复平面映射中有一个很有趣的数学分支——分形，下面三个习题均是这方面的内容。任意选定一个二维平面上的初始点坐标(x0,y0)，假设可以生成一个在[0,1]区间上均匀分布的随机数γi，那么根据其取值的大小，可以按下面的公式生成一个新的坐标点(x1,y1)。(x1,y1). ... .. x1=0,y1=y0/2,如果γi<0.05x1=0.42(x0. y0),y1=0.2 0.42(x0 y0),如果0.05.γi<0.45x1=0.42(x0 y0),y1=0.2. 0.42(x0. y0),如果0.45.γi<0.85x1=0.1x0,y1=0.2 0.1y0,其他情况试用该递推方法计算10000个点，并将这些计算点用圆点表示，观察得出的结果。（28）选定一个复数c，在对(xm,ym)到(xM,yM)平面区域内每个点z0=x0 jy0做如下映射zn 1=zn 2 c之后，如果变换的点仍为有界的，则再继续上述的映射。进行若干次映射后，可以得出一个新的复数，把它进行变换后赋给z(x0,y0)。试由该映射生成数据并绘制Julia图。（提示：可以参考本书提供的julia.m文件，区域.1.3.x,y.1.3，网格点300，迭代次数30，且.c0=0.27334 j0.00742，.c0=.0.7）
（29）在zn 1=zn 2 c映射中，如果z0遍取区域内(xm,ym)到(xM,yM)所有的点，则可以得出各个点处的测度，从而绘制出Mandelbrot图。（提示：可以参考本书提供的mandelbrot.m函数，选择.2.x.0.5，.1.25.y.1.35，网格点300，迭代次数200，c=0，收敛阈值a=2；若选择区域.0.74547.x..0.74538，0.11298.y.0.11304，重新绘制）
· 186·高等应用数学问题的MATLAB求解（第四版）
（30）试求解下面给出的差分方程模型。

.72y(t) 102y(t. 1) 53y(t. 2) 12y(t. 3) y(t. 4)=12u(t) 7u(t. 1)，u(t)为阶跃信号，且y(.3)=1,y(.2) = .1，y(.1)=y(0)=0，.y(t). 0.6y(t. 1) 0.12y(t. 2) 0.008y(t. 3)=u(t)，u(t)=e.0.1t，且y(t)初值为0。y(t. 2) 4y(t. 1) 2u(t)（31）试求解非线性差分方程y(t)=u(t) y(t. 2) 3y 2(t . 1) 1 y2(t . 2) y2(t . 1) ，且t.0时y(t)=0，u(t)=e.0.2t。（32）Fibonacci序列a(1)=a(2)=1，a(t 2)=a(t) a(t 1)，t=1,2,事实上是一个线性差分方程，试求出通项a(t)的解析解。··· ，（33）已知某离散系统的状态方程模型如下，且xT(0) = [1, .1]，试求该系统阶跃响应的解析解，并比较数值解。[. 01 .0.16.1]x(t) [1 1].x(t 1)=u(t)，... 11/6 .1/4 25/24 .2 11 .1 2 ...1..x(t) ..1/2 .3/8 1/4 ..u(t)。.x(t 1)=01 .1 01 .3/4 0 0

参考文献[1]ValsaJ，Bran.ikL.ApproximateformulaefornumericalinversionofLaplacetransforms.InternationalJournalofNumericalModelling:ElectronicNetworks,DevicesandFields，1998，11（3）：153–166.[2]ValsaJ.NumericalinversionofLaplacetransformsinMATLAB,MATLABCentralFileID:#32824，2011.[3]CallierFM，WinkinJ.In.nitedimensionalsystemtransferfunctions.CurtainRF，BensoussanA，LionsJL，eds.，AnalysisandOptimizationofSystems:StateandFrequencyDomainApproachesforIn.nite-DimensionalSystems.Berlin：Springer-Verlag，1993.[4]《数学手册》编写组.数学手册.北京：人民教育出版社，1979.
[5]宋叔尼，孙涛，张国伟.复变函数与积分变换.北京：科学出版社，2006.
[6]西安交通大学高等数学教研室编.复变函数（第二版）.北京：人民教育出版社，1981.
[7]郑大钟.线性系统理论（第2版）.北京：清华大学出版社，2002.