Transcript 點此下載

MATLAB 2012 workshop
Introduction of
Speaker: Hsing-Yu Wang (王馨右)
E-mail
[email protected]
Date:
09/06/2013
PSE Laboratory
Department of Chemical Engineering
National Taiwan University
©
MATLAB
Introduction
應用領域
高科技運算程式解決方案 Technical Computing
控制設計 Control Design
訊號處理及通訊 Signal Processing and
Communications
影像處理 Image Processing
測試及量測 Test Measurement
生物運算 Computational Biology
財務模型及分析 Financial Modeling and Analysis
教育專區 Education
2 /80
Outline
1
Fundamentals of MATLAB
2
Plotting
3
Implemetation
4
Summary
3 /80
軟體取得
1.計中網頁(http://www.cc.ntu.edu.tw/chinese/)
2.使用方式:第一次使用首先點選APPShare，非第一次使用
請直接跳至第13頁
4 /80
軟體取得
5 /80
軟體取得
6 /80
軟體取得
7 /80
軟體取得
8 /80
軟體取得
9 /80
軟體取得
10 /80
軟體取得
11 /80
軟體取得
10.非第一次使用
待下載啟動程式碼完成(右下角顯示100%),即會開
啟應用程式.
安裝軟體時,需使用擁有管理者(Administrator)權
限的帳戶.
注意事項:
1.此為“串流”(stream)程式碼到PC執行的方式,於
第一次使用新功能時,需下載程式碼, 請耐心等候.
2.若您的個人電腦已安裝相同應用程式, 建議先解除
該軟體的安裝.
3.若有安裝錯誤訊息出現,請跳過該訊息.
12 /80
MATLAB 簡介
• 什麼是 MATLAB
– MATrix LABoratory
– 以矩陣為運算基礎的數學軟體
– 直譯式的數學程式語言
– 提供完備的數學函數，且能讓使用者定義自己
的函數。
– 提供矩陣運算、數值分析、 統計分析、訊號處
理、圖形繪製(2D/3D)……等功能
13 /80
MATLAB Interface
Current
Folder
Details
Workspace
Command
Window
Command
History
14 /80
Scalar Arithmetic Operations
Symbol
Operation
MATLAB Form
^
exponentiation: ab
a^b
*
multiplication: ab
a*b
/
right division: a/b =
a
b
b
left division: a\b =
a
a/b
+
addition: a + b
a+b
-
subtraction: a – b
a-b
\
a\b
• Example:
>> (5*2+3.5)/5
ans =
2.7000
15 /80
Order of Precedence
Precedence
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
ans =
20
16 /80
矩陣的基本資料結構
• 直接列出矩陣內的元素產生矩陣：
– 用逗號 (,) 或空白區隔同一列的元素
– 用分號 (;) 或 Enter 鍵代表一列的結束
– 矩陣的開始和結束，使用中括號 ([]) 表示
• Example:
>> A=[1 2 3]
A=
1 2 3
>> B=[1 2 3; 4 5 6]
B=
1 2 3
4 5 6
17 /80
矩陣的基本資料結構 (cont’d)
• 直接列出矩陣內的元素產生矩陣：
– 或使用下列描述法構成：
X=起始值：增加值：結束值
• Example:
>> A=[1:10]
A=
1 2 3
4
5
>> A=[1:2:10]
A=
1 3 5
7
9
6
7
8
9
10
18 /80
矩陣的基本資料結構 (cont’d)
• 利用 MATLAB 的敘述或函式產生：
• Example:
>> A=ones(3,4)
A=
1 1 1 1
1 1 1 1
1 1 1 1
>> B=linspace(1,10,4)
B=
1 4 7 10
19 /80
矩陣的索引或下標
Ref: http://www.cs.nthu.edu.tw/~jang第九章:矩陣的處理與運算
20 /80
Predefined Constants
Command
ans
Description
Temporary variable containing the most recent
answer
eps
i,j
Specifies the accuracy of floating point precision
The imaginary unit 1
Inf
NaN
Infinity, ∞
Indicates an undefined numerical result (Not a
Number)
pi
The number π
21 /80
變數命名規則與使用
•
•
•
•
第一個字母必需是英文字母。
字母間不可留空格。
字母大小寫有異
最多只能有 31 個字母，MATLAB 會忽略多餘字母（在
MATLAB 第 4 版，則是 19 個字母）。
• MATLAB 在使用變數時，不需預先經過變數宣告（
Variable Declaration）的程序，而且所有數值變數均以預
設的 double 資料型式儲存。
• Example:
PSE123 （O）
3PSE21 （X）
PSE_123（O）
PSE-123 （X）
22 /80
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]
=[9,11]
-
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]
=[2, 2]
.^
Array exponent
A.^B
[3,5].^2
=[9,25]
2.^[3,5]
=[8,32]
[3,2].^[2,3]
=[9, 8]
23 /80
Relational operators
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
24 /80
矩陣的邏輯運算
• 矩陣的邏輯運算：
&：And 運算(=and)
|：Or 運算(=or)
~：Not 運算(=not)
any：任何一個元素不為零
all：所有元素都不為零
• Example:
>> X=[0 1 2 3]
X= 0 1 2
>> any(X)
ans = 1
>> any(X==4)
ans = 0
3
25 /80
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
26 /80
常見特殊矩陣
zeros(m,n)
ones(m,n)
eye(n)
rand(m,n)
randn(m,n)
[]
diag
linspace
%維度為mxn, 元素全為0的矩陣
%維度為mxn, 元素全為1的矩陣
%維度為nxn的identity矩陣
%維度為mxn,於[0,1]均勻分佈的亂數矩陣
%維度為mxn,於[0,1]高斯分佈的亂數矩陣
空矩陣
對角矩陣
產生線性的空間向量
27 /80
Special Characters
=：指定（Assign）運算 (ex: B=A)
' ： 矩陣轉置 (ex: B=A')
..：代表上一層目錄 (ex: cd ..)
...：繼續，代表指令未完，接到下一列
,：元素分隔符號
;：代表命令結束，且不印出執行結果
%：註解列
28 /80
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
29 /80
工作空間與變數的儲存及載入
• Workspace 中變數的保存方式：
– save：將 workspace 中的變數存到磁碟上
– load：由磁碟上的資料檔載入變數到
workspace 中
• 不加任何選項（Options）時，save 指令會將工作空間內
的變數以二進制（Binary）的方式儲存至副檔名為 mat 的
檔案
30 /80
工作空間與變數的儲存及載入
(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的文字檔中
31 /80
工作空間與變數的儲存及載入
(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』中載入變數
32 /80
Commands for Managing the
Work Session
Command
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
33 /80
MATLAB 下的程式寫作
• MATLAB 程式的基本：M-file
• M-file 的特性：
– 是純文字文件（只要延伸檔名維持 .m）
– 任何 MATLAB 指令都可以出現在 M-file 中
• M-file 的種類：
– Script （命令稿/命令巨集）
– Function （函式）
34 /80
開啟M-File
35 /80
MATLAB 下的程式寫作 (cont’d)
• Script：
– 用來連續執行一連串的指令
– 所有執行中使用或產生的變數皆存在於 workspace，
因此變數容易互相覆蓋而造成程式錯誤
– 沒有輸入與輸出變數
• Function：
– 用來定義一個自訂函式
– 除了全域變數(global variables)，所有執行中使用或產
生的變數皆不存在於 workspace，因此不會和
MATLAB 基本工作空間的變數相互覆蓋
– 可以有輸出與輸入變數
36 /80
MATLAB 下的程式寫作
• 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);
第一列為函數定義列（Function Definition Line）
•定義函數名稱（FindRoots，最好和檔案的檔名相同）
•輸入引數（ [x1 x2] ）
•輸出引數（a,b,c）
•function為關鍵字
第二列開始為函數主體（Function Body）
•描述函數運算過程，並指定輸出引數的值
37 /80
呼叫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 =
-1.0000 + 1.4142i
x2 =
-1.0000 - 1.4142i
38 /80
次函數
• 次函數(sub-function):一個M 檔案可以有一
個以上的次函數
Example:
TEST.m
function [x1 x2] = TEST(H)
[a b c]= Cal_abc(H)
[x1 x2] = FindRoots(a,b,c)
Cal_abc.m
function [a b c] = Cal_abc(H
a=H^3;
b=H/3;
c=H-3;
>> v=1;
>> [x1 x2] = TEST(v)
x1 =
1.2573
x2 =
-1.5907
39 /80
全域變數 (global variable)
• 全域變數(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))
40 /80
MATLAB 下的流程控制
• 程式的流程控制（Flow Control）：
– 迴圈指令（Loop Command）：
• for …… end
• while …… end
– 條件指令（Branching Command）：
• if … elseif … else … end
• switch … case … otherwise … end
41 /80
MATLAB 下的流程控制 - for
• for 迴圈：
– 語法 1：
for 變數=向量
運算式
end
– 語法 2：
for 變數=矩陣
運算式
end
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
42 /80
MATLAB 下的流程控制 - while
• while 迴圈：
– 語法：
while 條件式
運算式
end
Example:
% 1+ ε < 1?
e= 1;
while (1+e) > 1 %判斷式成立
e = e/2;
end
e = e*2
43 /80
MATLAB 下的流程控制 - if
• if 條件指令：
– 語法：
if 條件式一
運算式一
elseif 條件式二
運算式二
elseif 條件式三
運算式三
else
運算式四
end
Example:
if x > 10
y = log(x)
elseif x >= 0
y = sqrt(x)
else
y = exp(x) - 1
end
44 /80
MATLAB 下的流程控制 – switch
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)
End
Result:
Month 1 is Winter.
Month 2 is Winter.
Month 3 is Spring.
Month 4 is Spring.
…..
45 /80
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
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
46 /80
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));
end
Ans=B*C'
myfun.m
function C=myfun(B)
while B > 2
B=B-1;
end
C=B^2;
47 /80
Plotting
5
0
-5
-10
2
0
-2
-3
-2
-1
0
1
2
3
48 /80
Some MATLAB Plotting Commands
Command
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.
49 /80
二維基本繪圖簡介
• 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')
50 /80
Data Markers, Line Types, and
Colors
Data markers
Line types
Colors
Dot (．)
.
Solid line
-
Black
k
Asterisk (*)
*
Dashed line
--
Blue
b
Cyan
c
Green
g
Cross (×)
x
Dash-dotted line
Circle (。)
o
Dotted line
Plus sign (+)
+
Square (□)
-.
:
Magenta
m
s
Red
r
Diamond (◇)
d
White
w
Five-point star (★)
p
Yellow
y
Tips:
Command Line:help plot
51 /80
圖形標題與座標軸
Y vs X
1
Example:
0.5
f(X)
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)
sin(X)
0
-0.5
-1
0
2
4
X
6
8
52 /80
單一圖片呈現多個圖形
1. 使用 plot 的第三種形式，一次畫出多個數
對。
2. 使用 hold 將圖形「固定」住，重疊多個數
對的繪圖結果。
3. 使用 subplot 分割繪圖視窗，同時呈現多
個子圖形。
53 /80
單一圖片呈現多個圖形 - 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:')
54 /80
單一圖片呈現多個圖形 - 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)')
55 /80
單一圖片呈現多個圖形 - subplot
• 使用 subplot 分割視窗呈現數個子圖形：
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:')
>> xlabel('X')
>> ylabel('cos(X)')
56 /80
Plot tools
Show plot tool
Hide plot tool
57 /80
Enhance your figure
Generate Code
58 /80
三維基本繪圖簡介 – line plot
t = 0:pi/50:10*pi;
plot3(sin(t),cos(t),t)
grid on
axis square
40
30
20
10
0
1
1
0.5
0.5
0
0
-0.5
-0.5
-1
-1
59 /80
三維基本繪圖簡介 – 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);
5
axis([-3 3 -3 3 -10 5])
0
-5
-10
2
0
-2
-3
-2
-1
0
1
2
3
60 /80
Symbolic calculation
61 /80
subs函數
• subs(E,old,new)可以將算式 E 中 old 的部分以
new取代，其中 old 這個部分可以是符號變數或
者是符號表示式，並且 new 這個部分可以是符號
變數、符號表式、矩陣、數值的值，或者是數值
的矩陣。
Example:
>>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)
F =
62
62 /80
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
63 /80
limit函數
函數 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)
ans =
1/6
>>limit((sin(x+h)-sin(x))/h,h,0)
ans =
cos(x)
64 /80
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)
syms x
f=-2*x/(1+x^2)^2;
c=int(f)
d=int(f,0,1)
c=
c=
1/(1+x^2)
d=
-1/2
cos(x)
c2 =
-sin(x)
d=
-sin(y)
d2 =
-cos(y)
65 /80
拉氏轉換
Example:
>>syms b t
>>laplace(t^3)
ans =
6/s^4
>>laplace(exp(-b*t))
ans =
1/(s+b)
>>laplace(sin(b*t))
ans =
b/(s^2+b^2)
66 /80
反拉氏轉換
Example:
>>syms b s
>>ilaplace(1/s^4)
ans =
1/6*t^3
>>ilaplace(1/(s+b))
ans =
exp(-b*t)
>>ilaplace(b/(s^2+b^2)
ans =
sin(b*t)
67 /80
特徵方程式以及特徵根
Example:
>>syms k
>>A = [0 ,1;-k, -2];
>>poly(A)
ans =
x^2+2*x+k
>>solve(ans)
ans =
[ -1+(1-k)^(1/2) ]
[ -1-(1-k)^(1/2) ]
68 /80
Curve fitting (regression)
Input data
6
5
4
3
2
1
0
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
69 /80
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)
70 /80
Curve Fitting Tools
Choose model
Create data set
71 /80
Curve Fitting Tools (cont’d)
72 /80
Solving ODE by m-file
73 /80
用 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 的參數
74 /80
Example 1
y1  y2 y3
y2   y1 y3
y1 (0)  0
y2 (0)  1
y3  0.51y1 y2 y3 (0)  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);
75 /80
Example 1 (cont’d)
>>[t,Y] = ode45(@rigid,[0 12],[0 1 1]);
>>plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.')
1
0.8
0.6
0.4
Y
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
0
2
4
6
t
8
10
12
76 /80
Example 2
y   (1  y 2 ) y  y  0
先將二階化成標準格式，令
 y1  y2

2

y


(1

y
 2
1 ) y2  y1
vdp1.m
function dy = vdp1(t, y)
mu=1
dy = [
y(2)
mu*(1-y(1)^2)*y(2)-y(1)];
y1  y, y2  y
Command Lines:
ode45('vdp1', [0 25], [3 3]);
or
ode45(@vdp1, [0 25], [3 3]);
77 /80
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];
dC A F
 (C Af  C A )  k1C A  k3C A2
dt
V
dCB
F
  CB  k1C A  k2 CB
dt
V
Parameters:
5
5
1 mol
min 1 ; k2  min 1 ; k3 
6
3
6 L  min
Inputs:
k1 
F 4
 min 1 ; C Af  10mol/L
V 7
78 /80
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))
xlabel('t')
ylabel('cb')
79 /80
Summary of coding
•
•
•
•
•
Writing down before coding
Looking for help
Using recognizable name of variable
Checking matrix dimension
Checking loop and iteration
80 /80
Backup materials
81 /80
Numerical differentiation
True Slope
82 /80
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
83 /80
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.
84 /80
Numerical Differentiation
Polynomial Derivatives
p1  5 x  2
p2  10 x 2  4 x  3
dp2
 20 x  4
dx
p1 p2  50 x 3  40 x 2  7 x  6
dp1 p2
 150 x 2  80 x  7
dx
d ( p2 / p1 )
50 x 2  40 x  23

dx
25 x 2  20 x  4
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
85 /80
Numerical integration
86 /80
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.
87 /80
Numerical Integration


0 sin(x)dx   cos(x) 0  cos(0)  cos( )  2
x = linspace(0, pi, 10);
y = sin(x);
trapz(x, y)
ans =
1.9797
x = linspace(0, pi, 100);
y = sin(x);
trapz(x, y)
ans =
1.9998
x = linspace(0, pi, 1000);
y = sin(x);
trapz(x, y)
ans =
2.0000
88 /80
Interpolation
89 /80
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’
90 /80
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
91 /80
Minimization (optimization)
92 /80
Minimizing A Function of Variables
(fminsearch)
minimum of f  xe
 x2  y 2
, 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
93 /80
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
94 /80
Optimization Tool
95 /80
Solving algebraic equation
96 /80
Finding Roots of Equations
Equation Type
Polynomial roots
solver
roots
continuous function of one variable
nonlinear equations
fzero
fsolve (fminsearch…)
97 /80
roots
• roots: polynomial roots
– syntax：
x = roots([c1 c2 c3…. cn+1])
c1xn  c2 xn1  c3 xn2 ....  cn x  cn1
Example:
x3  6 x 2  72 x  27  0
p = [1 -6 -72 -27]
r = roots(p)
r=
12.1229
-5.7345
-0.3884
98 /80
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.
99 /80
Polynomial Roots
x3 - 7x2 + 40x -34 =0
roots = 1, 3 ± 5i
a = [1, -7, 40, -34],...
roots(a)
r = [1, 3+5i, 3-5i],...
poly(r)
a=
1 -7 40 -34
r=
1.00 3.00+5.00i 3.00-5.00i
ans =
3.000 + 5.000i
3.000 - 5.000i
1.000
ans =
1 -7 40 -34
100 /80
fzero 函數
• fzero: find root of continuous function of
one variable
– syntax：
x = fzero(@FUN, x0)
x：解
FUN：定義問題的 function M-file 檔名
x0：初始猜值矩陣
101 /80
func01.m
function f=func01(x)
f = sin(x) - cos(x);
Example - fzero
1.5
1
sin(x) - cos(x)
0.5
0
>> fzero('func01',[3 6])
ans =
3.9270
-0.5
>> fzero('func01',[0 3])
ans =
-1
0.7854
-1.5
0
1
2
3
4
5
x
6
7
8
9
10
102 /80
fsolve 函數
• fsolve: 解聯立方程式
– 語法：
x = fsolve(@FUN, x0)
x：解
FUN：定義問題的 function M-file 檔名
x0：初始猜值矩陣
– Reamark：fsolve() 是 optimization toolbox 的指令
103 /80
Example:
Example - fsolve
1 7 5   x1   2 
2 3 4  x2  4

   
9 6 8   x3   5 
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
 x1  1 7 5   2 
 x2  2 3 4 4
  
  
 x3  9 6 8   5 
x0=[1 1 1]
x=fsolve('func01',x0)
x= -0.4000
-1.1077
2.0308
% Use help fsolve for details
A=[1 7 5; 2 3 4;9 6 8];
b=[ 2 ;4;6];
c=inv(A)*b
1
x= -0.4000
-1.1077
2.0308
104 /80
Example
Van de Vusse reaction in an isothermal CSTR:
dC A F
 (C Af  C A )  k1C A  k3C A2
dt
V
dCB
F
  CB  k1C A  k2 CB
dt
V
Parameters:
5
5
1 mol
min 1 ; k2  min 1 ; k3 
6
3
6 L  min
Inputs:
k1 
F 4
 min 1 ; C Af  10mol/L
V 7
105 /80
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)];
106 /80
Simulink
107 /80
Simulink 簡介
• 什麼是 Simulink ?
– MATLAB 提供的一個 toolbox
– 應用於動態系統模擬
– 提供 GUI，使動態系統的建立如同畫方塊圖
（Block Diagram）
108 /80
啟動 Simulink
• 啟動 Simulink：
– 直接在 command window 鍵入 「simulink」
– 使用工具列的 Simulink 按鈕：
109 /80
啟動 Simulink
• Simulink 的外觀：
110 /80
動態模型的建立
• 動態模型的建立步驟：
– 開啟一個空白的模型（Model）
– 開啟適當的 Library 視窗，找出所有需要的
Block，將該 Block 拖入或複製至空白模型視窗
– 用滑鼠雙擊各 Block 開啟對話視窗，適當設定
各 Block 所需的參數
– 利用滑鼠連接各個 Block，建立完整的動態系
統
111 /80
動態模型的建立
112 /80
Solving ODE using dee
113 /80
dee 啟動指令
A window named dee will appear
Command Line:
dee
114 /80
模型建立
115 /80
執行指令與作圖
3
ca
k1=5/6;
k2=5/3;
k3=1/6;
fov=4/7;
caf=10;
x0=[2;1.117];
sim('vdv_sim_4',[0 5])
2.5
2
0.5
1
1.5
2
2.5
t
3
3.5
4
4.5
5
0
0.5
1
1.5
2
2.5
t
3
3.5
4
4.5
5
1.2
1.1
cb
subplot(2,1,1)
plot(t,ca)
xlabel('t')
ylabel('ca')
subplot(2,1,2)
plot(t,cb)
xlabel('t')
ylabel('cb')
0
1
0.9
116 /80
Debug
117 /80
Debug mode
breakpoint
Run vdv_sol_ae.m
118 /80
Debug mode (cont’d)
Step
Continue
119 /80
Solving ODE using S-function
120 /80
Example – VDV reaction
Enter search term
Constant
S-function
121 /80
Command Line:
open sfuntmpl.m
sfuntmpl.m will appear
122 /80
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
123 /80
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];
124 /80
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 = [];
125 /80
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)
sys = [];
126 /80
Example – VDV reaction (cont’d)
1
2
4
3
127 /80
Example – VDV reaction (cont’d)
Command Line:
x0=[2;1.117]
128 /80
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)
xlabel('t')
ylabel('cb')
129 /80
MV step change
Step
130 /80
Command Line:
x0=[3;1.117]
Start Simulation
Reach New Steady
State
131 /80
Refine Simulink file
132 /80
Refine Simulink file (cont’d)
Create Subsystem
in mouse right click menu
133 /80
Refine Simulink file (cont’d)
Mask Subsystem
in mouse right click menu
134 /80
Refine Simulink file (cont’d)
Add
135 /80
Refine Simulink file (cont’d)
Right double click on subsystem will appear the following window:
136 /80
Simulating LTI model
137 /80
LTI model 的建立
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)
[ys,ts]=step(vdv_tf);
plot(ts,ys)
0.6
0.5
0.4
0.3
0.2
0.1
0
-0.1
0
0.5
1
1.5
2
1.5
2
2.5
3
3.5
4
[yi,ti]=impulse(vdv_tf);
plot(ti,yi)
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
-1.2
0
0.5
1
2.5
3
3.5
4
pole(vdv_tf)
tzero(vdv_tf)
138 /80
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))
0
10
-1
10
-2
10
-1
10
0
10
1
10
2
10
400
300
200
100
[Gm,Pm,Wg,Wp]=margin(mag,phase,w)
0
-1
10
0
10
1
10
2
10
nyquist(kc*vdv_tf)
139 /80
Simulate LTI model using Simulink
140 /80
LTI model 的建立
Transfer Fcn
Nonlinear
Linear
1.19
1.18
1.17
1.16
cb
1.15
1.14
1.13
1.12
1.11
1.1
0
0.5
1
1.5
2
2.5
t
3
3.5
4
4.5
5
141 /80
Closed-loop simulation
142 /80
PID Controller
143 /80
Look Under
Mask
144 /80
Solving PDE
145 /80
Numerical Methods for PDE(pdepe)
ut  u xx
B.C. ffu (t , 0)  0, ffu (t ,1)  1
2x
I .C. ffu (0, x) 
1  x2
pdepe solves PDEs of the form
c( x, t , u, u x )ut  x
m
 m
( x b( x, t , u, u x ))  s ( x, t , u, u x )
x
with boundary conditions
p( xl , t , u)  q( xl , t ) b( xl , t , u, u x )  0
p( xr , t , u)  q( xr , t ) b( xr , t , u, u x )  0
146 /80
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 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);
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;
147 /80
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');
148 /80