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.

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)

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

多項式乘除: 範例 >>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

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

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

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)

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

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

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)

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

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

練習 求下列方程式的根，並用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)

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.

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.

時間，秒 溫度，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)’),…

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

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

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

Flow rate and fill time for a coffee pot.

最小的殘差平方和 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

Illustration of the least squares criterion.

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

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.

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)

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

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

評估配適曲線的品質 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 -

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

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

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

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

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

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

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

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

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

Using Residuals: Residual plots of four models.

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

一階模式： 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=-1.9307

二階模式： 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

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

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

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

A figure produced by the Basic Fitting interface