Matlab Ode

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 x

2

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之間是遞...