程式與科學計算 Lecture 2

Size: px

Start display at page:

Download "程式與科學計算 Lecture 2"

豚灾缪
6 years ago
Views:

1 多項式迴歸解非線性方程式數值積分微分解 ODE Matlab Math tools

2 多項式回歸 (regression):polyfit 找出最能夠代表 x, y 兩變數關係的 n 次多項式假設 x 向量存有 x 變數的觀測資料,y 向量有 y 變數的觀測資料, 想要找出 x 與 y 之間的關係, 統計上會使用回歸計算 (regression) p = polyfit(x, y, n) n=1, 猜測 x 與 y 應該是線性相關, 即 y=ax+b, 結果回傳 p 向量會有兩個元素,p(1) 為係數 a,p(2) 為係數 b n=2, 猜測 x 與 y 應該是二次項關係, 即 y=ax 2 +bx+c, 結果回傳 p 向量會有三個元素,p(1)=a,p(2)=b,p(3)=c 依此類推,p 向量是用 n 次多項式 ax n +bx n-1 +cx n-2 + 回歸所得的係數, 長度為 n+1 求出回歸的 n 次多項式之後, 可以用 polyval(p,x), 以儲存在 p 向量的係數, 算出此 n 次多項式在 x 向量各元素的值

3 y 多項式回歸 : 範例為了舉例方便, 先製造一組 x, y 資料, 各有 100 個資料點,x 為 0~10 之間的 100 個隨機值,y=cos(x) 但加上隨機的誤差值 : x=rand(1,100)*10; %data points of x (random between 0 and 10) y=cos(x)+2*rand(1,100); % data points of y (cos(x) with random noise) plot(x,y,'*') x

4 y 多項式回歸 : 範例接下來用線性 ( 一次 ) 二次三次多項式進行資料回歸, 並且把回歸結果與原始資料點畫在同張圖上 p1=polyfit(x,y,1); p2=polyfit(x,y,2); p3=polyfit(x,y,3); data n=1 fit n=2 fit n=3 fit xfit=0:0.1:10; yfit1=polyval(p1,xfit); yfit2=polyval(p2,xfit); yfit3=polyval(p3,xfit); plot(x,y,'*',xfit,yfit1,'r',xfit,yfit2,'g',xfit,yfit3,'k') legend('data','n=1 fit','n=2 fit', 'n=3 fit') x

5 解非線性方程式 :fsolve 找出非線性方程 F(x)=0 的解 x = fsolve(@fun, 求解的非線性方程 (F(x)) x0:initial guess of solution x:f(x)=0 的解另外可以用 optimset 設定顯示疊代求解的過程 opt=optimset( Display, iter ) x = fsolve(@fun, x0, opt) ( 關於 optimset 的眾多設定, 可參考指令說明 )

6 解非線性方程組 : 範例一求解猜測解為 x=1 x x 1 = 0 ( 直接傳入匿名函數 )( 只適用 Matlab 2011b 版本 ) >> y >> x = fsolve(y, 1) x =

7 求解猜測解為 x 1 = -5, x 2 = -5 解非線性方程組 : 範例二 2x 1 x 2 e x 1 = 0 x 1 + 2x 2 e x 2 = 0 ( 利用自訂函數 ) 開一個 m 檔, 用 function 設定求解的函數, 一個方程式為 F 的其中一個 column function F = myfun(x) F = [2*x(1) - x(2) - exp(-x(1)); -x(1) + 2*x(2) - exp(-x(2))]; 在主程式 ( 另一個 m 檔 ) 或指令視窗, 將自訂函數送入 fsolve 指令求解, 給定猜測值 -5, -5, 並設定顯示疊代過程 x0 = [-5; -5]; % starting guess opt=optimset('display','iter'); x = fsolve(@myfun,x0,opt)

8 求解結果 After several iterations, fsolve finds an answer: Norm of First-order Trust-region Iteration Func-count f(x) step optimality radius e e e e e e Equation solved. fsolve completed because the vector of function values is near zero as measured by the default value of the function tolerance, and the problem appears regular as measured by the gradient. x =

9 Matlab 提供幾種常見的數值積分方法 : 梯形法 trapz(x) 用一次多項式直線逼近數值積分適應式辛普森法 adaptive Simpson quad(x) 用二次多項式曲線逼近 global adaptive quadrature integral(x) 自動判斷高階多項式逼近

10 數值積分 (1): 梯形法 trapz Z = trapz(x, Y) X 向量 : 在 X 座標上的分隔位置, Y 向量 : 每個 X 位置的 f(x) 值 Z 變數 : 以梯形法計算 f(x) 在積分區間 (X(1)=a~X(n)=b) 的定積分結果範例一 (x,y 為資料點 ) >> x=[ ]; >> y=[ ]; >> area1=trapz(x,y) area1 = 範例二 (y 為 x 的函數 ) >> x=linspace(0,pi,25); >> y=sin(x); >> area2=trapz(x,y) area2 =

11 數值積分 (2): 辛普森法 quad Z = quad(fun, a, b) fun 是要積分的函數,a,b 是積分區間,Z 是積分結果範例二 ( 利用自訂函數 ) 範例一 ( 直接傳入匿名函數 ) ( 只適用 Matlab 2011b 版本 ) >> y=@(x) 1/(x^3-2x-5) >> Q = quad(y,0,2) Q = 開一個 m 檔, 用 function 設定要積分的函數 function y = myfun(x) y = 1./(x.^3-2*x-5); 在主程式 ( 另一個 m 檔 ) 或指令視窗, 將自訂函數送入 quad 指令 >> Q = quad(@myfun,0,2) Q =

12 數值積分 (3):integral Z = integral(fun, a, b) fun 是要積分的函數,a,b 是積分區間,Z 是積分結果可計算暇積分積分下限接近奇異點 singularity 等狀況範例 : e x2 dx 0 >> y exp(-x.^2); >> Q = integral(y,0,inf) Q =

13 數值微分 :diff, gradient dfdx = diff(f)./dx: forward difference 計算一次微分 diff(f) 計算 F 向量前後兩個元素的差值回傳結果 dfdx 向量會比 F 向量少一個元素 dx 是微分間距 ( 可能是純量或向量 ) dfdx = diff(f,n)./(dx.^n) forward difference 計算 n 次微分 diff(f,2) = diff(diff(f)) df(x) dx = F x i+1 F(x i ) x d n F(x) dx n dfdx = gradient(f, dx): central difference 一次微分 F 是一維向量, 回傳結果 dfdx 的長度與 F 相同 [dfdx, dfdy] = gradient(f, dx, dy) 計算曲面的梯度向量 F 是裝有曲面資料的二維矩陣 dfdx 是 X 方向的微分,dFdY 是 Y 方向的微分 df(x) dx = F x i+1 F(x i 1 ) 2 x F(x, y) = F F i + x y j

14 數值微分範例 : y=sin(x) 的一次二次微分 ni=100; x=linspace(-pi,pi,ni); dx=(2*pi)/((ni-1)); y=sin(x); dy=gradient(y,dx); ddy=gradient(dy,dx); plot(x,y) hold on plot(x,dy,'r-') plot(x,ddy,'k-') legend('y=sin(x)','y''','y"') hold off

15 數值解常微分方程組 (ODE):ode45 t: 產生時間向量 [t, [t0,t1], X0) X:solution 矩陣, 每一個 column 求解的 ODE 方程組 [t0 t1]: 積分時間區間 X0: 每個應變數初始值 X(t) = (x 1 t, x 2 t ) dx 1 (t) = f dt 1 (x 1 (t), x 2 (t), t) dx 2 (t) = f dt 2 (x 1 (t), x 2 (t), t) I. C. t = t0, X t0 = (x 1 t0, x 2 t0 )

16 解 ODE 方程組 : 範例求解下列常微分方程式, 畫出 t=0 ~10 時 x(t) 與 y(t) 的解 : F t = (x t, y t ) I. C. : t = 0, x 0 = 0, y 0 = 1 dx(t) = xcos t + y sin t dt dy(t) = xsin t + y cos t dt 步驟一 : 編輯 ODE 方程組的自訂函數 ( 假設取名 myode.m): function dfdt=myode(t,f) dfdt(1) = F(1)*cos(t)+F(2)*sin(t); dfdt(2) = -F(1)*sin(t)+F(2)*cos(t); dfdt = dfdt ; % change to column vector

17 y(t) x, y 解 ODE 方程組 : 範例步驟二, 求解, 畫出結果 : 回到指令視窗或主程式 ( 另一個.m 檔 ), 用 ode45 指令求出 ODE 方程組在時間 0~10 之間的解, 初始值設定 x=0, y=1 畫出 x, y 的時間變化序列, 以及 x-y 平面上的 phase diagram [t,f]=ode45(@myode, [0 10],[0 1]); x=f(:,1); y=f(:,2); plot(t,x,'r',t,y, g'); legend( x, y ); plot(x,y) y1 x y2 y t x(t)

18 Matlab 工具庫 (Toolbox) Matlab 有許多外掛或第三方開發的工具庫可以使用, 很多都需要額外購買安裝開啟 help window, 左邊的列表會顯示目前有安裝的工具庫, 點選之後會有使用說明函數功能介紹實際使用範例學校 / 系上的 Matlab 可能會安裝的工具庫有 : Bioinformatics Communications Control System Curve Fitting Filter Design Fixed-Point Image Processing Model Predictive Control Neural Network Optimization Parallel Computing Partial Differential Equation Robust Control Signal Processing Statistics Symbolic Math System Identification Wavelet

<4D6963726F736F667420576F7264202D2032303037C4EAC6D5CDA8B8DFB5C8D1A7D0A3D5D0C9FAC8ABB9FACDB3D2BBBFBCCAD4CEC4BFC6D7DBBACDCAD4BEEDBCB0B4F0B0B82DD6D8C7ECBEED2E646F63>

<4D6963726F736F667420576F7264202D2032303037C4EAC6D5CDA8B8DFB5C8D1A7D0A3D5D0C9FAC8ABB9FACDB3D2BBBFBCCAD4CEC4BFC6D7DBBACDCAD4BEEDBCB0B4F0B0B82DD6D8C7ECBEED2E646F63> 2007 年普通高等学校招生全国统一考试 ( 重庆卷 ) 文综试卷第一部分本部分共 35 题, 每题 4 分, 共 140 分在每题给出的四个选项中, 只有一项最符合题目的要求的读图 1, 回答 1-3 题 1. 某两洲面积之和与某大洋面积十分接近, 它们是