Transcript Learning Polynomials 台大生機系 方煒 MATLAB 工程應用 – Learning Polynomials

MATLAB 工程應用 – Learning Polynomials
Learning
Polynomials
台大生機系 方煒
MATLAB 工程應用 – Learning Polynomials
多項式乘除
Polynomial Multiplication and Division
The function conv(a,b) computes the product of the
two polynomials described by the coefficient arrays a
and b. The two polynomials need not be the same
degree. The result is the coefficient array of the
product polynomial.
The function [q,r] = deconv(num,den) computes
the result of dividing a numerator polynomial, whose
coefficient array is num, by a denominator polynomial
represented by the coefficient array den. The quotient
polynomial is given by the coefficient array q, and the
remainder polynomial is given by the coefficient array
r.
MATLAB 工程應用 – Learning Polynomials
Multiply Polynomials

The command conv multiplies two polynomial coefficient
arrays and returns the coefficient array of their product.
% a=x^2+0x+1 and b=x^2+0x-1
 a=[1,0,1];
 b=[1,0,-1];
 c=conv(a,b)

MATLAB 工程應用 – Learning Polynomials
Divide Polynomials



% a=x^2+x+1 and b=x+1
a=[1,1,1];b=[1,1];
[q,r]=deconv(a,b)
q=[1,0] and r=[0,0,1],
That is q =x + 0 = x and r = 0x2 + 0x + 1 = 1
x  x 1
1
 x
x 1
x 1
2
MATLAB 工程應用 – Learning Polynomials
多項式乘除: 範例
>>a = [9,-5,3,7];
>>b = [6,-1,2];
>>product = conv(a,b)
product =
54
-39
41
29
-1
14
>>[quotient, remainder] = deconv(a,b)
quotient =
1.5
-0.5833
remainder =
0
0
-0.5833
8.1667
MATLAB 工程應用 – Learning Polynomials
Polynomial Roots
The function roots(a)computes the roots of a polynomial
specified by the coefficient array a. The result is a
column vector that contains the polynomial’s roots.
For example,
>>r = roots([2, 14, 20])
r =
-2
-7
MATLAB 工程應用 – Learning Polynomials
Polynomial Coefficients
The function poly(r)computes the coefficients of the
polynomial whose roots are specified by the vector r.
The result is a row vector that contains the polynomial’s
coefficients arranged in descending order of power.
For example,
>>c = poly([-2, -7])
c =
1
7
10
MATLAB 工程應用 – Learning Polynomials
Roots

Roots of a Polynomial
% find the roots of a polynomial
 p=[1,2,-13,-14,24];
 r=roots(p)


Find the polynomial from the roots
r=[1,2,3];
 p=poly(r)

MATLAB 工程應用 – Learning Polynomials
Plotting Polynomials
The function polyval(a,x)evaluates a polynomial at
specified values of its independent variable x, which can
be a matrix or a vector. The polynomial’s coefficients of
descending powers are stored in the array a. The result
is the same size as x.
To plot the polynomial f (x) = 9x3 – 5x2 + 3x + 7 for
2 ≤ x ≤ 5, you type
>>a = [9,-5,3,7];
>>x = [-2:0.01:5];
>>f = polyval(a,x);
>>plot(x,f),xlabel(’x’),ylabel(’f(x)’)
–
MATLAB 工程應用 – Learning Polynomials
Plotting Polynomials with the
polyval Function
To plot the polynomial 3x5 + 2x4 – 100x3 + 2x2 – 7x + 90
over the range –6 x  6 with a spacing of 0.01, you type
>>x = [-6:0.01:6];
>>p = [3,2,-100,2,-7,90];
>>plot(x,polyval(p,x)),…
xlabel(’x’),ylabel(’p’)
MATLAB 工程應用 – Learning Polynomials
Evaluate a Polynomial

If you have an array of x-values and you want to
evaluate a polynomial at each one








% define the polynomial
a=[1,2,-13,-14,24];
% load the x-values
x=-5:.01:5;
% evaluate the polynomial
y=polyval(a,x);
% plot it
plot(x,y)
MATLAB 工程應用 – Learning Polynomials
Fitting Data to a Polynomial









x=linspace(0,pi,50);
% make a sine function with 1% random error on it
f=sin(x)+.01*rand(1,length(x));
% fit to the data
p=polyfit(x,f,4);
% evaluate the fit
g=polyval(p,x);
% plot fit and data together
plot(x,f,’r*’,x,g,’b-’)
MATLAB 工程應用 – Learning Polynomials
First Derivative of a polynomial array
a=[1,1,1,1]
 ap=polyder(a)

MATLAB 工程應用 – Learning Polynomials
練習
1.
求下列方程式的根，並用poly函數來驗證
X3+13X2+52X+6=0
2.
驗證 (20X3-7X2+5X+10)(4X2+12X-3)=
80X5+212X4-124X3+121X2+105X-30
3.
驗證 (12X3+5X2-2X+3)/(3X2-7X+4) =
4X+11 餘 (59X-41)
MATLAB 工程應用 – Learning Polynomials
The polyfit function
Command
Description
p =
polyfit(x,y,n)
Fits a polynomial of degree n to
data described by the vectors x and
y, where x is the independent
variable. Returns a row vector p of
length n + 1 that contains the
polynomial coefficients in order of
descending powers.
MATLAB 工程應用 – Learning Polynomials
The polyfit Function
Syntax: p = polyfit(x,y,n)
where x and y contain the data, n is the order of the
polynomial to be fitted, and p is the vector of
polynomial coefficients.
The linear function: y = mx + b. In this case the
variables w and z in the polynomial w = p1z+ p2 are
the original data variables x and y, and we can find
the linear function that fits the data by typing p =
polyfit(x,y,1). The first element p1 of the vector
p will be m, and the second element p2 will be b.
MATLAB 工程應用 – Learning Polynomials
詳見 coffee.m
範例：咖啡冷卻問題

時間，秒
溫度，oF
0
145
620
130
2266
103
3482
90
室溫 68oF下瓷杯內的咖啡由 145oF逐漸冷卻，溫度對應於
經過時間的數據如上，請建立溫度對時間的函數，並與
此模式預測，經過多少時間，杯內溫度會達120oF?
time=[0,620,2266,3482]; temp=[145,130,103,90]; temp=temp-68;
subplot(2,2,1);
plot(time,temp, time,temp,’o’),xlabel(‘Time (sec)’),…
ylabel(‘Relative Temperature (deg F)’)
subplot(2,2,2);
semilogy(time,temp, time,temp,’o’),xlabel(‘Time (sec)’),…
ylabel(‘Relative Temperature (deg F)’)
MATLAB 工程應用 – Learning Polynomials
範例：咖啡冷卻問題 (續)
P=polyfit(time, log10(temp),1);
m=p(1); b=10^p(2);
t_120=(log10(120-68)-log10(b))/m
%show thederived curve and estimated point
t=[0:10:4000];
T=68+b*10.^(m*t);
subplot(2,2,3);
semilogy(t,T-68,time,temp,'o',t_120,120-68,'+'),xlabel('Time (sec)'),...
ylabel('Relative Temperature (deg F)');
%show in rectilinear scales
subplot(2,2,4);
plot(t,T,time,temp+68,'o',t_120,120,'+'),xlabel('Time (sec)'),...
ylabel('Temperature (deg F)'););
MATLAB 工程應用 – Learning Polynomials
MATLAB 工程應用 – Learning Polynomials
An experiment to verify Torricelli’s principle.
f = c * sqrt(p)
流體的體積流量率
f 與流經管道的壓
力差 p 的平方根成
正比
f = r * sqrt(V)
流體的體積流量率
f 與水槽中液體體
積 V 的平方根成
正比
MATLAB 工程應用 – Learning Polynomials
An experiment to verify Torricelli’s principle.
(1). 右表數據是否符合托里
切利原理？如果不符合，請
求出下列公式之 r 與 n：
水槽中液體體積
(V)，杯
充滿一杯所需時
間 (t)，秒
15
6
12
7
f = r * Vn
9
8
6
9
(2). 廠商欲製造可容納 36
杯容量的咖啡壺，需要多長
可充滿一杯？
詳見 coffee2.m
MATLAB 工程應用 – Learning Polynomials
Flow rate and fill time for a coffee pot.
MATLAB 工程應用 – Learning Polynomials
最小的殘差平方和
used to fit a function f (x). It minimizes the sum of the
squares of the residuals, J. J is defined as
m
J 
S
i1
[ f (xi ) – yi ]2
有最小J值的函數有最好的配適結果.
MATLAB 工程應用 – Learning Polynomials
Illustration of the least squares criterion.
MATLAB 工程應用 – Learning Polynomials
The least squares fit for the example data
x
y
0
2
5
6
10
11
m
J 
S
i1
[ f (xi ) – yi ]2
J = (0m+b-2)2+(5m+b-6)2+(10m+b-11)2
J
 250m  30b  280  0
m
J
 30m  6b  38  0
b
MATLAB 工程應用 – Learning Polynomials
The polyfit function is based on the least-squares
method. Its syntax is
p =
polyfit(x,y,n)
Fits a polynomial of degree n to
data described by the vectors x
and y, where x is the independent
variable. Returns a row vector p of
length n+1 that contains the
polynomial coefficients in order of
descending powers.
MATLAB 工程應用 – Learning Polynomials
Script file
x=[1:9];
y=[5,6,10,20,28,33,34,36,42];
xp=[1:0.01:9];
for k=1:4
coeff=polyfit(x,y,k);
yp(k,:)=polyval(coeff,xp);
J(k)=sum((polyval(coeff,x)-y).^2);
end
subplot(2,2,1); plot(xp,yp(1,:),x,y,’o’),axis([0 10 0 50]);
subplot(2,2,2); plot(xp,yp(2,:),x,y,’o’),axis([0 10 0 50]);
subplot(2,2,3); plot(xp,yp(3,:),x,y,’o’),axis([0 10 0 50]);
subplot(2,2,4); plot(xp,yp(4,:),x,y,’o’),axis([0 10 0 50]);
disp(J)
MATLAB 工程應用 – Learning Polynomials
Regression using polynomials of first through fourth degree.
J = 72
J = 57
J = 42
J = 4.7
MATLAB 工程應用 – Learning Polynomials
多項式迴歸的次數越高，當然越準，但是也越不能用來作外插
六個數據點，取五次式，可得 perfect fit.
MATLAB 工程應用 – Learning Polynomials
評估配適曲線的品質
y值到平均值（y）差異量的平方和為 S, 公式如下：
m
S
S
(yt – y
)2
i1
殘差的平方和為 J, 公式如下：
m
J
S
[f(xt) – yi]2
i1
how much the data is
spread around the
mean
how much of the
data spread is
unaccounted for by
the model
判定係數 Coefficient of determination，又稱 r平方值
r2 1 J
S
J / S indicates the
fractional variation
unaccounted for by the
model
MATLAB 工程應用 – Learning Polynomials
判定係數
For a perfect fit, J 0 and thus r 2 1. Thus the closer
r 2 is to 1, the better the fit. The largest r 2 can be is 1.
It is possible for J to be larger than S, and thus it is
possible for r 2 to be negative. Such cases, however,
are indicative of a very poor model that should not be
used.
As a rule of thumb, a good fit accounts for at least 99
percent of the data variation. This value corresponds
to r 2  0.99.
MATLAB 工程應用 – Learning Polynomials
迴歸與數據正確度



使用高次的多項式時，要留意求出的係數不可隨
意減少位數
次數用到三次，多數時候仍可保留其強健性
範例：數據如下表
MATLAB 工程應用 – Learning Polynomials
使用 format long

假設各項係數小數點下只取八位數，其餘四捨五入
MATLAB 工程應用 – Learning Polynomials
Effect of coefficient accuracy on a sixth-degree polynomial
上圖:小數點下取14位數. 下圖:小數點下只取8位數的係數.
MATLAB 工程應用 – Learning Polynomials
為避免高次項的誤差，寧可使用兩
個三次式來替代六次式

0 < x < 15

15 < x < 40
MATLAB 工程應用 – Learning Polynomials
HW #3-1



請在上圖繪出前述兩個三次式，另外也繪
出當各項係數只取小數點下前八位時的另
兩條曲線。
計算原來的J, S 與 r2
計算後來的J, S 與 r2
MATLAB 工程應用 – Learning Polynomials
Avoiding high degree polynomials:
Use of two cubics to fit data.
MATLAB 工程應用 – Learning Polynomials
HW#3-2 細菌繁殖模型


某種細菌的濃度對時
間的繁殖數據如右。
請使用一次、二次與
三次式及指數配適
(fitting) ，並檢查殘差
以決定何者最適合。
其中，指數配適使用 y
= 10(mt+a)的公式。
0min
1
2
3
4
5
6
7
8
9
6 ppm
13
23
33
54
83
118
156
210
282
10min
11
12
13
14
15
16
17
18
19
350ppm
440
557
685
815
990
1170
1350
1575
1830
MATLAB 工程應用 – Learning Polynomials
Script file 未完成
x=[0:19];
y=[6,13,23,33,54,83,118,156,210,282,…
350,440,557,685,815,990,1170,1350,1575,1830];
%curve fitting
p1=polyfit(x,y,1);
p2=polyfit(x,y,2);
p3=polyfit(x,y,3);
p4=polyfit(x,log10(y),1);
%residuals
res1=polyval(p1,x)-y;
res2=polyval(p2,x)-y;
res3=polyval(p3,x)-y;
res4=polyval(p4,x)-y;
subplot(2,2,1);plot(x,res1); subplot(2,2,2);plot(x,res2);
subplot(2,2,3);plot(x,res3); subplot(2,2,4);plot(x,res4);
MATLAB 工程應用 – Learning Polynomials
Using Residuals: Residual plots of four models.
MATLAB 工程應用 – Learning Polynomials
HW #3-3

某醫學工程儀器可以很快量出電壓反應的
響應曲線數據如下表：
t(s)
v(V)

0
0
0.3
0.6
0.8
1.28
1.1
1.5
1.6
1.7
2.6
1.75
3
1.8
當 t -> inf 時電壓會達穩態常數值，T 是電壓達
穩態值 95%所需花費的時間，假設 T為 3 sec.
已知一階與二階指數公式如下：
v(t)=a1+a2e(-3t/T)
v(t)=a1+a2e(-3t/T)+a3te(-3t/T)
MATLAB 工程應用 – Learning Polynomials

一階模式：
Xa=y’
t=[0,0.3,0.8,1.1,1.6,2.3,3];
y=[0,0.6,1.28,1.5,1.7,1.75,1.8];
X=[ones(size(t));exp(-t)];
a=X\y’

得到 a1=2.0258, a2=-1.9307
MATLAB 工程應用 – Learning Polynomials

二階模式：
Xa=y’
t=[0,0.3,0.8,1.1,1.6,2.3,3];
y=[0,0.6,1.28,1.5,1.7,1.75,1.8];
X=[ones(size(t)); exp(-t); t.exp(-t)];
a=X\y’

得到 a1=1.7496, a2=-1.7682, a3=0.8885
MATLAB 工程應用 – Learning Polynomials
Result of HW #3 look like this. Please write the script.
MATLAB 工程應用 – Learning Polynomials
Basic Fitting Interface







使用三次樣條 cubic spline 法或最多十次的
多項式法來作數值配適
對同一組數據允許同時繪出多組配適圖形
繪出殘差圖形
允許檢查數值配適的結果
允許對數據作內差與外差
可繪出數值配適後圖形與殘差的範數
可將配適與計算結果存入工作區 (workspace)
MATLAB 工程應用 – Learning Polynomials
The Basic Fitting interface

在工具列中
選取 Basic
Fitting 選項

按下右箭頭
幾次看看各
頁面
MATLAB 工程應用 – Learning Polynomials
A figure produced by the Basic Fitting interface