现代科学运算-matlab语言与应用

$\qquad \qquad$ 东北大学 http://www.icourses.cn/home/ (薛定宇)
《高等应用数学问题的MATLAB求解（第四版）》
代码可在matlab r2016b 运行。

05 积分变换与复变函数问题的计算机求解

05.01 Laplace变换

Laplace变换及其反变换
Laplace变换是科学与工程运算中的重要数学变换
很多学科都源于此变换控制系统的传递函数模型直接由Laplace变换定义

Lplace变换及反变换定义与性质
Laplace变换(拉氏变换)的数学性质
$\mathscr{L}[f(t)] = \int_0^{\infty} f(t) e^{-st} dt = F(s)$
Laplace变换性质
线性性质：$\mathscr{L}[af(t) \pm bg(t)] = a \mathscr{L}[f(t)] \pm b \mathscr{L}[g(t)]$；其中，a与b均为标量
时域平移性质：$\mathscr{L}[f(t-a)] = e^{-as}F(s)$
s-域平移性质：$\mathscr{L}[e^{-at}f(t)] = F(s+a)$
微分性质：$\mathscr{L}[df(t)/dt] = sF(s) - f(0^+)$
n阶微分：$\mathscr{L}[\dfrac{d^n}{d t^n} f(t)] = s^nF(s) - s^{n-1} f(0^+) - s^{n-2} \dfrac{d f(0^+)}{dt} - \cdots - \dfrac{d^{n-1} f(0^+)}{d t^{n-1}}$
当初值为0时，则 $\mathscr{L}[\dfrac{d^n f(t)}{d t^n}] = s^n F(s)$

积分性质：
零初始条件: $\mathscr{L}[\int_0^t f(\tau) d \tau] = \dfrac{F(s)}{s}$
多重积分: $\mathscr{L}[\int_0^t \cdots \int_0^t f(\tau) d \tau^n] = \dfrac{F(s)}{s^n}$
初值性质： $\lim \limits_{t \to 0} f(t) = \lim \limits_{s \to \infty} s F(s)$
终值性质：$如果F(s)没有s \geq 0 的几点，\lim \limits_{t \to \infty} f(t) = \lim \limits_{s \to 0} sF(s)$
卷积性质: $\mathscr{L}[f(t)*g(t)] = \mathscr{L}[f(t)] \mathscr{L}[g(t)]$
其中卷积算子 * 的定义：
$f(t)*g(t) = \int_0^t f(\tau) g(t-\tau) d\tau = \int_0^t f(t - \tau) g(\tau) d\tau$
其它性质：
$\mathscr{L}[t^n f(t)] = (-1)^n \dfrac{d ^n F(s)}{d s^n}$
$\mathscr{L}[\dfrac{f(t)}{t^n}] = \int_s^{\infty} \cdots \int_s^{\infty} F(S) d s^n$
Laplace反变换：
$f(t) = \mathscr{L}^{-1}[F(s)] = \dfrac{1}{j2\pi} \int_{\sigma-j\infty}^{\sigma+j\infty} F(s) e^{st} ds$
其中，$\sigma$ 大于函数F(s)奇点的实部

Laplace变换的计算机求解
Laplace变换问题的求解步骤:
$\mathscr{L}[f(t)] = \int_0^{\infty} f(t) e^{-st} dt = F(s)$
符号定义变量t, 再定义时域函数
直接调用laplace()函数
采用默认的t为时域变量
F = laplace(fun)
用户指定时域变量v和复域变量名u
F = laplace(fun, v, u)
调用pretty()或latex()进一步处理结果

Laplace反变换
Laplace反变换计算
函数调用格式：
采用默认的t为时域变量
f = ilaplace(fun)
用户指定时域变量u和复域变量名v
f = ilaplace(fun, u, v)

例5-1 Laplace变换
给定函数 $f(t) = t^2 e^{-2t} \sin (t + \pi)$

% Laplace变换
syms t;
f = t^2* exp(-2*t)*sin(t + pi);
F = laplace(f)
% 化简得出的结果
simplify(F),
pretty(ans)

例5-4 函数及导数的变换
给定 $f(t) = t^2 e^{-2t} \sin(t+\pi)$
推导 $s^5\mathscr{L}[f(t)]$ 和 $\mathscr{L}[d^5 f(t) / d t^5]$ 之间的关系

syms t s;
f = t^2*exp(-2*t)*sin(t+pi);
F = simplify(laplace(diff(f, t, 5)))
F0 = laplace(f);
simplify(F - s^5 * F0)
% 误差是什么？

初值的影响
考虑到初值条件：

ss = 0;
f1 = f;
for i = 4: -1 : 0
  ss = ss - s^i * subs(f1, t, 0);
  f1 = diff(f1, t);
end;
ss

回忆公式$\mathscr{L}[\dfrac{d^n}{d t^n} f(t)] = s^nF(s) - s^{n-1} f(0^+) - s^{n-2} \dfrac{d f(0^+)}{dt} - \cdots - \dfrac{d^{n-1} f(0^+)}{d t^{n-1}}$

例5-5 Laplace变换导数公式推导
前面给出的微分公式不易记忆
试推导出 $\mathscr{L}[d^2 f(t) /d t^2]$ 的微分公式

syms t y(t);
laplace(diff(y, t, 2))
% 函数8阶导数的laplace变换
laplace(diff(y, t, 8))

例5-6 导数的Laplace变换
已知 $f(t) = e^{-5t} \cos(2t + 1) + t$
求 $\mathscr{L}\Big[\dfrac{d^5 f(t)}{d t^5} \Big]$

syms t;
f = exp(-5*t)*cos(2*t+1) + t;
F = laplace(diff(f, t, 5));
F = simplify(F);
pretty(F)

05.02 数值Laplace变换

Laplace变换问题的数值求解
为什么要研究数值求解？
原函数过于复杂，没有解析解
无需解析表达式，得到图形即可
Laplace反变换数值解析工具
Juraj Valsa 函数: Math Works File Exchange 下载
调用方法: [t, y] = INVLAP(f, $t_0$, $t_f$, N, other pars)
一般不建议给出other pars，采用默认
注意：函数本身bug，初始时刻不能为0

% INVLAP ?Numerical Inversion of Laplace Transforms 
function [radt,ft]=INVLAP(Fs,tini,tend,nnt,a,ns,nd);
% Fs is formula for F(s) as a string
% tini, tend are limits of the solution interval
% nnt is total number of time instants
% a, ns, nd are parameters of the method 
% if not given, the method uses implicit values a=6, ns=20, nd=19
% it is recommended to preserve a=6
% increasing ns and nd leads to lower error
% an example of function calling  
% [t,ft]=INVLAP('s/(s^2+4*pi^2)',0,10,1001);
% to plot the graph of results write plot(t,ft), grid on, zoom on
FF=strrep(strrep(strrep(Fs,'*','.*'),'/','./'),'^','.^');
if nargin==4
  a=6; ns=20; nd=19;  end;    % implicit parameters
radt=linspace(tini,tend,nnt); % time vector
if tini==0  radt=radt(2:1:nnt);  end;  % t=0 is not allowed
tic                 % measure the CPU time
for n=1:ns+1+nd               % prepare necessary coefficients
   alfa(n)=a+(n-1)*pi*j;
   beta(n)=-exp(a)*(-1)^n;
end;
n=1:nd;
bdif=fliplr(cumsum(gamma(nd+1)./gamma(nd+2-n)./gamma(n)))./2^nd;
beta(ns+2:ns+1+nd)=beta(ns+2:ns+1+nd).*bdif;
beta(1)=beta(1)/2;
for kt=1:nnt                  % cycle for time t
   tt=radt(kt);
   s=alfa/tt;                 % complex frequency s
   bt=beta/tt;
   btF=bt.*eval(FF);          % functional value F(s)
   ft(kt)=sum(real(btF));     % original f(tt)
end;
toc

复杂系统的输出计算
系统框图
$U(s) \to G(s) \to Y(s)$
用数值方法计算输出
输出信号 Y(s) = G(s)U(s)
系统复杂，如分数阶系统，解析解不能求出
输入为复杂信号，不能求其Laplace变换
输入信号只给出数据点，没有函数

Laplace变换与反变换的数值解

% 扩展代码
function [t,y]=INVLAP_new(G,t0,tn,N,H,tx,ux)
% INVLAP_new - updated version of INVLAP for closed-loop system with any inputs
%
%   [t,y]=INVLAP_new(G,t0,tn,N)
%   [t,y]=INVLAP_new(G,t0,tn,N,H)
%   [t,y]=INVLAP_new(G,t0,tn,N,H,u)
%   [t,y]=INVLAP_new(G,t0,tn,N,H,tx,ux)
%
%   G - the string representation of the open-loop model
%   t0, tn, N - the time interval and number of points in the interval
%   H - the string representation of the feedback model
%   u - the function handle of the input signal
%   tx, ux - the samples of time and input signal
%   y, t - the output and time vectors

% Copyright (c) Dingyu Xue, Northeastern University, China
% Last modified 13 June, 2017 
G=add_dots(G); if nargin<=5, tx='1'; end, if nargin<=4, H=0; end
if ischar(H), H=add_dots(H); end
if ischar(tx), tx=add_dots(tx); end
a=6; ns=20; nd=19; t=linspace(t0,tn,N);
if t0==0, t=t(2:N); N=N-1; end, 
n=1:ns+1+nd; alfa=a+(n-1)*pi*j; bet=-exp(a)*(-1).^n;
n=1:nd; bet(1)=bet(1)/2;
bdif=fliplr(cumsum(gamma(nd+1)./gamma(nd+2-n)./gamma(n)))./2^nd;
bet(ns+2:ns+1+nd)=bet(ns+2:ns+1+nd).*bdif;
if isnumeric(H), H=num2str(H); end
for i=1:N
   tt=t(i); s=alfa/tt; bt=bet/tt; sG=eval(G); sH=eval(H);
   if ischar(tx), sU=eval(tx);
   else
      if isnumeric(tx),
         f=@(x)interp1(tx,ux,x,'spline').*exp(-s.*x);
      else, f=@(x)tx(x).*exp(-s.*x); end
      sU=integral(f,t0,tn,'ArrayValued',true);
   end
   btF=bt.*sG./(1+sG.*sH).*sU; y(i)=sum(real(btF));
end
% remove and add back dots uniformly
function F=add_dots(F)
F=strrep(strrep(strrep(F,'.*','*'),'./','/'),'.^','^');
F=strrep(strrep(strrep(F,'*','.*'),'/','./'),'^','.^');

新函数的调用方法
调用格式
[t, y] = INVLAP_new(G, $t_0, t_n,$ N)
[t, y] = INVLAP_new(G, $t_0, t_n,$ N, H)
[t, y] = INVLAP_new(G, $t_0, t_n,$ N, H, u)
[t, y] = INVLAP_new(G, $t_0, t_n, N, H, t_x, u_x$)
负反馈系统

例5-8 分数阶模型的数值Laplace反变换
分数阶模型
$G(s) = \dfrac{s^{0.4} + 0.4s^{0.2} + 0.5}{\sqrt{s}(s^{0.2} + 0.02s^{0.1} + 0.6)^{0.4}(s^{0.3} + 0.5)^{0.6}}$
Laplace反变换解析解不存在，只能求数值解

G = ['(s^0.4 + 0.4*s^0.2 + 0.5)/' ...
  'sqrt(s)/(s^0.2+0.02*s^0.1+0.6)^0.4/(s^0.3+0.5)^0.6'];
[t, y] = INVLAP_new(G, 0, 1, 10000);
plot(t, y)

例5-9 分数阶模型时域响应
分数阶模型
$G(s) = \dfrac{s^{0.4} + 0.4s^{0.2} + 0.5}{\sqrt{s}(s^{0.2} + 0.02s^{0.1} + 0.6)^{0.4}(s^{0.3} + 0.5)^{0.6}}$
输入信号: $u(t) = e^{-0.3t} \sin t^2$

f = @(t)exp(-0.3*t).* sin(t.^2);
G = ['(s^0.4 + 0.4*s^0.2 + 0.5)/' ...
  'sqrt(s)/(s^0.2+0.02*s^0.1+0.6)^0.4/(s^0.3+0.5)^0.6'];
tic,
[t, y] = INVLAP_new(G, 0, 15, 400, 0, f);
toc
plot(t, y)
figure();
x0 = 0: 0.2: 15;
u0 = exp(-0.3*x0).* sin(x0.^2);
tic,
[t, y] = INVLAP_new(G, 0, 15, 200, 0, x0, u0);
toc
plot(t, y)

例5-11 闭环系统阶跃响应
开环无理传递函数模型
$G(s) = \Big[ \dfrac{\sinh(0.1\sqrt{s})}{0.1 \sqrt{s}} \Big]^2 \dfrac{1}{\sqrt{s} \sinh(\sqrt{s})}$
单位负反馈闭环响应计算

G = ['(sinh(0.1*sqrt(s))/0.1/sqrt(s))^2' ...
  '/sqrt(s)/sinh(sqrt(s))'];
[t, y] = INVLAP_new(G, 0, 10, 1000, 1, '1/s');
plot(t, y)

05.03 Fourier变换

Fourier变换及反变换定义与性质
Fourier变换定义:$\mathscr{F}[f(t)] = \int_{-\infty}^{\infty} f(t) e^{-j\omega t} dt$
Fourier反变换定义:$f(t) = \mathscr{F}^{-1}[F(\omega)] = \dfrac{1}{2 \pi} \int_{-\infty}^{\infty} F(\omega) e^{j\omega t} d \omega$

Fourier变换性质
线性性质:其中a与b均为标量
$\mathscr{F}[af(t) \pm bg(t)] = a\mathscr{F}[f(t)] \pm b\mathscr{F}[g(t)]$
平移性质:$\mathscr{F}[f(t \pm a)] = e^{\pm ja \omega} F(\omega)$
复域平移性质：$\mathscr{F}[e^{\pm j a t} f(t)] = F(\omega \mp a)$
微分性质：$\mathscr{F}[\dfrac{df(t)}{dt}] = j \omega F(\omega)$
n阶微分的Fourier变换:
$\mathscr{F}\Big[\dfrac{d^n}{d t^n} f(t) \Big] = (j\omega)^n \mathscr{F}[f(t)]$
积分性质：$\mathscr{F}\Big[\int_{-\infty}^{\infty} f(\tau) d\tau \Big] = \dfrac{F(\omega)}{(j \omega)}$
n重积分的Fourier变换:$\mathscr{F}\Big[\int_{-\infty}^{\infty} \cdots \int_{-\infty}^{\infty} f(\tau) d \tau^n \Big] = \dfrac{\mathscr{F}[f(t)]}{(j \omega)^n}$

Fourier变换的计算机求解
Fourier变换的函数调用格式
默认变量：F = fourier(fun)
$v \to u$: F = fourier(fun, v, u)
Fourier反变换的函数调用格式
默认变量: f = ifourier(fun)
$u \to v$: f = ifourier(fun, u, v)

例5-11 Fourier变换
给定 $f(t) = 1/(t^2 + a^2), a > 0$

% Fourier变换
syms t w;
syms a positive,
f = 1/(t^2+a^2);
F = fourier(f, t, w)
% Fourier反变换
f1 = ifourier(F)

例5-12 Fourier变换
给定 $f(t) = \sin^2(at) / t, a > 0$

% Fourier变换
syms t w; 
syms a positive;
f = sin(a*t)^2/t;
F = fourier(f, t, w)

得出结果：
$F = -\dfrac{\pi heaviside(-2a - w) li}{2} - \dfrac{\pi heaviside(2a - w) li}{2} + \pi heaviside(-w) li$
手工化简结果：
$\mathscr{F}[f(t)] = \left \{ \begin{array}{l} 0, \qquad |\omega| > 2a \\ -j \pi \mathbb{sign}(\omega) /2, |\omega| < 2a \end{array} \right.$

例5-13 Fourier变换
给定 $f(t) = e^{-a|t|}/\sqrt{|t|}$
Fourier变换的数学手册结果
$F = \dfrac{\sqrt{\sqrt{\omega^2 + a^2} + a}}{\sqrt{\omega^2 + a^2}}$
使用fourier()命令(早期版本不可用)

syms w t;
syms a positive;
f = exp(-a*abs(t))/sqrt(abs(t));
F = fourier(f, t, w),
pretty(F)

Fourier 正弦和余弦变换
Fourier正弦正反变换的一般定义
$\mathscr{F}[f(t)] = \int_0^{\infty} f(t) \sin \omega t dt = F_s(\omega)$
$\mathscr{F}[f(t)] = \int_0^{\infty} f(t) \cos \omega t dt = F_c(\omega)$
$\mathscr{F}_s^{-1}[F_s(t)] = \dfrac{2}{\pi} \int_{-\infty}^{\infty} F_s(\omega) \sin(\omega t) d \omega$
$\mathscr{F}_c^{-1}[F_c(t)] = \dfrac{2}{\pi} \int_{-\infty}^{\infty} F_c(\omega) \cos(\omega t) d \omega$

例5-13 余弦Fourier变换
给定 $f(t) = t^n e^{-at}, a > 0, n = 1, 2, \cdots 8$
试求出其余弦Fourier变换

syms t; syms w real;
syms a positive;
for i = 1:8
  f = t^i*exp(-a*t);
  F = simplify(int(f*cos(w*t), t, 0, inf)),
end

例5.14 分段函数的变换
分段函数 $f(t) = \left \{ \begin{array}{l}\cos(t), 0 < t < a \\ 0, others \end{array} \right.$
试求其Fourier余弦变换

syms t w;
syms a positive;
f = cos(t);
F = simplify(int(f*cos(w*t), t, 0, a))

离散Fourier正弦、余弦变换
离散Fourier正、余弦变换
$F_s(k) = \int_0^a f(t) \sin \dfrac{k \pi t}{a} dt$
$F_c(k) = \int_0^a f(t) \cos \dfrac{k \pi t}{a} dt$
离散Fouriere正、余弦反变换
$f(t) = \dfrac{2}{a} \sum \limits_{k=1}^{\infty} F_s(k) \sin \dfrac{k \pi t}{a}$
$f(t) = \dfrac{1}{a}F_c(0) + \dfrac{2}{a} \sum \limits_{k=1}^{\infty} F_c(k) \cos \dfrac{k \pi t}{a}$

例5-15 分段函数
给定 $f(t) = \left \{ \begin{array}{l}t, t \leq a/2, \\ a - t, t > a/2 \end{array} \right.$
其中a > 0,计算其离散Fourier正弦变换

syms t k;
assumeAlso(k, 'integer');
syms a positive;
f1 = t;
f2 = a-t;
Fs = int(f1*sin(k*pi*t/a), t, 0, a/2) + ...
  int(f2*sin(k*pi*t/a), t, a/2, a);
simplify(Fs)
% 由分段函数直接求
f = piecewise('t <= a/2', 't', 't > a/2', 'a-t');
Fs = simplify(int(f*sin(k*pi*t/a), t, 0, a))

快速Fourier变换
离散Fourier变换
$X(k) = \sum \limits_{i=1}^{N} x_i e^{-2 \pi j (k-1)(i-1)/N},$ where $1 \leq k \leq N$
反变换
$x(k) = \dfrac{1}{N} \sum \limits_{i=1}^{N} X(i) e^{2 \pi j (k-1)(i-1)/N},$ where $1 \leq k \leq N$

快速Fourier变换(Fast Fourier transform,FFT)技术是求离散Fourier变换最实用、也最通用的方法
matlab直接求解: f = fft(x) $\qquad \hat x$ = ifft(f)
特点：高效、快速；任意序列长度，长度不要求为$2^n$

例5-16 给定信号的FFT
原型函数：$x(t) = 12 \sin (2 \pi \times t + \pi /4) + 5 \cos (2 \pi \times 4 t)$
先FFT，再FFT反变换，看看能否还原

h = 0.01;
t = 0: h:10;
x = 12*sin(2*pi*t + pi/4) + 5*cos(2*pi*4*t);
X = fft(x);
f = t/h/10;
L = floor(length(f)/2);
stem(f(1:L), abs(X(1:L))), xlim([0, 10])
figure();
ix = real(ifft(X));
plot(t, x, t, ix, ':');
xlim([0, 10])

多为FFT
二维FFT
正变换: fft2()
反变换: ifft2()
多维FFT
正变换: fftn()
反变换: ifftn()

05.04 Mellin变换与Hankel变换

Mellin变换
Mellin变换定义：$\mathscr{M}[f(x)] = \int_0^{\infty} f(x) x^{z-1} dx = M(z)$
Mellin反变换定义：$f(x) = \mathscr{M}^{-1}[M(z)] = \dfrac{1}{j2\pi} \int_{c-j\infty}^{c+j\infty} M(z) x^{-z} dz$

例5-19 Mellin变换
给定 $f(t) = \ln t / (t + a) , a > 0$
试求其Mellin变换

syms t z;
syms a positive;
f = log(t)/(t+a);
M = int(f*t^{z-1}, t, 0, inf)

结果：$\left \{ \begin{array}{l}\dfrac{x^2\cos(\pi z)}{cos(\pi z)^2 - 1}, \quad a = 1 \land z \neq \dfrac{1}{2} \land z \in (0, 1) \\ \int_0^{\infty} \dfrac{t^{z-1} \ln (t)}{a + t} dt, \quad a \neq 1 \lor z = \dfrac{1}{2} \lor z \notin (0, 1) \end{array} \right.$

大部分Mellin变换的解析解不存在

数值Mellin变换
Mellin变换：$\mathscr{M}[f(x)] = \int_0^{\infty} f(x) x^{z-1} dx = M(z)$
数值Mellin函数: F = mellin_trans(f, z, pars)

function F = mellin_trans(f, z, varargin)
  f1 = @(x)f(x).*x.^(z-1);
  F = integral(f1, 0, Inf, 'ArrayValued', true, varargin{:});

例5-20 给定函数的Mellin变换
原函数:$f(x) = \sin(3x^{0.8})/(x+2)^{1.5}$
Mellin变换的解析解不存在，只能求数值解

f = @(x)sin(3*x.^0.8)./(x+2).^1.5;
z = 0:0.05:1;
F = mellin_trans(f, z);
plot(z, F)

Hankel变换及求解
v阶Hankel变换的数学表达式为
$\mathscr{H}[f(t)] = \int_0^{\infty} tf(t) J_v(\omega t) dt = H_v(\omega)$
其中，$J_v(\cdot)$为Bessel函数，J = besselj(v, z)
Hankel反变换公式：
$\mathscr{H}^{-1}[H(\omega)] = \int_0^{\infty} \omega H_v(\omega) J_v(\omega t) d \omega$

例5-21 直接积分求Hankel变换
已知函数 $f(t) = e^{-a^2t^2/2}$
求 0 阶 Hankel 变换

syms t w a positive;
f = exp(-a^2*t^2/2);
F = int(f*t*besselj(0, w*t), t, 0, inf);
F = simplify(F)
f1 = int(w*F*besselj(0, w*t), w, 0, inf)

数值Hankel变换
Hankel变换公式
$\mathscr{H}[f(t)] = \int_0^{\infty} t f(t) J_v(\omega t) d t = H_v(\omega)$
matlab求解：H = hankel_trans(f, w, mu, pars)

function H = hankel_trans(f, w, mu, varargin)
  F = @(t)t.*f(t).*besselj(mu, w*t);
  H = integral(F, 0, Inf, 'ArrayValued', true, varargin{:});

例5-22 Hankel数值变换
已知函数: $f(t) = e^{-a^2t^2/2}, a = 2$
选择不同阶次v，用数值方法求出Hankel变换

syms t w a positive;
f = exp(-a^2*t^2/2);
F = int(f*t*besselj(0, w*t), t, 0, inf);
F = simplify(F)
F1 = subs(F, a, 2);
ezplot(F1, [0, 10]);
a = 2;
f = @(t)exp(-a^2*t.^2/2);
w = 0: 0.4: 10;
for i = 0:4
  H = hankel_trans(f, w, i);
  line(w, H);
end
% v = 0与解析解曲线重合

05.05 z变换

z变换及其反变换
z变换并不是积分变换，只是因为其求解方法类似于其它积分变换，故归并到一起
在控制领域，z变换用于离散系统

z变换及反变换定义与性质
离散序列信号 $f(k), k = 1, 2, \cdots$
z变换定义：$\mathscr{Z}[f(k)] = \sum \limits_{k=0}^{\infty} f(k) z^{-k} = F(z)$
z变换性质：
线性性质：其中，a与b均为标量
$\mathscr{Z}[af(k) \pm bg(k)] = a\mathscr{Z}[f(k)] \pm b\mathscr{Z}[g(k)]$
时域后向平移性质：$\mathscr{Z}[f(k-n)] = z^{-n}F(z)$
前向平移性质: $\mathscr{Z}[f(k+n)] = z^n F(z) - \sum \limits_{i=0}^{n-1} z^{n-i} f(i)$
零初值：$\mathscr{Z}[f(k+n)] = z^n F(z)$
z域比例性质：$\mathscr{Z}[r^{-k}f(k)] = F(rz)$
频域微分性质：$\mathscr{Z}[kf(k)] = -z \dfrac{dF(z)}{dz}$
频域积分性质：$\mathscr{Z}[\dfrac{f(k)}{k}] = \int_0^{\infty} \dfrac{F(\omega)}{\omega} d \omega$
终止性质：
$\lim \limits_{k \to 0} f(k) = \lim \limits_{z \to \infty} F(z)$
$\lim \limits_{k \to \infty} f(k) = \lim \limits_{k \to 1} (z - 1) F(z)$
其中，F(z)无单位圆外的极点
卷积性质:$\mathscr{Z}[f(k)*g(k)] = \mathscr{Z}[f(k)]\mathscr{Z}[g(k)]$
式中离散信号的卷积算子*定义为
$f(k) * g(k) = \sum \limits_{l = 0}^{\infty} f(k) g(k-l)$
函数F(z) 的 z 反变换定义为
$f(k) = \mathscr{Z}^{-1}[f(k)] = \dfrac{1}{j2\pi} \oint F(z) z^{k-1} dz$

z变换求解
z变换求解类似于Laplace变换
函数调用格式：F = ztrans(fun) $\qquad$ F = ztrans(fun, k, z)
z反变换的函数调用格式
F = iztrans(fun) $\qquad$ F = iztrans(fun, z, k)

例5-23 给定函数求z变换
给定$f(kT) = akT - 2 + (akT + 2) e^{-akT}$
试求其z变换

syms a T k;
f = a*k*T - 2 + (a*k*T + 2)*exp(-a*k*t);
F = ztrans(f)
pretty(F)

例5-24 总结z反变换公式
给定 $F(z) = \dfrac{q}{(z^{-1} - p)^m}$
对不同的m值进行z反变换，总结一般规律

syms p q z;
assume(p ~= 0);
for i = 1:8
  disp(simplify(iztrans(q/(1/z-p)^i))),
end

一般规律：
$\mathscr{Z}^{-1} \Big[\dfrac{q}{(z^{-1} - p)^m} \Big] = \dfrac{(-1)^m q}{(m-1)!p^{n+m}} \prod \limits_{i=1}^{m-1}(n+i)$

双边z变换
k拓展到整个整数空间
$\mathscr{Z}[f(k)] = \sum \limits_{k=-\infty}^{\infty} f(k) z^{-k} = F(z)$
matlab没有现成函数
可以用底层求和命令直接分段求解
F = symsum(f*z^(-k), k, 0, inf) + symsum(f*z^(-k), k, -inf, -1)

例 5-25 双边z变换
分段函数 $f(n) = \left \{ \begin{array}{l}2^n, n \geq 0 \\ -3^n, n < 0 \end{array} \right.$

syms z n;
F = symsum(2^n*z^(-n), n, 0, inf) + symsum(-3^n*z^(-n), n, -inf, -1)

结果:$F = \dfrac{z}{z-2} + \dfrac{z}{z-3}$

有理函数z反变换的数值求解
有理函数的通解$F(z^{-1})$
$z^{-d} \dfrac{b_0 + b_1z^{-1} + b_2z^{-2} + \cdots + b_{m-1}z^{-(m-1)} + b_mz^{-m}}{a_0 + a_1z^{-1} + a_2z^{-2} + \cdots + a_{n-1}z^{-(n-1)} + a_nz^{-n}}$
$z^{-1}$ 的幂级数展开
$F(z^{-1}) = f_0 + f_1z^{-1} + f_2z^{-2} + \cdots = \sum \limits_{k=0}^{\infty} f_kz^{-k}$
自定义函数的调用: y = inv_z(num, den, d, N)

function y = inv_z(num, den, varargin)
  [d, N] = default_vals({0,10}, varargin{:});
  num(N) = 0;
  for i = 1:N-d
    k = num(1) / den(1);
    y(d+i) = k;
    if length(num) > 1
      ii = 2:length(den);
      if length(den) > length(num)
        num(length(den)) = 0;
      end
      num(ii) = num(ii)-k*den(ii);
      num(1) = [];
    end;
  end;

function varargout = default_vals(vals, varargin)
  if nargout ~= length(vals)
    error('number of arguments mismatch'); 
  else
    nn = length(varargin) + 1;
    varargout=varargin;
    for i = nn:nargout
      varargout{i} = vals{i};
    end;
  end;
end

例 5-26 数值z反变换
有理数 $G(z) = \dfrac{z^2 + 0.4}{z^5 - 4.1 z^4 + 6.71 z^3 - 5.481z^2 + 2.2356z - 0.3645}$
变换成
$F(z^{-1}) = z^{-3} \dfrac{1 + 0.4z^{-2}}{1 - 4.1z^{-1} + 6.71 z^{-2} - 5.481z^{-3} + 2.2356z^{-4} - 0.3645z^{-5}}$
求解与绘图

N = 50;
num = [1 0 0.4];
den = [1 -4.1 6.71 -5.481 2.2356 -0.3645];
y = inv_z(num, den, 3, N);
t = 0: (N-1);
stem(t, y)

05.06 复数映射与Riemann曲面

复数矩阵及其变换
函数调用格式(已知复数矩阵z)
共轭复数矩阵: $Z_1$ = conj(Z)
实部、虚部提取: R = real(Z) $\qquad$ I = imag(Z)
幅值、相位表示: A = abs(Z) $\qquad$ P = angle(Z)

复变函数映射

例5-28 已知复变函数 $f(z) = \dfrac{z^2 + 3z + 1}{(z-1)^5}$
求 $f^{(3)}(-j\sqrt{5})$

syms z;
f = (z^2 + 3*z + 4)/(z-1)^5;
f3 = diff(f, z, 3);
d3 = subs(f3, z, -sqrt(-5))

复函数的映射(变量替换)
平移: $z = \omega + \gamma$
反演: $z = 1/\omega$
双线性: $z = (a\omega + b)/(c \omega + d)$

例5-29 双线性变换映射
已知函数$f(z) = \dfrac{z^2 + 3z + 4}{(z-1)^5}$
双线性变换： $z = \dfrac{s - 1}{s + 1}$
直接映射:

syms z s;
f = (z^2 + 3*z + 4)/(z-1)^5;
F = simplify(subs(f, z, (s-1)/(s+1)))

例5-30 单位圆内样本点映射
左半平面点映射

[x, y] = meshgrid(-1: 0.1: 1);
ii = find(x.^2+y.^2 <= 1);
x = x(ii);
y = y(ii);
z = x + sqrt(-1)*y;
plot(z, '+');
hold on;
ezplot('x^2+y^2=1')
figure;
s = (z-1)./(z+1);
plot(s, 'x')

Riemann曲面绘制
复变函数映射图形绘制步骤
生成网格 z = cplxgrid(n)
计算数据：通过点运算计算函数值，如
f = z.^3.*sin(z.^2)
绘图(Riemann曲面): cplxmap(z, f)

function cplxmap1(z,w,B)
%CPLXMAP Plot a function of a complex variable.
%   CPLXMAP(z,f(z),(optional bound))
%   Used by CPLXDEMO.
%
%   See also CPLXGRID.

%   Copyright 1984-2002 The MathWorks, Inc. 
%   $Revision: 5.8 $  $Date: 2002/04/15 03:30:08 $
%   modified by Dingyu Xue for multi-valued complex functions
blue = 0.2;
x = real(z);
y = imag(z);
u = real(w);
v = imag(w);

if nargin > 2
   k = find((abs(w) > B) | isnan(abs(w)));
   if length(k) > 0
      u(k) = B*sign(u(k));
      v(k) = zeros(size(k));
      v = v/max(max(abs(v)));
      v(k) = NaN*ones(size(k));
   end
end

M = max(max(u));
m = min(min(u));
axis([-1 1 -1 1 m M]);
caxis([-1 1]);
s = ones(size(z));
surf(x,y,u,v);
colormap(hsv(64))

例5-31 Riemann曲面
复变函数 $f(z) = z^3 \sin z^2$

% 绘制Riemann曲面
z = cplxgrid(50);
f = z.^3.*sin(z.^2);
cplxmap(z, f)

多值 $f(z) = \sqrt[n]{z}$ 函数的Riemann曲面
方根函数的绘制 cplxroot(n)

例5-32 绘制Riemann面 $\sqrt[3]{z}, \sqrt[4]{z}$

cplxroot(3), figure, cplxroot(4)

该函数只能绘制方根函数，对其它多值复变函数无能为力

例5-32 重绘三次方跟曲面
三次方根函数的三个分支$\sqrt[3]{z}$
一个分支：$f_1(z) = \sqrt[3]{z}$
其余两个分支: $f_1(z) e^{-2j\pi/3}, f_1(z) e^{-4j \pi/3}$
直接绘制：

z = cplxgrid(30);
f1 = z.^(1/3);
a = exp(-2i*pi/3);
cplxmap1(z, f1), hold on;
cplxmap1(z, a*f1);
cplxmap1(z, a^2*f1);
zlim([-1, 1]), hold off;

优点：果能求出所有的分支，可以扩展到其它多值函数

05.07 奇点、极点与留数

留数的概念与计算
若函数f(z)在复平面的区域内各点处均为单值，且其导数为有限值，则称f(z)在复平面内为解析的
单支函数不解析的点称为奇点
使得f(z)分母多项式等于零的奇点又称为极点
如果 z = a 为f(z)函数的单奇点，则留数的定义为
$Res [f(z), z=a] = \lim \limits _{z \to a} (z-a) f(z)$
单极点留数: c = limit(F*(z-a), z, a)

重奇点的留数
若 z = a 为函数 f(z) 的 m 重奇点，则该点的留数定义为
$Res[f(z), z] = \lim \limits_{z \to a} \dfrac{1}{(m-1)!} \dfrac{d^{m-1}}{dz^{m-1}} [f(z)(z-a)^{m}]$
m重奇点的计算: c = limit(diff(F*(z-a)^m, z, m-1)/prod(1:m-1), z, a)

奇点与重数的计算
matlab调用格式: [p, m] = poles(f)
直接计算极点与留数 [r, p, m] = residuesym(f, a, b)

function [r, p, m] = residuesym(f, a, b)
  z = symvar(f);
  if nargin == 1
    [p, m] = poles(f);
  else
    [p, m] = poles(f, a, b);
  end;
  for k = 1: length(p)
    r(k) = limit(diff(f*(z-p(k))^m(k), z, m(k)-1) ...
       /factorial(m(k)-1), z, p(k));
  end;

例5-34 求奇点与留数
计算下式留数
$f(z) = \dfrac{1}{z^3(z-1)} sin(z + \dfrac{\pi}{3})e^{-2z}$
找出奇点z = 0, z = 1

syms z;
f = sin(z + pi/3)*exp(-2*z)/(z^3*(z-1))
F1 = limit(diff(f*z^3, z, 2)/prod(1:2), z, 0),
F2 = limit(f*(z-1), z, 1)
[r, p, m] = residuesym(f)

例5-35 奇点的重数
函数$f(z) = \dfrac{\sin z - z}{z^6}$
找出奇点 z = 0, 重数 m = 6?

syms z;
f = (sin(z) - z)/z^6;
limit(diff(f*z^6, z, 5)/prod(1:5), z, 0)

真正的重数
$\lim \limits_{z \to a} \dfrac{d^{k-1}}{dz^{k-1}} [(z-a)^k f(z)] < \infty$

例5-36 推导留数公式
函数$f(z) = \dfrac{1}{z \sin z}, 找出奇点 z = 0, z = \pm k \pi$
z = 0 点的留数

syms z;
f = 1/(z*sin(z));
c0 = limit(diff(f*z^2, z, 1), z, 0)

$z = \pm k \pi$ 留数计算, residuesym函数不能使用

syms k;
assume(k, 'integer'),
assumeAlso(k ~= 0);
R = simplify(limit(f*(z-k*pi), z, k*pi))

05.08 部分分式展开与封闭曲线积分计算

有理函数的部分分式展开
有理函数
$G(x) = \dfrac{B(x)}{A(x)} = \dfrac{b_1x^m + b_2x^{m-1} + \cdots + b_mx + b_{m+1}}{x^n + a_1x^{n-1} + a_2x^{n-2} + \cdots + a_{n-1}x + a^n}$
其中，$a_i$ 和 $b_i$ 均为常数
分子分母多项式最大公约式(GCD)，判定互质
C = gcd(A, B)
得出互质多项式，约简 A/C, B/C

例5-37 多项式互质
$A(x) = x^4 + 7x^3 + 13x^2 + 19x + 20$
$B(x) = x^7 + 16x^6 + 103x^5 + 346x^4 + 655x^3 + 700x^2 + 393x + 90$
判定它们是否互质

syms x;
A = x^4 + 7*x^3 + 13*x^2 + 19*x +20;
B = x^7 + 16*x^6 +103*x^5 + 346*x^4 + 655*x^3 + 700*x^2 + 393*x + 90;
d = gcd(A, B)
% 分式化简
simplify(A/d), simplify(B/d)

部分分式展开
部分分式展开
假设 A(x) 和 B(x) 互质
并且 A(x) = 0的根都为单根 $-p_i, i = 1, 2, \cdots, n$
则部分分式展开为
$G(x) = \dfrac{r_1}{x + p_1} + \dfrac{r_2}{x + p_2} + \cdots + \dfrac{r_n}{x + p_n}$
留数的计算如下：
$r_i = Res[G(-p_i)] = \lim \limits_{x \to -p_i} G(s) (x + p_i)$

有重根时的部分分式展开
假设 $-p_i$ 是k重根
$G(x) = \dfrac{r_i}{x + p_i} + \dfrac{r_{i+1}}{(x+p_i)^2} + \cdots + \dfrac{r_{i+k-1}}{(x+p_i)^{k}}$
各个系数如下: j = 1, 2, …, k
$r_{i+j-1} = \dfrac{1}{(k-1)!} \lim \limits_{x \to -p_i} \dfrac{d^{j-1}}{dx^{j-1}} [G(x) (x+p_i)^k]$
部分分式展开的函数–数值计算
[r, p, K] = residue(b, a)

部分分数展开的符号运算
早期版本编写的重载函数 residue
新版本，不支持重载函数
MuPAD内核 partfrac函数，底层命令
F = feval(symengine, ‘partfrac’, f)
编写一个接口
F = partfrac(f) or F = partfrac(f, s)

function Y = partfrac(varargin)
  Y = feval(symengine, 'partfrac', varargin{:});

例5-38 解析方法的展开
部分分式展开
$G(s) = \dfrac{s^3 + 2s^2 + 3s + 4}{s^6+11s^5+48s^4+106s^3 + 125s^2 + 75s + 18}$

syms s;
f = (s^3 + 2*s^2 + 3*s + 4)/(s^6 + 11*s^5 + 48*s^4 + 106*s^3 + 125*s^2 + 75*s+18);
G1 = partfrac(f)
pretty(G1)

封闭曲线积分问题计算
封闭曲线积分的数学表达式 $\oint_{\Gamma} f(z) dz$
其中，$\Gamma$ 是一个逆时针方向的闭曲线，那么
$\oint_{Gamma} f(z) dz = j 2 \pi \sum \limits_{i=1}^{m} Res [f(p_i)]$
该封闭曲线包围m个奇点
$p_i, (i = 1, 2, \cdots, m)$

例5-34 封闭曲线积分
$f(z) = \dfrac{2z^7 + 2z^3 + 8}{z^8 + 30z^7 + 386 z^6 + 2772z^5 + 12093z^4 + 32598z^3 + 52520z^2 + 45600z + 16000}$
计算在 |z| = 6 上的封闭曲线积分
$\dfrac{13041}{(s+5)^3} + \dfrac{341863}{12(s+5)^2}+\dfrac{7198933}{144(s+5)} - \dfrac{16444}{3(s+4)^3} + \dfrac{193046}{9(s+4)^2} - \dfrac{1349779}{27(s+4)} + \dfrac{11}{9(s+2)} + \dfrac{1}{432(s+1)}$
$\oint_{|z|=6} f(z) dz = 2 \pi j [\dfrac{7198933}{144} - \dfrac{1349779}{27} + \dfrac{11}{9} + \dfrac{1}{432}] = 4\pi j$

直接积分方法
圆的方程(局限性)
$z = 6 \cos t + j 6 \sin t, t \in (0, 2 \pi)$
封闭曲线积分直接计算

syms z t;
G = (2*z^7 + 2*z^3 + 8) ...
  /(z^8 + 30*z^7 + 386*z^6 + 2772*z^5 + 12093*z^4 + ...
  32598*z^3 + 52520*z^2 + 45600*z + 16000);
F = subs(G, z, 6*cos(t) + 6*sin(t)*sqrt(-1));
I = int(F*diff(6*cos(t)+6*sin(t)*sqrt(-1), t), t, 0, 2*pi)

例5.44 封闭曲线积分
$\oint_{\Gamma} \dfrac{1}{(z+j1)^{10}(z-1)(z-3)} dz$
其中$|Gamma}$ 是|z| = 2(逆时针圆封闭曲线)
曲线内奇点 z = 1, z = -j1

i = sym(sqrt(-1));
syms z;
f = 1/((z+i)^10*(z-1)*(z-3));
r1 = limit(diff(f*(z+i)^10, z, 9)/prod(1:9), z, -i);
r2 = limit(f*(z-1), z, 1);
a = 2*pi*i*(r1 + r2)

05.09 差分方程求解

差分方程求解
常系数线性差分方程的一般形式为:
$y[(k+n)T] + a_1y[(k+n-1)T] + \cdots + a_ny(kT)$
$\qquad = b_1u[(k-d)T] + \cdots + b_mu((k-d-m+1) T]$
差分方程简单记号:
$y(t+n) + a_1 y(t+n-1) + \cdots + a_ny(t)$
$\qquad = b_1u(t+m-d) + \cdots + b_{m+1}u(t-d)$
T采样周期，常系数

一般差分方程的解析解求解方法
差分方程求z变换
$z^nY(z) - \sum \limits_{i=0}^{n-1} z^{n-i}y(i) + a_1z^{n-1}Y(z) -a_1 \sum \limits_{i=0}^{n-2} z^{n-i}y(i) + \cdots + a_nY(z)$
$=z^{-d}[b_1z^mU(z) - b_1\sum limits_{i=0}^{m-1} z^{n-i}u(i) + \cdots + b_{m+1}U(z)]$
得出：
$Y(z) = \dfrac{(b_1z^m + b_2z^{m-1} + \cdots + b_{m+1})z^{-d}U(z) + E(z)}{z^n + a_1z^{n-1} + z^2z^{n-2} + \cdots + a_n}$
其中
$E(z) = \sum \limits_{i=0}^{n-1} z^{n-i} y(i) - a_1 \sum \limits_{i=0}^{n-2} z^{n-i} y(i) - a_2 \sum \limits_{i=0}^{n-3} z^{n-i}y(i) - \cdots - a_{n-z}y(0) + \hat u(n)$
$\hat u(n) = -b_1 \sum \limits_{i=0}^{m-1} z^{n-i}u(i) - \cdots - b_m zu(0)$

算法实现

function y = diff_eq(A, B, y0, U, d)
  E = 0;
  n = length(A) - 1;
  syms z;
  if nargin == 4
    d=0;
  end
  m = length(B)-1;
  u = iztrans(U);
  u0 = subs(u, 0:m-1);
  for i = 1:n
    E = E + A(i)*y0(1:n+1-i)*[z.^(n+1-i:-1:1)].';
  end
  for i = 1:m
    E = E-B(i)*u0(1:m+1-i)*[z.^(m+1-i:-1:1)].';
  end
  Y = (poly2sym(B,z)*U*z^(-d)+E)/poly2sym(A,z);
  y = iztrans(Y); 

求解语法：y = diff_eq(A, B, y0, U, d)

例5-45 差分方程数值解
差分方程
$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$

syms z n;
u = (1/5)^n;
U = ztrans(u);
y = diff_eq([48 -76 44 -11 1], [2 3 1], [1 2 0 -1], U)
n0 = 0:20;
y0 = subs(y, n, n0);
stem(n0, y0)

线性时变差分方程的数值解法
线性时分差分状态方程
$\left \{ \begin{array}{l}x(k+1) = F(k)x(k) + G(k)u(k) \\ y(k) = C(k)x(k) + D(k)u(k) \end{array} \right.$
其中, $x(0) = x_0$
递推:
$x(1) = F(0)x_0 + G(0)u(0)$
$x(2) = F(1)x(1) +G(1)u(1)$
$\qquad = F(1)F(0)x_0 + F(1)G(0)u(0) + G(1)u(1)$
$x(k) = F(k-1)F(k-2) \cdots F(0) x_0 + G(k-1)u(k-1) + F(k-1)G(k-2)u(k-2) + \cdots + F(k-1) \cdots F(0)G(0)u(0)$
$\qquad = \prod \limits_{j=0}^{k=1} F(j) x_0 + \sum \limits_{i=0}^{k-1} \Big[ \prod \limits_{j=i+1}^{k-1} F(j) \Big] G(i)u(i)$
可通过递推方法求解

例5-46 时变差分方程
求离散线性时变差分方程
$\begin{bmatrix}x_1(k+1) \\ x_2(k+1) \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ 1 & \cos(k \pi) \end{bmatrix} \begin{bmatrix} x_1(k) \\ x_2(k) \end{bmatrix} + \begin{bmatrix} \sin(k \pi /2) \\ 1 \end{bmatrix} u(k)$
其中
$\begin{bmatrix}x_1(0) \\ x_2(0) \end{bmatrix} = \begin{bmatrix} 1 \\ 1 \end{bmatrix}, and u(k) = \left \{ \begin{array}{l} 1, k = 0, 2, 4, \cdots \\ -1, k = 1, 3, 5, \cdots \end{array} \right.$
可以通过递推方法直接求解

x0 = [1; 1];
x = x0;
for k = 0: 100
  if rem(k, 2) == 0
    u = 1;
  else
    u = -1;
  end;
  F = [0 1; 1 cos(k*pi)];
  G = [sin(k*pi/2); 1];
  x1 = F*x0 + G*u;
  x0 = x1;
  x = [x x1];
end;
subplot(211), stairs(x(1, :)),
subplot(212), stairs(x(2, :))

线性时不变差分方程的解法
若线性时不变系统有
$F(k) = \cdots = F(0) = F, G(k) = \cdots = G(0) = G$
则有
$x(k) = F^kx_0 + \sum limits_{i=0}^{k-1} F^{k-i-1} Gu(i)$
线性时不变状态方程
$\left \{ \begin{array}{l} x(k+1) = Fx(k) + Gu(k) \\ y(k) = Cx(k) + Du(k) \end{array} \right. , x(0) = x_0$

直接求解
两边取z变换
$X(z) = (zI -F)^{-1}[zx_0 + GU(z) - Gzu_0]$
由此得出解析解
$x(k) = \mathscr{Z}^{-1}[(zI-F)^{-1}z]x_0 + \mathscr{Z}^{-1}\lbrace(zI -F^{-1})[GU(z) - Gzu_0]\rbrace$
由iztrans()函数可以直接求解

例5-47 差分方程解析解
离散系统的状态方程如下
$x(k+1) = \begin{bmatrix}11/6 & -5/4 & 3/4 & -1/3 \\ 1 & 0 & 0 & 0 \\ 0 & 1/2 & 0 & 0 \\ 0 & 0 & 1/4 & 0 \end{bmatrix} x(k) + \begin{bmatrix}4 \\ 0 \\ 0 \\ 0 \end{bmatrix} u(k)$
零初始状态: $x_0 = 0$

F = sym([11/6, -5/4, 3/4, -1/3; 1, 0, 0, 0; 0, 1/2, 0, 0; 0, 0, 1/4, 0]);
G = sym([4; 0; 0; 0]);
syms x k;
U = ztrans(sym(1));
x = iztrans(inv(z*eye(4)-F)*G*U, z, k)

一般非线性离散系统的求解方法
假设已知差分方程的显示形式，即
$y(t) = f(t, y(t-1), \cdots, y(t-n), u(t), \cdots, u(t-m))$

例5-46 非线性系统
$y(t) = \dfrac{y(t-1)^2 + 1.1y(t-2)}{1 + y(t-1)^2 + 0.2y(t-2) + 0.4y(t-3)} + 0.1u(t)$
输入信号: $u(t) = \sin(t)$
采样周期: T = 0.05

y0 = zeros(1,3);
T = 0.05;
t = 0: T: 4*pi;
u = sin(t);
for i = 1: length(t)
  y(i) = (y0(3)^2 + 1.1* y0(2))/(1+y0(3)^2+...
    0.2*y0(2) + 0.4*y0(1)) + 0.1*u(i);
  y0 = [y0(2:3), y(i)];
end;
plot(t, y, t, u)

非线性系统的输出会产生畸变