Learning Polynomials 台大生機系 方煒.

Presentation on theme: "Learning Polynomials 台大生機系 方煒."— Presentation transcript:

1 Learning Polynomials 台大生機系 方煒

2 多項式乘除 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.

3 Multiply Polynomials % a=x^2+0x+1 and b=x^2+0x-1 a=[1,0,1];
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)

4 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

5 多項式乘除: 範例 >>a = [9,-5,3,7]; >>b = [6,-1,2];
>>product = conv(a,b) product = >>[quotient, remainder] = deconv(a,b) quotient = remainder =

6 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

7 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 =

8 Find the polynomial from the 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)

9 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)’)

10 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’)

11 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)

12 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-’)

13 First Derivative of a polynomial array
ap=polyder(a)

14 練習 求下列方程式的根，並用poly函數來驗證 X3+13X2+52X+6=0
驗證 (20X3-7X2+5X+10)(4X2+12X-3)= 80X5+212X4-124X3+121X2+105X-30 驗證 (12X3+5X2-2X+3)/(3X2-7X+4) = 4X+11 餘 (59X-41)

15 The polyfit function Command p = polyfit(x,y,n) Description
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.

16 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.

17 時間，秒 溫度，oF 145 620 130 2266 103 3482 90 詳見 coffee.m 範例：咖啡冷卻問題 室溫 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)’),…

18 範例：咖啡冷卻問題 (續) 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)'););

19

20 An experiment to verify Torricelli’s principle.
f = c * sqrt(p) 流體的體積流量率 f 與流經管道的壓力差 p 的平方根成正比 f = r * sqrt(V) 流體的體積流量率 f 與水槽中液體體積 V 的平方根成正比

21 An experiment to verify Torricelli’s principle.
15 6 12 7 9 8 (1). 右表數據是否符合托里切利原理？如果不符合，請求出下列公式之 r 與 n： f = r * Vn (2). 廠商欲製造可容納 36 杯容量的咖啡壺，需要多長可充滿一杯？ 詳見 coffee2.m

22 Flow rate and fill time for a coffee pot.

23 最小的殘差平方和 used to fit a function f (x). It minimizes the sum of the squares of the residuals, J. J is defined as 有最小J值的函數有最好的配適結果. m S J = [ f (xi ) – yi ]2 i=1

24 Illustration of the least squares criterion.

25 The least squares fit for the example data
y 2 5 6 10 11 m S J = [ f (xi ) – yi ]2 i=1 J = (0m+b-2)2+(5m+b-6)2+(10m+b-11)2

26 The polyfit function is based on the least-squares method
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.

27 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([ ]); subplot(2,2,2); plot(xp,yp(2,:),x,y,’o’),axis([ ]); subplot(2,2,3); plot(xp,yp(3,:),x,y,’o’),axis([ ]); subplot(2,2,4); plot(xp,yp(4,:),x,y,’o’),axis([ ]); disp(J)

28 Regression using polynomials of first through fourth degree.
J = 72 J = 57 J = 42 J = 4.7

29 多項式迴歸的次數越高，當然越準，但是也越不能用來作外插 六個數據點，取五次式，可得 perfect fit.

30 評估配適曲線的品質 S S y值到平均值（y）差異量的平方和為 S, 公式如下： m S = (yt – y )2 i=1
how much the data is spread around the mean S S = (yt – y )2 i=1 殘差的平方和為 J, 公式如下： how much of the data spread is unaccounted for by the model m S J = [f(xt) – yi]2 i=1 判定係數 Coefficient of determination，又稱 r平方值 J / S indicates the fractional variation unaccounted for by the model J S r 2 = 1 -

31 判定係數 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.

32 迴歸與數據正確度 使用高次的多項式時，要留意求出的係數不可隨意減少位數 次數用到三次，多數時候仍可保留其強健性 範例：數據如下表

33 使用 format long 假設各項係數小數點下只取八位數，其餘四捨五入

34 Effect of coefficient accuracy on a sixth-degree polynomial
上圖:小數點下取14位數. 下圖:小數點下只取8位數的係數.

35 為避免高次項的誤差，寧可使用兩個三次式來替代六次式
0 < x < 15 15 < x < 40

36 HW #3-1 請在上圖繪出前述兩個三次式，另外也繪出當各項係數只取小數點下前八位時的另兩條曲線。 計算原來的J, S 與 r2

37 Avoiding high degree polynomials: Use of two cubics to fit data.

38 HW#3-2 細菌繁殖模型 某種細菌的濃度對時間的繁殖數據如右。
請使用一次、二次與三次式及指數配適(fitting) ，並檢查殘差以決定何者最適合。其中，指數配適使用 y = 10(mt+a)的公式。 0min 6 ppm 10min 350ppm 1 13 11 440 2 23 12 557 3 33 685 4 54 14 815 5 83 15 990 6 118 16 1170 7 156 17 1350 8 210 18 1575 9 282 19 1830

39 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);

40 Using Residuals: Residual plots of four models.

41 HW #3-3 某醫學工程儀器可以很快量出電壓反應的響應曲線數據如下表：
當 t -> inf 時電壓會達穩態常數值，T 是電壓達穩態值 95%所需花費的時間，假設 T為 3 sec. 已知一階與二階指數公式如下： v(t)=a1+a2e(-3t/T) v(t)=a1+a2e(-3t/T)+a3te(-3t/T) t(s) 0.3 0.8 1.1 1.6 2.6 3 v(V) 0.6 1.28 1.5 1.7 1.75 1.8

42 一階模式： Xa=y’ 得到 a1=2.0258, a2=-1.9307 t=[0,0.3,0.8,1.1,1.6,2.3,3];
X=[ones(size(t));exp(-t)]; a=X\y’ 得到 a1=2.0258, a2=

43 二階模式： 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= , a3=0.8885

44 Result of HW #3 look like this. Please write the script.

45 Basic Fitting Interface
使用三次樣條 cubic spline 法或最多十次的多項式法來作數值配適 對同一組數據允許同時繪出多組配適圖形 繪出殘差圖形 允許檢查數值配適的結果 允許對數據作內差與外差 可繪出數值配適後圖形與殘差的範數 可將配適與計算結果存入工作區 (workspace)

46 The Basic Fitting interface
按下右箭頭幾次看看各頁面

47 A figure produced by the Basic Fitting interface