MATLAB 程式設計:進階篇 多項式的處理與分析 張智星 (Roger Jang) jang@mirlab.org http://mirlab.org/jang 台大資工系 多媒體檢索實驗室
7-1多項式的加減乘除 n 次多項式的表示法 MATLAB 中,用 n+1 的列向量表示 p(x) Ex. p = [1,2,3,1] 表示一個三次多項式
多項式的加減 (I) 計算其和與差 多項式的加減,直接由向量的加減推出 兩多項式, 及 範例7-1:polyPlus.m 兩多項式, 及 計算其和與差 範例7-1:polyPlus.m p1 = [1,0,1,1]; p2 = [0,1,-1,2]; p1 + p2 p1 - p2
多項式的加減 (II) 矩陣 p1 與 p2 的長度要一致,否則 MATLAB 會產生運算錯誤的訊息 ans = 1 1 0 3 1 1 0 3 1 -1 2 -1 矩陣 p1 與 p2 的長度要一致,否則 MATLAB 會產生運算錯誤的訊息
多項式的乘除 (I) 多項式的乘除,使用 conv 及 deconv 指令達成 欲求多項式 與 的乘積 以多項式的表示法,為 欲求多項式 與 的乘積 範例7-2:conv01.m 以多項式的表示法,為 p1 = [1, 0, 1, 1]; p2 = [0, 1, -1, 2]; p3 = conv(p1, p2) p3 = 0 1 -1 3 0 1 2
多項式的乘除 (II) 欲求 除以 的商式與餘式 代表 除以 ,得到商式 ,餘式 範例7-3:deconv01.m 欲求 除以 的商式與餘式 範例7-3:deconv01.m 代表 除以 ,得到商式 ,餘式 p1 = [1, 0, 1, 1]; p2 = [1, -1, 2]; [q, r] = deconv(p1, p2) q = 1 1 r = 0 0 0 -1
多項式的加減乘除 函數 說明 p1 + p2 多項式 與 的和 p1 - p2 多項式 與 的差 conv(p1, p2) 多項式 與 的和 p1 - p2 多項式 與 的差 conv(p1, p2) 多項式 與 的乘積 [q, r] = deconv(p1, p2) 多項式 除以 ,得到商式為 ,餘式為
7-2 多項式的求值、求根 、微分與積分 計算多項式的值,使用 polyval 指令 範例7-4:polyval01.m x = 0:0.1:3; y = polyval(p, x); plot(x, y, '-o');
計算多項式的值-1 x 和 y 是長度為 31 的向量,y(i) 的值為 在 x = x(i) 的函數值
計算多項式的值-2 計算p(A),A為方陣,可用polyvalm指令 範例7-5:polyvalm01.m 結果和 B = A^2 + 2*A + 1 相同 p = [1 2 1]; A = [1 2; 3 4]; B = polyvalm(p, A) B = 10 14 21 31
計算多項式的根 (I) 多項式的根,可用 roots 指令求得 多項式 的根 範例7-6:roots01.m 多項式 的根 範例7-6:roots01.m p = [1, 3, 1, 5, -1]; % 多項式 r = roots(p); % 求多項式的根 r = -3.2051 0.0082 + 1.2862i 0.0082 - 1.2862i 0.1886
計算多項式的根 (II) 為驗證此四根,可輸入如下 範例7-7:roots02.m p = [1, 3, 1, 5, -1]; % 多項式 r = roots(p); % 求多項式的根 polyval(p, r); % 將根帶入多項式求值 ans = 1.0e-012 * 0.1465 -0.0135 + 0.0271i -0.0135 - 0.0271i 0.0002
提示 fzero指令: Roots指令: 用於一般函數的求根 一次只找一個根 使用牛頓法 只用於多項式求根 一次找到全部的根 將多項式表示成「伴隨矩陣」(Companion Matrix),再用解特徵值方法求根
多項式的微分 求多項式的微分,使用指令polyder 範例7-8:polyder01.m 表示 微分後的結果為 p = [1 3 3 1]; 表示 微分後的結果為 p = [1 3 3 1]; q = polyder(p) q = 3 6 3
多項式的積分 (I) MATLAB 6.x 以後已經提供 polyint 指令 範例7-9:polyint02.m 表示 積分後為 表示 積分後為 p = [4 3 2 1]; k = 8; % 積分後的不定常數 q = polyint(p, k) % 積分後的多項式 q = 1 1 1 1 8
多項式的積分 (II) MATLAB 5.x 並無對多項式積分的指令 用其它方法達成積分 範例7-10:polyint01.m t = length(p):-1:1; k = 8; % 積分後的不定常數 q = [p./t, k] % 積分後的多項式 q = 1 1 1 1 8
多項式的值、根、積分、微分 函數 說明 q = polyval(p, x) 計算 p(x) 的值。 q = polyvalm(a, A) 計算 p(A),A 為一方陣。 q = polyder(p) q(x) 為 p(x) 的微分。 q = [p./length(p):-1:1, k] q(x) 為 p(x) 的積分,其中 k 為任意常數。 q = polyint(p, k) q(x) 為 p(x) 的積分,其中 k 為任意常數。(本指令只適用於 MATLAB 6.x 以後的版本。)
矩陣的特徵多項式 給定一方陣 A,其特徵多項式為 使用 poly 指令計算特徵多項式 即為矩陣 A 的特徵多項式 範例7-11:poly01.m 即為矩陣 A 的特徵多項式 A = [1 3 4; 2 4 1; 1 6 2]; p = poly(A) p = 1.0000 -7.0000 -2.0000 -25.0000
矩陣的特徵值 特徵多項式的根即為矩陣 A 的特徵值 範例7-12:poly02.m A = [1 3 4; 2 4 1; 1 6 2]; p = poly(A); % 方陣的特徵多項式 r = roots(p); % 特徵方程式的根,亦即固有值 det(A-r(1)*eye(3)) det(A-r(2)*eye(3)) det(A-r(3)*eye(3)) ans = -9.5707e-014 ans = -2.5292e-014 - 1.5899e-015i ans = -2.5292e-014 +1.5899e-015i
提示 特徵方程式的根會讓 ,也是方陣 A 的特徵值或固有值(Eigenvalues)的定義 可用eig指令來計算方陣 A 的特徵值及特徵向量
7-4 部份分式展開 指令 residue 計算一分式的部份分式展開 A(s) 和 B(s) 為多項式,且 B(s) 無重根 分式 可以表示為 , ,…, 為 的根(或 的極點) , ,…, 為常數 為一多項式
多項式部份分式展開 (I) 欲求 的部份分式展開 範例7-13:residue01.m b = [3 8]; a = [1 5 6]; 欲求 的部份分式展開 範例7-13:residue01.m b = [3 8]; a = [1 5 6]; [r, p, k] = residue(b, a)
多項式部份分式展開 (II) 由以結果得知, 。 r = 1.0000 2.0000 p = -3.0000 -2.0000 k = []
提示 部份分式展開特別適用於線性系統之轉換函數(Transfer Function)的分析 若 為一系統之轉換函數,經由上述部份分式展開後,可知其脈衝響應(Impulse Response)為
多項式部份分式還原 residue 指令將部份分式轉回原來的形式 範例7-14:residue02.m 將部份分式結合成單一分式: p = [-3 -2]; k = []; [b2, a2] = residue(r, p, k) b2 = 3 8 a2 = 1 5 6
7-5多項式擬合 曲線擬合(Curve Fitting) 資料分析上常用的方法 以一參數化的曲線逼近一組給定的資料點 若參數化的曲線是一多項式,此曲線擬合稱為多項式擬合(Polynomial Fitting) polyfit 指令用以進行多項式擬合
多項式擬合範例 (I) 美國人口在1790至1990(10年一期)資料 範例7-15:plotCensus.m load census.mat plot(cdate, pop, 'o'); xlabel('年份'); ylabel('人口(單位:百萬)');
多項式擬合範例 (II)
多項式擬合範例 (III) 欲預測美國在西元 2002 年的總人口,使用 polyfit 找出多項式逼近已知的資料點 範例7-16:polyfit01.m load census.mat plot(cdate, pop, 'o'); p3 = polyfit(cdate, pop, 3); cdate2 = 1790:2002; pop2 = polyval(p3, cdate2); plot(cdate, pop, 'o', cdate2, pop2, '-'); xlabel('年份'); ylabel('人口(單位:百萬)'); legend('實際人口', '預測人口'); popAt2002 = polyval(p3, 2002)
多項式擬合範例 (IV) 此曲線預測在 2002 年美國人口約為 282.5 百萬人 population = 282.5332
不同次數的比較 (I) 不同的次數進行多項式擬合,「外插」(Extrapolation)時得到的差異相當大 範例7-17:polyfit02.m load census.mat cdate2 = min(cdate):(max(cdate)+30); curve = zeros(7, length(cdate2)); for i = 1:7 curve(i,:) = polyval(polyfit(cdate,pop,i), cdate2); end plot(cdate, pop,'o', cdate2, curve); legend('實際資料', 'order=1', 'order=2','order=3', 'order=4', 'order=5', 'order=6', 'order=7'); xlabel('年份'); ylabel('人口(單位:百萬)');
不同次數的比較 (II) 次數越高,越能逼近給定的資料 不表示預測的準確度會提高 高次多項式對雜訊(Noise)敏感度較高,容易產生不準確的外插預測
多項式擬合的誤差及錯誤 上述範例中,常出現 MATLAB 警告訊息 表示實際資料的不足 或資料數值範圍的過大或過小 在進行多項式擬合時,最好先將資料進行必要的平移、放大或縮小範圍