1. 程式人生 > >matlab多項式及其運算

matlab多項式及其運算

0 建立多項式

多項式的一般形式如下:

f(x) = a_{0}x^{n} + a_{1}x^{n-1} + a_{2}x^{n-2} + ... + a_{n-1}x + a_{n}

我們可以使用它的係數向量來表示,P = [a_{0}, a_{1}, ... , a_{n-1}, a_{n}]

matlab中,提供了poly2sym函式實現多項式的構造。

  • r = poly2sym(c):c為多項式的係數向量
  • r = poly2sym(c, v):c為多項式的係數向量,v為其變數
>> poly2sym([1 3 2])
 
ans =
 
x^2 + 3*x + 2

>> poly2sym(sym([1 0 1 -1 2]),sym('y'))
 
ans =
 
y^4 + y^2 - y + 2

1.多項式的求根

多項式的根:

matlab使用roots函式求解多項式的根,即求解函式等於0的根

  • r = roots(c):其中c為多項式的係數向量,r為求解多項式的根
>> p = [1 -12 0 25 116];
>> r = roots(p);
>> p

p =

     1   -12     0    25   116

>> r

r =

  11.7473 + 0.0000i
   2.7028 + 0.0000i
  -1.2251 + 1.4672i
  -1.2251 - 1.4672i

由根建立多項式:

matlab中規定,多項式是行向量,根是列向量。給出一個多項式的根,也可以構造相應的多項式。

  • p = poly(A):如果A為方陣,則多項式p為該方陣的特徵多項式;如果A為向量,則A的元素為該多項式p的根。n階方陣的特徵多項式存放在行向量中,並且特徵多項式最高次的係數一定為1
>> pp = poly(r)

pp =

    1.0000  -12.0000   -0.0000   25.0000  116.0000

2.多項式的四則運算

多項式的加法:

如果兩個多項式向量大小相同,相加時就與標準的陣列加法相同

>> p1 = [5 40 6 21 9 3];
>> p2 = [4 0 3 72 1 8];
>> p3 = p1 + p2

p3 =

     9    40     9    93    10    11

>> r1 = poly2str(p3,'x')   % 顯示多項式

r1 =

   9 x^5 + 40 x^4 + 9 x^3 + 93 x^2 + 10 x + 11

>> p4 = p1 - p2

p4 =

     1    40     3   -51     8    -5

>> r2 = poly2str(p4,'x')   % 顯示多項式

r2 =

   x^5 + 40 x^4 + 3 x^3 - 51 x^2 + 8 x - 5

注意:當兩個多項式階次不同時,低階的多項式用首0填補,使其與高階多項式有同樣的階次。要求首零而不是尾0,是因為相關的係數像x冪一樣,必須整齊 

多項式的乘法:

conv函式實現多項式的乘運算,deconv函式實現多項式的除運算

  • c = conv(a,b):執行a,b兩個向量的卷積運算
  • c = conv(a,b,'shape'):按形參‘shape’返回卷積運算,shape取值如下:
    • full:為返回完整的卷積,是預設值
    • same:為返回部分卷積,其大小與向量a大小相同
    • valid:只返回無填充0部分的卷積,此時輸出向量c的最大值為max(length(a) - max(0,length(b) - 1),0).
>> f = [1 4 -2 7 11];
>> g = [9 -11 5 0 8];
>> c = conv(f,g)

c =

     9    25   -57   105    20   -54    39    56    88

注意:conv函式只能進行兩個多項式的乘法,兩個以上的多項式的乘法需要重複使用conv

  • [q,r] = deconv(v,u):求多項式v, u的除法運算,其中q為返回多項式v除以u的商式,r為返回v除以u的餘式。
>> c = [1 5 15 35 69 100 118 110 72];
>> b = [1 2 3 6 8];
>> [a,r] = deconv(c,b)

a =

     1     3     6     8     9


r =

     0     0     0     0     0    -2    -5    -8     0

3.多項式的導數

  • k = polyder(p):求多項式的導函式多項式
  • k = polyder(a,b):求多項式a與多項式b乘積的導函式多項式
  • [q,d] = polyder(b,a):求多項式b與多項式a相除的導函式,導函式的分子存入q,分母存入d
>> a = [3 6 9];
>> b = [1 2 0];
>> k = polyder(a,b)

k =

    12    36    42    18

>> K = poly2str(k,'x')

K =

   12 x^3 + 36 x^2 + 42 x + 18
>> [q,d] = polyder(b,a)

q =

    18    18


d =

     9    36    90   108    81

4.多項式的積分

  • polyint(p,k):返回以向量p為係數的多項式積分,積分的常數項為k
  • polyint(p):返回以向量p為係數多項式的積分,積分的常數項為預設值0
>> p = [1 -1 2];
>> k = 1/2;
>> F = polyint(p,k)

F =

    0.3333   -0.5000    2.0000    0.5000

>> df = poly2sym(F)
 
df =
 
x^3/3 - x^2/2 + 2*x + 1/2

5.多項式的估值

matlab提供了polyval函式與polyvalm函式用於求多項式p(x)在x=a的取值。輸入可以是標量或矩陣

  • y = polyval(p,x):p為多項式的係數向量,x為矩陣,它是按陣列運算規則來求多項式的值
  • [y,delta] = polyval(p,x,S):使用可選的結構陣列S產生由polyfit函式輸出的估計引數值;delta是預測未來的觀測估算的誤差標準偏差
  • y = polyval(p,x,[],mu)或[y,delta] = polyval(p,x,S,mu):使\hat{x} = (x - \mu _{1}) / \mu _{2}替代x,\mu _{1} = mean(x), \mu _{2} = std(x),其中心點與座標值mu = [\mu _{1},\mu _{2}]可由polyfit函式計算得出

polyvalm函式的輸入引數只能是N階方陣,這時可以將多項式看作矩陣函式

  • Y = polyvalm(p,X):p為多項式的係數向量,X為方陣,其實按矩陣運算規則來求多項式的值。
>> X = pascal(4)

X =

     1     1     1     1
     1     2     3     4
     1     3     6    10
     1     4    10    20

>> p = poly(X)

p =

    1.0000  -29.0000   72.0000  -29.0000    1.0000

>> P = poly2str(p,'x')

P =

   x^4 - 29 x^3 + 72 x^2 - 29 x + 1

>> y = polyval(p,X)

y =

   1.0e+04 *

    0.0016    0.0016    0.0016    0.0016
    0.0016    0.0015   -0.0140   -0.0563
    0.0016   -0.0140   -0.2549   -1.2089
    0.0016   -0.0563   -1.2089   -4.3779

>> y = polyvalm(p,X)

y =

   1.0e-10 *

   -0.0003   -0.0036   -0.0052   -0.0143
   -0.0021   -0.0136   -0.0179   -0.0464
   -0.0059   -0.0330   -0.0400   -0.1047
   -0.0130   -0.0639   -0.0750   -0.1962

6.有理多項式

matlab中,有理多項式由它們的分子多項式和分母多項式表示。對有理多項式進行運算的兩個函式是residue和polyder。redidue執行部分分式展開的運算

  • [r,p,k] = residue(b,a):b,a分別為分子和分母多項式係數的行向量,r為留數行向量
  • [b,a] = residue(r,p,k):p為極點行向量,k為直項行向量
>> b = [5 3 -2 7];
>> a = [-4 0 8 3];
>> [r,p,k] = residue(b,a)

r =

   -1.4167
   -0.6653
    1.3320


p =

    1.5737
   -1.1644
   -0.4093


k =

   -1.2500

>> [b,a] = residue(r,p,k)

b =

   -1.2500   -0.7500    0.5000   -1.7500


a =

    1.0000   -0.0000   -2.0000   -0.7500

7.多項式的微分

  • k = polyder(p):p,k分別為原多項式及微分多項式的多項式表示
  • k = polyder(a,b):求多項式a與多項式b乘積的導函式多項式
  • [q,b] = polyder(b,a):求多項式b與多項式a相除的導函式,導函式的分子存入q,分母存入d
>> a = [3 6 9];
>> b = [1 2 0];
>> k = polyder(a,b)

k =

    12    36    42    18

>> K = poly2str(k,'x')

K =

   12 x^3 + 36 x^2 + 42 x + 18
>> [q,d] = polyder(b,a)

q =

    18    18


d =

     9    36    90   108    81