Matlab Ode

Only available on StudyMode
  • Topic: Trigraph, Gh, Numerical ordinary differential equations
  • Pages : 12 (987 words )
  • Download(s) : 137
  • Published : October 9, 2012
Open Document
Text Preview
MATLAB
數值微積分與微分方程式求解

數值積分


b

a

f ( x) dx 等於由界限範圍 x = a 到 x = b 之間

曲線 f(x) 底下的面積

(a) 矩形以及 (b) 梯形
矩形以及
梯形
數值積分的圖解說明

數值積分
已知數據點的積分,
已知數據點的積分
已知數據點的積分,不知函數 f(x):trapz
I = trapz(x, y)

(梯形積分法)

x : 數據點之 x 值所構成的向量
y : 數據點之 f(x) 值所構成的向量

Ex:
>>x=[0 10 20 30 40];
>>y=[0.5 0.7 0.9 0.6 0.4];
>>area=trapz(x,y) %梯形法
area =
26.5000

數值積分
已知函數
之形式:
之形式
已知函數 f(x) 之形式:quad , quadl
I = quad(@fun, a, b)

(適應性辛普森法)

I = quadl(@fun, a, b)

(羅伯特二次式)

fun:定義函數的 function m-file 檔名
a :積分下限
b :積分上限

數值積分
Ex:



1

0

e − x cos( x) dx

1. edit fun.m
function y=fun(x)
y=exp(-x).*cos(x);
2. 求積分(回到Matlab Command Window)
求積分
area=quadl(@fun,0,1)
亦可使用
area=quadl(‘exp(-x).*cos(x)’,0,1)
NOTE: 函數內之數學運算必須使用向量個別元素之運算 函數內之數學運算必須使用向量個別元素之運算 (.* ./ .^)
(註:比較此結果與利用trapz指令計算之結果)

數值微分
已知數據點的微分
已知數據點的微分

在 x2 之微分

數值微分
可利用
可利用 diff 函數

Ex:
>>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)

數值微分
Ex: f ( x) = sin( x), f ′( x) = ? x ∈ [ 0, π ]
x = linspace(0,pi,20);
y = sin(x);
d = diff(y)./diff(x); % backward or forward difference
dc = (y(3:end)-y(1:end-2))./(x(3:end)-x(1:end-2)); % central difference dy = cos(x); % 實際微分值
plot(x, dy, x(2:end), d,'o', x(1:end-1), d,'x', x(2:end-1), dc,'^') xlabel('x'); ylabel('Derivative')
1
0.8
0.6
0.4
Derivative

>>
>>
>>
>>
>>
>>
>>

0.2
0
-0.2
-0.4
-0.6
-0.8
-1

0

0.5

1

1.5

2
x

2.5

3

3.5

工程問題中常微分方程式的解

常微分方程式
常微分方程式之形式:
dy
= f (t , y)
dt
一般解之形式:

yi +1 = yi + φ h
φ 是斜率 增量函數 (increment function),被用來自舊值yi 斜率或增量函數
斜率
外插到新值yi+1,h為步長大小 (step size)。
步長大小

此方法稱為單步方法 (one-step method),因為增量函數的 單步方法
值是根據單一點 i 的資訊。

歐拉法
在ti處的一次微分形式提供了直接的斜率估計: dy
= f (t i , yi )
dt ti

φ = f ( ti , yi )
yi +1 = yi + f ( ti , yi ) h

此公式即稱為歐拉法
歐拉法
(Euler's method) 。新的值y
利用斜率(等於在t處的一
階微分)在步長大小h上做
線性外插來得到預測值。

倫基-庫達法
倫基-
倫基
倫基-庫達(RK)法(Runge-Kutta methods)可以達到泰勒級數 方式的正確度,而不需計算高階的微分。

φ = a1k1 + a2 k2 + ⋯ + an kn
φ稱為增量函數,代表整個區間的斜率。
n 表示項數(階數),a是常數,p和q也都是常數。k之間是遞 迴的關係。
k1 = f (t i , yi )
k2 = f (t i + p1h, yi + q11k1h)
k3 = f (t i + p2 h, yi + q21k1h + q22 k2 h)

kn = f (t i + pn−1h, yi + qn−1,1k1h + qn−1,2 k2 h + ⋯ + qn−1,n−1kn−1h)

古典四階倫基-庫達法
最普遍使用的倫基-庫達法是古典四階RK法:
古典四階
1
yi +1 = yi + (k1 + 2k 2 + 2k3 + k 4 ) h
(4個斜率之加權平均)
6

k1 = f ( ti , yi )
h
h

k2 = f  ti + , yi + k1 
2
2

h
h

k3 = f  ti + , yi + k2 
2
2

k4 = f ( ti + h, yi + k3h )

方程式系統
許多實際的工程及科學問題需要求解的是聯立常微分方 程式系統,而不只是單一方程式。
dy1
= f1 (t , y1 , y2 ,⋯, yn )
dt
dy2
= f2 (t , y1, y2 ,⋯, yn )
dt

dyn
= fn (t , y1, y2 ,⋯, yn )
dt

求此系統的解需要n個在開始值t 時已知的起始條件。 對於單一方程式之數值方法皆可推展:對於每一個方程 式都使用單步方法,然後再進行下一步。

適應性倫基-庫達法
自動調整步長大小可以避免過
度使用計算量,因此有很大的
好處。
因為此方法對於解的軌跡具有
「適應性」,所以又稱為適應
適應
性步長大小控制 (adaptive stepsize control)。
要利用這個方法,每一步皆需
要估計局部截斷誤差。此誤差
估計可以當作縮短或增長步長
大小的依據。

適應性方法
有兩種主要方法可以合併使用適應性步長大小控制與 單步方法。
步長減半
步長減半 (step halving) 的方式是每一步做兩次,一次 是一整步,接下來是兩個半步。這兩者之間的差異代 表局部截斷誤差的估計值,我們根據這個誤差估計值 調整步長大小。
第二個方式稱為嵌入式RK法 (embedded RK method),
嵌入式
局部截斷誤差的估計值是根據兩個不同階數RK法預測 之間的差異。這是目前主流的使用方式,因為此方法 比步長減半法更有效率。

Matlab求解常微分方程式
dy (t )
= f ( t , y ) : ode45, ode23
dt
[t, y] = ode45(@fun, tspan, y0)

求解微分方程式
求解微分方程式

t:時間 (行向量)
y:解所形成之矩陣(每一行為一個應變數,每一列都對應到 行向量 t 中的時間)
fun:描述 ode function m-file 檔檔名
tspan:積分時間 ([起始時間 結束時間])
y0:初始值

ode定義檔格式 : 輸入為 t 與 y ,而輸出為表示 dy/dt 的行向量
function dydt = fun(t, y)
dydt(1) = < Insert a function of t and/or y here. >
dydt(2) = …………..
dydy = dydt’; % 轉換成行向量

Matlab求解常微分方程式
Ex:  dy1 = y + y e −t
1
2


 dt

 dy2 = − y y + cos(t )
12
 dt


y1 (0) = 0,

y2 (0) = 1

1. edit fun.m
function dydt=fun(t,y)
dydt(1) = y(1)+y(2)*exp(-t);
dydt(2) = -y(1)*y(2)+cos(t);
dydt = dydt’;
2.求解ode45 ode23 % 回到matlab command window
回到
[t,y]=ode45(@fun, [0 10],[0 1]) %[0 10]: t上下限, [0 1]: x的起始值 y1=y(:,1);
y2=y(:,2);...
tracking img