Introduction of MATLAB© MATLAB 2014 workshop Introduction of MATLAB© Speaker: Hsing-Yu Wang (王馨右) E-mail f00524011@ntu.edu.w Date: 09/12/2014 Perhaps we should begin. Good afternoon, ladies and gentlemen. Thanks for coming. My name is Shih-Wei Liu. This afternoon I would like to present A Dynamic Model with Structured Recurrent Neural Network to Predict Blood Glucose in a patient with Type 1 Diabetes Mellitus. PSE Laboratory Department of Chemical Engineering National Taiwan University

Introduction 應用領域 高科技運算程式解決方案 Technical Computing 控制設計 Control Design 訊號處理及通訊 Signal Processing and Communications 影像處理 Image Processing 測試及量測 Test Measurement 生物運算 Computational Biology 財務模型及分析 Financial Modeling and Analysis 教育專區 Education

Fundamentals of MATLAB Outline 1 Fundamentals of MATLAB 2 Plotting Here is my outline. This presentation is about 15 mins and will be divided into four parts. First of all, I will start off by filling you on the background to this field. Then I will elaborate on the methodology of this research. Next, I would like to report on the results of a case study. Finally, I will make a brief conclusion of my presentation today. 3 Implemetation 4 Summary

軟體取得 1.計中網頁(http://www.cc.ntu.edu.tw/chinese/) 2.使用方式:第一次使用首先點選APPShare，非第一次使用請直接跳至第13頁

軟體取得

軟體取得

軟體取得

軟體取得

軟體取得

軟體取得

軟體取得

軟體取得 10.非第一次使用 待下載啟動程式碼完成(右下角顯示100%),即會開 啟應用程式. 安裝軟體時,需使用擁有管理者(Administrator)權 限的帳戶. 注意事項: 1.此為“串流”(stream)程式碼到PC執行的方式,於 第一次使用新功能時,需下載程式碼, 請耐心等候. 2.若您的個人電腦已安裝相同應用程式, 建議先解除 該軟體的安裝. 3.若有安裝錯誤訊息出現,請跳過該訊息.

MATLAB 簡介 什麼是 MATLAB MATrix LABoratory 以矩陣為運算基礎的數學軟體 直譯式的數學程式語言 提供完備的數學函數，且能讓使用者定義自己的函數。 提供矩陣運算、數值分析、 統計分析、訊號處理、圖形繪製(2D/3D)……等功能

MATLAB Interface Workspace Current Folder Command Window Command History Details

Scalar Arithmetic Operations Symbol Operation MATLAB Form ^ exponentiation: ab a^b * multiplication: ab a*b / right division: a/b = a/b \ left division: a\b = a\b + addition: a + b a+b - subtraction: a – b a-b Example: >> (5*2+3.5)/5 ans = 2.7000

Order of Precedence Example: >> 5^2*2 ans = 50 >> 5*2^2 20 Operation First Parentheses, evaluated starting with the innermost pair Second Exponentiation, evaluated from left to right Third Multiplication and division with equal precedence (left to right) Fourth Addition and subtraction with equal precedence (left to right) Example: >> 5^2*2 ans = 50 >> 5*2^2 20

矩陣的基本資料結構 直接列出矩陣內的元素產生矩陣： Example: 用逗號 (,) 或空白區隔同一列的元素 用分號 (;) 或 Enter 鍵代表一列的結束 矩陣的開始和結束，使用中括號 ([]) 表示 Example: >> A=[1 2 3] A = 1 2 3 >> B=[1 2 3; 4 5 6] B = 4 5 6

矩陣的基本資料結構 (cont’d) 直接列出矩陣內的元素產生矩陣： Example: 或使用下列描述法構成： X=起始值：增加值：結束值 1 2 3 4 5 6 7 8 9 10 >> A=[1:2:10] 1 3 5 7 9

矩陣的基本資料結構 (cont’d) 利用 MATLAB 的敘述或函式產生： Example: >> A=ones(3,4) 1 1 1 1 >> B=linspace(1,10,4) B = 1 4 7 10

矩陣的索引或下標 Ref: http://www.cs.nthu.edu.tw/~jang第九章:矩陣的處理與運算

Predefined Constants Command Description ans Temporary variable containing the most recent answer eps Specifies the accuracy of floating point precision i,j The imaginary unit Inf Infinity, ∞ NaN Indicates an undefined numerical result (Not a Number) pi The number π

變數命名規則與使用 Example: 第一個字母必需是英文字母。 字母間不可留空格。 字母大小寫有異 最多只能有 31 個字母，MATLAB 會忽略多餘字母（在 MATLAB 第 4 版，則是 19 個字母）。　 MATLAB 在使用變數時，不需預先經過變數宣告（Variable Declaration）的程序，而且所有數值變數均以預設的 double 資料型式儲存。 Example: PSE123 （O） 3PSE21 （X） PSE_123（O） PSE-123 （X）

Matrix operations Symbol Operation Form Example + Scalar-array addition A+b [6,3]+ 2 =[8, 5] - Scalar-array subtraction A-b [6,3]- 2 =[9,11] Array addition A+B [6,3]+[3,8] Array subtraction A-B [6,3]-[3,8] =[3,-5] * Matrix multiplication A*B [1,2]*[3,4]' =[11] \ Matrix left division x=A-1B A\B [1,2]\[10] =[0 5]' .* Array multiplication A.*B [3,2].*[2,4] =[6, 8] ./ Array right division A./B [4,8]./[2,4] =[2, 2] .\ Array left division A.\B [2,4].\[4,8] .^ Array exponent A.^B [3,5].^2 =[9,25] 2.^[3,5] =[8,32] [3,2].^[2,3] =[9, 8]

Relational operators Example: Relational Operator Meaning < Less than <= Less than or equal to > Greater than >= Greater than or equal to == Equal to ~= Not equal to Example: >> x = [ 6, 3, 9]; >> y = [14, 2, 9]; >> z = (x < y) z = 1 0 0

矩陣的邏輯運算 Example: 矩陣的邏輯運算： &：And 運算(=and) |：Or 運算(=or) ~：Not 運算(=not) any：任何一個元素不為零 all：所有元素都不為零 Example: Xor Logical exclusive-OR >> X=[0 1 2 3] X = 0 1 2 3 >> any(X) ans = 1 >> any(X==4) ans = 0

Some Useful Array Functions Command Description find(x) Computes an array containing indices of nonzero elements of array x length(A) Computes either the number of elements of A if A is a vector or the largest value of m or n if A is an m × n matrix size(A) Returns a row vector [m n] containing the sizes of the m × n array max(A) Returns a row vector containing largest element in each column of A min(A) Same as max(A) but returns minimum values Inv(A) Return the inverse of the square matrix A det(A) Return the determinant of the square matrix A eig(A) Returns a vector of the eigenvalues of matrix A sort(A) Sorts each column of array A in ascending order and returns an array the same size as A sum(A) Sums the elements in each column of array A and returns a row vector containing the sums

常見特殊矩陣 zeros(m,n) %維度為mxn, 元素全為0的矩陣 ones(m,n) %維度為mxn, 元素全為1的矩陣 eye(n) %維度為nxn的identity矩陣 rand(m,n) %維度為mxn,於[0,1]均勻分佈的亂數矩陣 randn(m,n) %維度為mxn,於[0,1]高斯分佈的亂數矩陣 [] 空矩陣 diag 對角矩陣 linspace 產生線性的空間向量

Special Characters =：指定（Assign）運算 (ex: B=A) ' ： 矩陣轉置 (ex: B=A') ..：代表上一層目錄 (ex: cd ..) ...：繼續，代表指令未完，接到下一列 ,：元素分隔符號 ;：代表命令結束，且不印出執行結果 %：註解列

Pre-allocation 在 MATLAB 運算過程，可以隨時改變每一個矩陣的維度，以容納新的資料。但是每一次增大矩陣的容量時，MATLAB 就必須要向電腦的作業系統索取記憶體，因此會造成程式執行效率的降低，而且造成記憶體的分散使用現象（Memory Segmentation），此種由於記憶體的動態配置而造成的分散使用現象，可能導致實際雖仍有充足的記憶體，但卻未有連續空間以處理較大矩陣的現象。 基於以上種種原因，若預先知道使用矩陣的維度大小，則最好實施矩陣的預先配置（Pre-allocation），可用的指令有 zeros、ones、cell（用於配置異值陣列）及 struct（用於配置結構陣列）等。例如，欲計算矩陣 A 的 n 次方的行列式，其中 n = 1~100，可輸入如下： Example: A = [1 2 3; 4 5 6; 7 8 9]; det_value = zeros(1, 100); for i = 1:100 det_value(i) = det(A^i); end

工作空間與變數的儲存及載入 Workspace 中變數的保存方式： save：將 workspace 中的變數存到磁碟上 load：由磁碟上的資料檔載入變數到 workspace 中 不加任何選項（Options）時，save 指令會將工作空間內的變數以二進制（Binary）的方式儲存至副檔名為 mat 的檔案

工作空間與變數的儲存及載入 (cont’d) 使用 save 保存變數： save：將所有變數以二進位格式存至 “matlab.mat”中 save fname：將所有變數存至”fname.mat”中 save fname X Y：將變數 X 與 Y 存至”fname.mat”中 save fname X -ascii：使用 8 位數文字格式將變數X 存至”fname”中 save fname X -ascii -double：使用 16 位數文字格式將變數 X 存至”fname”中 save name.txt X -ascii -double –tabs ：使用 16 位數文字格式將變數 X 存至”fname”.txt的文字檔中

工作空間與變數的儲存及載入 (cont’d) 使用 load 載入變數： load：由『matlab.mat』中載入所有變數 load fname：由『fname.mat』中載入所有變數 load fname X Y：由『fname.mat』中載入變數 X 與 Y load fname -ascii：由文字檔『fname』中載入變數 load fname.txt：由文字檔『fname.txt』中載入變數

Commands for Managing the Work Session Description clc Clears the command window clear Removes all variables from memory clear var1 var2 Removes var1 and var2 from memory exist(’name’) Determine if a file or variable exists having the name ’name’ quit Stop MATLAB : Colon; generates an array having regularly spaced elements , Comma; separates elements of an array ; Semicolon; suppresses screen printing; also denotes a new row in an array . . . Ellipsis; continues a line to delay execution

MATLAB 下的程式寫作 MATLAB 程式的基本：M-file M-file 的特性： M-file 的種類： Script （命令稿/命令巨集） Function （函式）

開啟M-File

MATLAB 下的程式寫作 (cont’d) Script： 用來連續執行一連串的指令 所有執行中使用或產生的變數皆存在於 workspace，因此變數容易互相覆蓋而造成程式錯誤 沒有輸入與輸出變數 Function： 用來定義一個自訂函式 除了全域變數(global variables)，所有執行中使用或產生的變數皆不存在於 workspace，因此不會和 MATLAB 基本工作空間的變數相互覆蓋 可以有輸出與輸入變數

MATLAB 下的程式寫作 Example: 第一列為函數定義列（Function Definition Line） FindRoots.m function [x1 x2] = FindRoots(a,b,c) x1=(-b+(b^2-4*a*c)^0.5)/(2*a); x2=(-b-(b^2-4*a*c)^0.5)/(2*a); 第一列為函數定義列（Function Definition Line） 定義函數名稱（FindRoots，最好和檔案的檔名相同） 輸入引數（ [x1 x2] ） 輸出引數（a,b,c） function為關鍵字 第二列開始為函數主體（Function Body） 描述函數運算過程，並指定輸出引數的值

呼叫Function [Output1, Output2, …] = CommandName(input1, input2, …) CommandName: built-in or user-defined function Example: FindRoots.m function [x1 x2] = FindRoots(a,b,c) x1=(-b+(b^2-4*a*c)^0.5)/(2*a); x2=(-b-(b^2-4*a*c)^0.5)/(2*a); >> a=1; >> b=3; >> c=1; >> [x1 x2] =FindRoots(a,b,c) x1 = -0.3820 x2 = -2.6180

次函數 次函數(sub-function):一個M 檔案可以有一個以上的次函數 Example: TEST.m >> v=1; function [x1 x2] = TEST(H) [a b c]= Cal_abc(H) [x1 x2] = FindRoots(a,b,c) >> v=1; >> [x1 x2] = TEST(v) x1 = 1.2573 x2 = -1.5907 Cal_abc.m function [a b c] = Cal_abc(H) a=H^3; b=H/3; c=H-3;

全域變數 (global variable) 此變數可在不同的函數及底稿(m-file)使用這在工程計算相當有用 Example: TEST.m global a b c a=18.3036 b=3816.44 c=-46.13 for i=1:11 t=373.2+10*(i-1) %t in Kelvin psat(i)=antoine(t) % psat in mmHg end Antoine.m function p_vap=antoine(t) global a b c p_vap=exp(a-b/(t+c))

MATLAB 下的流程控制 程式的流程控制（Flow Control）： 迴圈指令（Loop Command）： for …… end while …… end 條件指令（Branching Command）： if … elseif … else … end switch … case … otherwise … end

MATLAB 下的流程控制 - for for 迴圈： 語法 1： 語法 2： for 變數=向量 運算式 end for 變數=矩陣 Example: >>x=1; >> for i=1:5 >> x=x*i; >> end >> disp(x) 120 >> x=[1 3; 2 4]; >> for i=x >> fprintf('%d\t',i) >> end 1 2 3 4

MATLAB 下的流程控制 - while while 迴圈： 語法： while 條件式 運算式 end % 1+ ε < 1? Example: % 1+ ε < 1? e= 1; while (1+e) > 1 %判斷式成立 e = e/2; end e = e*2

MATLAB 下的流程控制 - if if 條件指令： 語法： if 條件式一 運算式一 elseif 條件式二 運算式二 運算式三 else 運算式四 end Example: if x > 10 y = log(x) elseif x >= 0 y = sqrt(x) else y = exp(x) - 1 end

MATLAB 下的流程控制 – switch switch 條件指令： 語法： switch 運算式 case 值一 運算式一 Example: switch 條件指令： 語法： switch 運算式 case 值一 運算式一 case 值二 運算式二 otherwise 運算式三 end for month=1:12 switch month case {3,4,5} season='Spring'; case {6,7,8} season='Summer'; case {9,10,11} season='Autumn'; otherwise season='Winter'; end fprintf('Month %d is %s.\n',month,season) Result: Month 1 is Winter. Month 2 is Winter. Month 3 is Spring. Month 4 is Spring. …..

MATLAB 下的流程控制 break: Terminate execution of for or while loop return: Return to invoking function Example: for i = 1:1000 if prod(1:i) > 1e10 break;% %跳出迴圈 end disp(i) 14 det.m function d = det(A) %DET det(A) is the determinant of A. if isempty(A) d = 1; return else ... end

Exercise clear all A=[1 2 3 4]'; B(1)=A(1);C(1)=A(1); for i=2:length(A) if A(i)>=3 B(1,i)=A(i)+C(i-1); else B(1,i)=A(i)+B(i-1); end C(1,i)=myfun(B(i)); Ans=B*C' myfun.m function C=myfun(B) while B > 2 B=B-1; end C=B^2;

Plotting

Some MATLAB Plotting Commands Description axis([xmin xmax ymin ymax]) Sets the minimum and maximum limits of x- and y-axes. fplot(’string’, [xmin xmax]) Performs intelligent plotting of functions, where ’string’ is a text string that describes the function to be plotted and [xmin xmax] specifies the minimum and maximum values of the independent variable. The range of the dependent variable can also be specified. In this case the syntax is fplot (’string’, [xmin xmax ymin ymax]). grid Displays gridlines at the tick marks corresponding to the tick labels. Type grid on to add gridlines; type grid off to stop plotting gridlines. When used by itself, grid switched this feature on or off. plot(x,y) Generates a plot of the array x on rectilinear axes. legend('leg1','leg2', ...) Creates a legend using the strings leg1, leg2, and so on and specifies its placement with the mouse. print Prints the plot in the figure window. title('text') Puts text in a title at the top of a plot. xlabel('text') Adds a text label to the x-axis(the abscissa). ylabel('text') Adds a text label to the y-axis(the ordinate). hold Freezes the current plot for subsequent graphics commands. subplot(m,n,p) Splits the figure window into an array of subwindows with m rows and n columns and directs the subsequent plotting commands to the pth subwindow. text(x,y,'text') Places the string text in the figure window at point specified by coordinates x, y.

二維基本繪圖簡介 plot -- 繪圖指令的根本： 語法：plot(X,Y, format)：將向量 Y 對向量 X 做圖 Example: >> X=[0:0.1:2*pi]; >> Y=sin(X); >> plot(X, Y, 'bx-') ‘color marker line-type')

Data Markers, Line Types, and Colors Dot (．) . Solid line - Black k Asterisk (*) * Dashed line -- Blue b Cross (×) x Dash-dotted line -. Cyan c Circle (。) o Dotted line : Green g Plus sign (+) + Magenta m Square (□) s Red r Diamond (◇) d White w Five-point star (★) p Yellow y Tips: Command Line:help plot

圖形標題與座標軸 X=[0:0.1:2*pi]; Y=sin(X); plot(X, Y, 'b--','linewidth',3) Example: X=[0:0.1:2*pi]; Y=sin(X); plot(X, Y, 'b--','linewidth',3) title('Y vs X','fontsize',15) xlabel('X','fontsize',15) ylabel('f(X)','fontsize',15) legend('sin(X)','fontsize',15) set(gca,'fontsize',15)

單一圖片呈現多個圖形 使用 plot 的第三種形式，一次畫出多個數對。 使用 hold 將圖形「固定」住，重疊多個數對的繪圖結果。 使用 subplot 分割繪圖視窗，同時呈現多個子圖形。

單一圖片呈現多個圖形 - plot plot -- 繪圖指令的根本： 語法：plot(X1,Y1, format1, X2,Y2, format2...)：將 Yi 對 Xi 以 formati 做圖 Example: >> X=[0:0.1:2*pi]; >> Y=sin(X); >> Z=cos(X); >> plot(X, Y, 'bx-', X, Z, 'r:')

單一圖片呈現多個圖形 - hold 使用 hold 重疊多個數對圖形： Example: >> X=[0:0.1:2*pi]; >> Y=sin(X); >> Z=cos(X); >> plot(X, Y, 'b--') >> hold on; >> xlabel('X') >> ylabel('f(X)') >> plot(X, Z, 'r:') >> text(0.2,0,'sin(X)') >> text(0.5*pi+0.2,0,'cos(X)')

單一圖片呈現多個圖形 - subplot 使用 subplot 分割視窗呈現數個子圖形： >> X=[0:0.1:2*pi]; Example: >> X=[0:0.1:2*pi]; >> Y=sin(X); >> Z=cos(X); >> subplot(211); >> plot(X, Y, 'b--') >> xlabel('X') >> ylabel('sin(X)') >> subplot(2,1,2); >> plot(X, Z, 'r:') >> ylabel('cos(X)')

Plot tools Show plot tool Hide plot tool

Enhance your figure Generate Code

三維基本繪圖簡介 – line plot t = 0:pi/50:10*pi; plot3(sin(t),cos(t),t) grid on axis square

三維基本繪圖簡介 – mesh plot % meshgrid:Generate X and Y arrays for 3-D plots [X,Y] = meshgrid(-3:.125:3); Z = peaks(X,Y); % mesh, meshc, meshz Mesh plots mesh(X,Y,Z); axis([-3 3 -3 3 -10 5])

Symbolic calculation

subs函數 subs(E,old,new)可以將算式 E 中 old 的部分以 new取代，其中 old 這個部分可以是符號變數或者是符號表示式，並且 new 這個部分可以是符號變數、符號表式、矩陣、數值的值，或者是數值的矩陣。 >> syms x y >> E = x^2+6*x+7; >> F = subs(E,x,y) F = y^2+6*y+7 >> F = subs(E,x,5) 62 Example:

solve 函數 有三種使用solve函數的方式。例如，要求解方程式 x + 5 = 0 >> eq1 = 'x+5=0'; >> solve(eq1) ans = -5 >> solve('x+5=0') ans = -5 >> syms x >> solve(x+5) ans = -5

limit函數 >> syms a x >> limit(sin(a*x)/x) ans = a 函數 limit(E) 可以用來求出當 x → 0的時候 E 的極限。 >> syms a x >> limit(sin(a*x)/x) ans = a 另外limit(E,v,a)這個形式可以用來求出當u → a 時的極限。舉例來說， >> syms h x >> limit((x-3)/(x^2-9),3) 1/6 >> limit((sin(x+h)-sin(x))/h,h,0) cos(x)

Symbolic Calculation (diff, int) syms x f=sin(x); c=diff(f) d=diff(f,2) c_value=subs(c,2) d_value=subs(d,2) c = cos(x) d = -sin(x) c_value = -0.4161 d_value = -0.9093 syms x y f=sin(x)+cos(y); c=diff(f,x) c2=diff(f,x,2) d=diff(f,y) d2=diff(f,y,2) c = cos(x) c2 = -sin(x) d = -sin(y) d2 = -cos(y) syms x f=-2*x/(1+x^2)^2; c=int(f) d=int(f,0,1) c = 1/(1+x^2) d = -1/2

拉氏轉換 Example: >> syms b t >> laplace(t^3) ans = 6/s^4 >> laplace(exp(-b*t)) 1/(s+b) >> laplace(sin(b*t)) b/(s^2+b^2)

反拉氏轉換 Example: >> syms b s >> ilaplace(1/s^4) ans = t^3/6 >> ilaplace(1/(s+b)) exp(-b*t) >> ilaplace(b/(s^2+b^2)) sin(b*t)

特徵方程式以及特徵根 Example: >> syms k >> A = [0 ,1;-k, -2]; >> poly(A) ans = x^2+2*x+k >> solve(ans) [ -1+(1-k)^(1/2) ] [ -1-(1-k)^(1/2) ]

Curve fitting (regression)

Functions for Polynomial Regression 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. [p,s, mu] = 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 and a structure s for use with polyval to obtain error estimates for predictions. The optional output variable mu is two-element vector containing the mean and standard deviation of x. Example: The linear function y = mx + b gives a straight line when plotted on rectilinear axes. Its slope is m and its intercept is b. → p = polyfit(x, y, 1)

Curve Fitting Tools Choose model Create data set

Curve Fitting Tools (cont’d)

Solving ODE by m-file

用 MATLAB 解 ODE 解 ode 可用的 function：ode23(), ode113(), ode15s(), ode23s(), ode23t(), ode23tb(), ode45() 基本語法（以 ode45() 為例）： [T,Y] = ode45('F',TSPAN,Y0,OPTIONS) or [T,Y] = ode45(@F,TSPAN,Y0,OPTIONS) T：時間 F：ode function M-file 定義檔檔名 TSPAN：積分時間 Y0：初始值 OPTIONS：特殊選項 P1, P2 …：其餘欲傳給 ode function M-file 的參數

Example 1 rigid.m function dy = rigid(t,y) dy = zeros(3,1); % a column vector dy(1) = y(2) * y(3); dy(2) = -y(1) * y(3); dy(3) = -0.51 * y(1) * y(2);

Example 1 (cont’d) >> [T,Y] = ode45(@rigid,[0 12],[0 1 1]); >> plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.')

Example 2 先將二階化成標準格式，令 vdp1.m function dy = vdp1(t, y) mu=1 dy = [ y(2) mu*(1-y(1)^2)*y(2)-y(1)]; Command Lines: ode45('vdp1', [0 25], [3 3]); or ode45(@vdp1, [0 25], [3 3]);

Example 3 – VDV reaction vdv.ode.m function xdot = vdv_ode(t,x) ca=x(1); cb=x(2); k1=5/6; k2=5/3; k3=1/6; fov=4/7; caf=10; dcadt=fov*(caf-ca)-k1*ca-k3*(ca^2); dcbdt=-fov*cb+k1*ca-k2*cb; xdot=[dcadt;dcbdt];

Example 3 (cont’d) x0=[2;1.117]; tspan=[0 5]; [t,x]=ode45(@vdv_ode,tspan,x0); subplot(2,1,1) plot(t,x(:,1)) xlabel('t') ylabel('ca') subplot(2,1,2) plot(t,x(:,2)) ylabel('cb')

Summary of coding Writing down before coding Looking for help Using recognizable name of variable Checking matrix dimension Checking loop and iteration

Backup materials

Numerical differentiation True Slope

Numerical Differentiation(diff) d = diff(x) = [x(2) − x(1), x(3) − x(2), . . . , x(n) − x(n − 1)] Example: x=0:0.1:1; y=[0.5 0.6 0.7 0.9 1.2 1.4 1.7 2.0 2.4 2.9 3.5]; dx=diff(x); dy=diff(y); dydx=diff(y)./diff(x) dydx = 1.0000 1.0000 2.0000 3.0000 2.0000 3.0000 3.0000 4.0000 5.0000 6.0000

Numerical Differentiation Functions Command Description b = polyder(p) Returns a vector b containing the coefficients of the derivative of the polynomial represented by the vector p. b=polyder(p2,p1) Returns a vector b containing the coefficients of the derivative of the polynomial that is the product of the polynomials represented by p2 and p1. [n,d]=polyder(p2,p1) Returns the vectors n and d containing the coefficients of the numerator and denominator polynomials of the derivative of the quotient p2/p1, where p2 and p1 are polynomials.

Numerical Differentiation Polynomial Derivatives p1 = [5,2]; p2 = [10,4,-3]; der2 = polyder(p2) prod = polyder(p2,p1) [num,den] = polyder(p2,p1) der2 = 20 4 prod = 150 80 -7 num = 50 40 23 den = 25 20 4

Numerical integration

Numerical Integration Command Description trapz(x,y) Uses trapezoidal integration to compute the integral of y with respect to x, where the array y contains the function values at the points contained in the array x quad(’function’,a,b,tol) Uses an adaptive Simpson’s rule to compute the integral of the function ’function’ with a as the lower integration limit and b as the upper limit. The parameter tol is optional. tol indicates the specified error tolerance. quad1(’function’,a,b,tol) Uses Lobatto quadrature to compute the integral of the function ’function’. The rest of the syntax is identical to quad.

Numerical Integration x = linspace(0, pi, 10); ans = y = sin(x); 1.9797 trapz(x, y) x = linspace(0, pi, 100); ans = y = sin(x); 1.9998 x = linspace(0, pi, 1000); ans = y = sin(x); 2.0000

Interpolation

Polynomial Interpolation Functions(1D) y_est=interp1(x,y,x_est,'method') This command specifies alternate methods. The default is linear interpolation. Available methods are: nearest - nearest neighbor interpolation linear - linear interpolation spline - piecewise cubic spline interpolation (SPLINE) pchip - piecewise cubic Hermite interpolation (PCHIP) cubic - same as ’pchip’

Polynomial Interpolation Functions(2D) z_est=interp2(x,y,z,x_est,y_est,’method’) This command specifies alternate methods. The default is linear interpolation. Available methods are: nearest - nearest neighbor interpolation linear - bilinear interpolation (or bilinear) spline - cubic spline interpolation (SPLINE) cubic - bicubic interpolation

Minimization (optimization)

Minimizing A Function of Variables (fminsearch) minimum of , near (0,0) function f = f146(x) f = x(1).*exp(-x(1).^2-x(2).^2); %% end, stored in f145.m fminsearch('f146',[0,0]) % [f, fval] = fminsearch(’f146’,[0,0]) ans = -0.7071 0.0000

Other Built-in Functions fmincon: find minimum of constrained nonlinear multivariable function fminbnd: find minimum of single-variable function on fixed interva Isqlin: constrained linear least-squares problems Equation lsqnonlin: Solve nonlinear least-squares (nonlinear data-fitting) problems quadprog: Solve quadratic programming problems

Optimization Tool

Solving algebraic equation

Finding Roots of Equations Equation Type solver Polynomial roots roots continuous function of one variable fzero nonlinear equations fsolve (fminsearch…)

roots roots: polynomial roots syntax： x = roots([c1 c2 c3…. cn+1]) Example: p = [1 -6 -72 -27] r = roots(p) r = 12.1229 -5.7345 -0.3884

Some Polynomial Functions Command Description poly(r) Computes the coefficients of the polynomial whose roots are specified by the array r. The resulting coefficients are arranged in descending order of powers. polyval(a,x) Evaluates a polynomial at specified values of its independent variable x. The polynomial’s coefficients of descending powers are stored in the array a. The result is the same size as x. roots(a) Computes an array containing the roots of a polynomial specified by the coefficient array a.

Polynomial Roots x3 - 7x2 + 40x -34 =0 a = [1, -7, 40, -34],... roots(a) a = 1 -7 40 -34 ans = 3.000 + 5.000i 3.000 - 5.000i 1.000 roots = 1, 3 ± 5i r = [1, 3+5i, 3-5i],... poly(r) r = 1.00 3.00+5.00i 3.00-5.00i ans = 1 -7 40 -34

fzero 函數 fzero: find root of continuous function of one variable syntax： x = fzero(@FUN, x0) x：解 FUN：定義問題的 function M-file 檔名 x0：初始猜值矩陣

Example - fzero >> fzero('func01',[3 6]) ans = 3.9270 func01.m function f=func01(x) f = sin(x) - cos(x); >> fzero('func01',[3 6]) ans = 3.9270 >> fzero('func01',[0 3]) ans = 0.7854

fsolve 函數 fsolve: 解聯立方程式 語法： x = fsolve(@FUN, x0) x：解 FUN：定義問題的 function M-file 檔名 x0：初始猜值矩陣 Reamark：fsolve() 是 optimization toolbox 的指令

Example - fsolve A=[1 7 5; 2 3 4;9 6 8]; b=[ 2 ;4;6]; c=inv(A)*b Method 1 Method 2 function F=func01(X) A=[1 7 5; 2 3 4;9 6 8]; b=[ 2 ;4;6]; F=A*X-b A=[1 7 5; 2 3 4;9 6 8]; b=[ 2 ;4;6]; c=inv(A)*b x= -0.4000 -1.1077 2.0308 x0=[1 1 1] x=fsolve('func01',x0) x= -0.4000 -1.1077 2.0308 % Use help fsolve for details

Example Van de Vusse reaction in an isothermal CSTR:

Example vdv_sol_ad.m x0=[2;1.117]; k1=5/6; k2=5/3; k3=1/6; fov=4/7; caf=10; options=optimset('Display','iter'); [x]=fsolve(@vdv_ae,x0,options,k1,k2,k3,fov,caf); cas=x(1) cbs=x(2) vdv_ad.m function F = vdv_ae(x,k1,k2,k3,fov,caf) F = [fov*(caf-x(1))-k1*x(1)-k3*(x(1)^2); -fov*x(2)+k1*x(1)-k2*x(2)];

Simulink

Simulink 簡介 什麼是 Simulink ? MATLAB 提供的一個 toolbox 應用於動態系統模擬 提供 GUI，使動態系統的建立如同畫方塊圖（Block Diagram）

啟動 Simulink 啟動 Simulink： 直接在 command window 鍵入 「simulink」

啟動 Simulink Simulink 的外觀：

動態模型的建立 動態模型的建立步驟： 開啟一個空白的模型（Model） 開啟適當的 Library 視窗，找出所有需要的 Block，將該 Block 拖入或複製至空白模型視窗 用滑鼠雙擊各 Block 開啟對話視窗，適當設定各 Block 所需的參數 利用滑鼠連接各個 Block，建立完整的動態系統

動態模型的建立

Solving ODE using dee

dee 啟動指令 A window named dee will appear Command Line: dee

模型建立

執行指令與作圖 k1=5/6; k2=5/3; k3=1/6; fov=4/7; caf=10; x0=[2;1.117]; sim('vdv_sim_4',[0 5]) subplot(2,1,1) plot(t,ca) xlabel('t') ylabel('ca') subplot(2,1,2) plot(t,cb) ylabel('cb')

Debug

Debug mode breakpoint Run vdv_sol_ae.m

Debug mode (cont’d) Step Continue

Solving ODE using S-function

Example – VDV reaction Enter search term Constant S-function

Command Line: open sfuntmpl.m sfuntmpl.m will appear

vdv_sfun.m function [sys,x0,str,ts] = vdv_sfun(t,x,u,flag,x0) switch flag, case 0, [sys,x0,str,ts]=mdlInitializeSizes(x0); case 1, sys=mdlDerivatives(t,x,u); case 2, sys=mdlUpdate(t,x,u); case 3, sys=mdlOutputs(t,x,u); case 4, sys=mdlGetTimeOfNextVarHit(t,x,u); case 9, sys=mdlTerminate(t,x,u); otherwise DAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag)); end

vdv_sfun.m (cont’d) function [sys,x0,str,ts]=mdlInitializeSizes(x0) sizes = simsizes; sizes.NumContStates = 2; sizes.NumDiscStates = 0; sizes.NumOutputs = 2; sizes.NumInputs = 2; sizes.DirFeedthrough = 0; sizes.NumSampleTimes = 1; % at least one sample time is needed sys = simsizes(sizes); str = []; ts = [0 0];

vdv_sfun.m (cont’d) function sys=mdlDerivatives(t,x,u) ca=x(1); cb=x(2); fov=u(1); caf=u(2); k1=5/6; k2=5/3; k3=1/6; dcadt=fov*(caf-ca)-k1*ca-k3*(ca^2); dcbdt=-fov*cb+k1*ca-k2*cb; sys = [dcadt;dcbdt]; function sys=mdlUpdate(t,x,u) sys = [];

vdv_sfun.m (cont’d) function sys=mdlOutputs(t,x,u) sys = [x(1);x(2)]; function sys=mdlGetTimeOfNextVarHit(t,x,u) sys = []; function sys=mdlTerminate(t,x,u)

Example – VDV reaction (cont’d) 1 2 4 3

Example – VDV reaction (cont’d) Command Line: x0=[2;1.117]

Example – VDV reaction (cont’d) sfun_plot.m subplot(2,1,1) plot(t,ca) xlabel('t') ylabel('ca') subplot(2,1,2) plot(t,cb) ylabel('cb')

MV step change Step

Command Line: x0=[3;1.117] Start Simulation Reach New Steady State

Refine Simulink file

Refine Simulink file (cont’d) Create Subsystem in mouse right click menu

Refine Simulink file (cont’d) Mask Subsystem in mouse right click menu

Refine Simulink file (cont’d) Add

Refine Simulink file (cont’d) Right double click on subsystem will appear the following window:

Simulating LTI model

LTI model 的建立 [ys,ts]=step(vdv_tf); plot(ts,ys) [yi,ti]=impulse(vdv_tf); plot(ti,yi) pole(vdv_tf) tzero(vdv_tf) A=[-2.4048 0;0.83333 -2.2381]; B=[7;-1.117]; C=[0 1]; D=[0]; vdv_ss=ss(A,B,C,D) num=[-1.117 3.1472]; den=[1 4.6429 5.3821]; vdv_tf=tf(num,den) vdv_zpk=zpk(2.817,[-2.238 -2.405],-1.117) vdv_tf1=tf(vdv_ss) vdv_zpk1=zpk(vdv_tf1) vdv_tfd=c2d(vdv_tf,1,'zoh') vdv_tf2=tf(num,den,'Inputdelay',1)

LTI model 進階分析 kc=1; figure(1) bode(kc*vdv_tf) figure(2) [mag,phase,w]=bode(kc*vdv_tf) subplot(2,1,1),loglog(w,squeeze(mag)) subplot(2,1,2),semilogx(w,squeeze(phase)) [Gm,Pm,Wg,Wp]=margin(mag,phase,w) nyquist(kc*vdv_tf)

Simulate LTI model using Simulink

LTI model 的建立 Transfer Fcn

Closed-loop simulation

PID Controller

Look Under Mask

Solving PDE

Numerical Methods for PDE(pdepe) pdepe solves PDEs of the form with boundary conditions

Numerical Methods for PDE(pdepe) function [c,b,s] = eqn1(x,t,u,DuDx) %EQN1: MATLAB function M-le that %species %a PDE in time and one space %dimension. c = 1; b = DuDx; s = 0; function [pl,ql,pr,qr] = bc1(xl,ul,xr,ur,t) %BC1: MATLAB function M-le that %species boundary conditions %for a PDE in time and one space %dimension. pl = ul; ql = 0; pr = ur-1; qr = 0; function value = initial1(x) %INITIAL1: MATLAB function M-le %that species the initial condition %for a PDE in time and one space %dimension. value = 2*x/(1+x^2);

Numerical Methods for PDE(pdepe) %PDE1: MATLAB script M-le that solves and plots %solutions to the PDE stored in eqn1.m m = 0; x = linspace(0,1,20); t = linspace(0,2,10); %Solve the PDE u = pdepe(m,@eqn1,@initial1,@bc1,x,t); %Plot solution surf(x,t,u); title('Surface plot of solution.'); xlabel('Distance x'); ylabel('Time t');