matlab/Matlab多项式运算与数据可视化

Matlab多项式表示与运算

MATLAB习惯将降幂排列的多项式 $\boldsymbol{a}(\boldsymbol{x})=\boldsymbol{a}_{1} \boldsymbol{x}^{\boldsymbol{n}}+\boldsymbol{a}_{2} \boldsymbol{x}^{\boldsymbol{n}-\mathbf{1}}+$$\cdots+a_{n} x+a_{n+1}$ 存储为系数行向量 $a=\left[a_{1}, a_{2}, \ldots, a_{n+1}\right]$

  • 多项式的乘法,基于MATLAB的多项式系数定义,有:

注意:乘积多项式的系数恰为两个原始系数的卷积

Matlab函数

  • c=conv (a,b):一种快速获取多项式相乘系数的方法
  • [q,r]=deconv(b,a) :多项式带余除法, q为商,r为余式
  • [r,p,k]=residue(b,a):表示 $\frac{b(x)}{a(x)}=\sum_{i=1}^{n} \frac{r_{i}}{x-p_{i}}+\boldsymbol{k}(\boldsymbol{x})$
    其中r为留数向量,p为极点向量,k为余式系数行向量

  • r=roots(a):求多项式a(x)的所有根(含重根和任意复根),注意到该函数返回的根用列向量存储。

  • a=poly(r):当r为一个行向量时,a为根等于r的多项式(即函数roots的逆运算);而当r为方阵时,poly函数将获得r的特征多项式(用于计算特征值)

  • V=polyval(p,X)p为一个多项式的系数行向量,X可以是任意规模的矩阵,此函数将X每一个元素代入p(x)并获得代入后的数值矩阵p(X)

  • V=polyvalm(p,X):其含义略有不同,此函数意味计算以矩阵X为整体的多项式代入的结果矩阵,即

  • poly2str(a,'s'):以a向量元素为系数,s为自变量生成多项式字符串,按降幂顺序排列。

实例应用

例:多项式带余除法

1
2
3
4
5
6
7
8
9
format rat
p1 = conv([1, 0, 2], conv([1, 4], [1, 1]));
p2 = [1, 0, 1, 1];
[q, r] = deconv(p1, p2);
cq='商多项式为 '; cr='余多项式为 ';
disp([cq,poly2str(q,'s')]),disp([cr,poly2str(r,'s')])
qp2=conv(q,p2); %回代pp1=商X分母+余式
pp1=qp2+r;
pp1==p1

例:矩阵的特征多项式与特征根

1
2
3
4
5
6
7
8
9
10
format short
A=[11 12 13;14 15 16;17 18 19];
PA=poly(A);
PPA=poly2str(PA,'s');
s = eig(A);
r = roots(PA);
n = length(PA);
AA = diag(ones(1, n-2, class(PA)), -1);
AA(1,:) = -PA(2:n) ./ PA(1); %第一行减掉了特征多项式后三个系数
sr = eig(AA); %该矩阵AA的特征多项式与PA相同,特征值也相同

我们可以看到PA的结果为1.0000 -45.0000 -18.0000 -0.0000,这是该矩阵的特征多项式,即$\lambda^3-45\lambda^{2}-18\lambda = 0$,打印PPA显示特征多项式的字符串结果:s^3 - 45 s^2 - 18 s - 9.8142e-15',最后一项为机器误差可以忽略。由高等代数的知识我们知道该特征多项式的根即为eig的打印结果,因此我们可以发现sr是相同的。

AA为用长度为n-2的向量生成一个向下平移一行的(n-1)×(n-1)维矩阵,打印结果为:

1
2
3
0   0   0
1 0 0
0 1 0

例:构造特征根的多项式

1
2
3
4
R=[-0.5,-0.3+0.4*i,-0.3-0.4*i];	%三个特征根,复根需有共轭复数对
P=poly(R) %生成对应的多项式
PR=real(P)
PPR=poly2str(PR,'x')

PPR的打印结果为:

1
2
PPR =
x^3 + 1.1 x^2 + 0.55 x + 0.125

例:polyvalpolyvalm函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
clear
p=[1,2,3];
poly2str(p,'x')
X=[1,2;3,4];
va = X.^2+2*X+3 %对应元素代入多项式p(x)
Va = polyval(p,X) %等价的MATLAB函数polyval

vm = X^2+2*X+3*eye(2) %将矩阵整体代入p(x)
Vm = polyvalm(p,X) %等价的MATLAB函数

cp=poly(X);
poly2str(cp,'x') %获取并显示X的特征多项式
cpXa = polyval(cp,X) %代入X的元素,得到结果矩阵
cpX = polyvalm(cp,X) %整体代入验证H-C定理

MATLAB有限长向量卷积函数conv

已知 $a=\left[a_{1}, a_{2}, \ldots, a_{n}\right], b=\left[b_{1}, b_{2}, \ldots, b_{m}\right],$ 类似于前面多项式乘法的例子, $\operatorname{conv}(a, b)$ 函数将返回 $a * b$ 的向量, 即

这里补充定义 $\forall i>n, a_{i}=\mathbf{0}, \forall j>m, b_{j}=0$

  • conv(a,b,'same')将返回一个长度等于 $n$ 的卷积结果, 即截取上述卷积结果中间n个元素的结果, 常用于信号或图像处理问题。若3个选2个,则默认选择后两个。
  • conv(a,b,'valid')将返回一个长度等于 $n-m+1$ 的卷积 结果, 即仅包含 $b$ 的全部元素参与卷积的中间项。
  • deconv为反卷积函数。卷积运算之前对于边界条件的处理有时需要跟实际问题相结合。具体方法有0填充(默认),对称延拓,复制延拓,周期延拓等。

例:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
% 方法1:
N1=3;N2=12;a=ones(1,N2+1);a(1:N1)=0;%相同的向量定义,用0填充左侧
M1=2;M2=9;b=ones(1,M2+1);b(1:M1)=0;
c=conv(a,b); %直接conv函数计算
kc=0:(N2+M2);

% 方法2
N1=3;N2=12;M1=2;M2=9;
A=ones(1,(N2-N1+1)); B=ones(1,(M2-M1+1)); %去掉左边的0,与前法类似
C=conv(A,B);
Nc1=N1+M1;Nc2=N2+M2; %左右端点下标计算
KC=Nc1:Nc2;

subplot(2,1,1),stem(kc,c), text(20,6,'0 起点法') %方法1直接绘图
CC=[zeros(1,KC(1)),C]; %方法2的结果左侧填充5个0
subplot(2,1,2),stem(kc,CC),text(18,6,'非平凡区间法')
xlabel('n') %注:计算单点的卷积值使用sum函数速度更快

  • 版权声明: 本博客所有文章除特别声明外,著作权归作者所有。转载请注明出处!

扫一扫,分享到微信

微信分享二维码
  • Copyrights © 2020-2021 chenk
  • 由 帅气的CK本尊 强力驱动
  • 访问人数: | 浏览次数:

请我喝杯咖啡吧~

支付宝
微信