在"北太天元學(xué)習(xí)18"中，我們以離散的時間步長對現(xiàn)象進(jìn)行了建模。今天我們學(xué)習(xí)連續(xù)型模型，它在連續(xù)的時間尺度上跟蹤正在發(fā)生的事情。連續(xù)模型通常采用微分方程的形式，而不是差分方程。由于這些微分方程中的許多都沒有解析解，因此也有必要理解數(shù)值方法。

我們舉一個物理上的例子，一個質(zhì)點在下落的過程中收到重力和阻力，其中阻力大小和速度成正比， 方向與速度方向相反。 我們使用牛頓第二定律建立數(shù)學(xué)模型

設(shè)質(zhì)點的質(zhì)量是m, 在0時刻速度是 v0,? 我們可以得到下面的方程
?? ??? ? v'(t) =? g - k v(t)
其中g(shù) 是重力常數(shù)，k 是阻力系數(shù)。
這是一個常微分方程， 要求的未知量(unknow)是函數(shù) v(t)。
這是一個連續(xù)型的模型，不再像我們在“北太天元學(xué)習(xí)18”中那樣僅僅給出了速度在兩個不同時刻的關(guān)系，而是使用了導(dǎo)數(shù)，給出速度在每個時刻的變化率（對時間的導(dǎo)數(shù)) 和其它的量之間的關(guān)系。

我們以前學(xué)習(xí)的 3x+2= 5 這樣的一元一次方程，這是一個代數(shù)方程，它的未知數(shù)不是一個函數(shù), 而僅僅是一個數(shù)。 再看一個方程: 3x(t)+2 = 5*cos(t), 這個方程的未知數(shù)x(t)是一個函數(shù)， 但是我們求解這個方程和求解 3x+2=5 相比 一點都不困難：例如t=0, 我們輕松地求解 x(0), 在求解的過程中，會看到 x(0) 的求解 和 t=1,t=1.5 等t在別的地方的值是沒有關(guān)系的，因此這還是一個代數(shù)方程。

但是再看我們剛才說的常微分方程: v'(t) = g- k v(t)? 要求出 v(t) 顯然不僅僅依賴于 t 在一點的值，也就是說t的所有點信息是耦合在一起的。

簡單地說, 常微分方程的定義包含了三個要素，

1. ? 它是一個方程,

2. ?它的unknown(未知量)是一個一元函數(shù)（只有一個自變量），?

3. ?它涉及到這個函數(shù)的導(dǎo)數(shù)或者高階導(dǎo)數(shù)。
   
   在"北太天元學(xué)習(xí)18"，我們用離散模型研究人口問題，我們給出了兩個時刻的人口之間的關(guān)系?? x_{n+1} = x_n + k *x_n ?
   這被稱為離散馬爾薩斯人口模型(discrete Malthusian population model)，表示兩個時刻的人口變化量和上一個時刻的人口成正比，比例系數(shù)是k. 我們可以假設(shè)人口增長率和當(dāng)前時刻的人口成正比，這里的人口增長率是指人口對時間的導(dǎo)數(shù), 這樣就給出了馬爾薩斯人口模型， 是一個常微分方程
   ?? ??? ??? ??? ?x'(t) = r * x(t),
   其中 r 是增長常數(shù)。
   
   這個方程是可以寫出解析解的，但是很多常微分方程是無法解析求解，需要進(jìn)行數(shù)值求解。我們可以把 x'(t) = r*x(t) 離散成差分方程來求解，例如
   假設(shè)我們對 x(t) , 0<=t <= T 感興趣，那么我們可以把[0,T]分成N 個小區(qū)間,其中小區(qū)間的端點x_n = n h, h=T/N 是時間步長，n=0,1,...N. 我們尋找 x(t) 在t_n 上近似值 x_n,? 因為我們是從 x'(t) = r*x(t) 出發(fā)找到了 x_n 滿足的一個近似的方程
   ????? ( x_{n+1} - x_{n} ) / h =? r* x_n
   這個被稱為差分方程，因為導(dǎo)數(shù)又被稱為微商（微分之商), 我們剛剛得到的這個差分方程和微分方程的區(qū)別是把微商 x'(t_n)換成了差商? ( x_{n+1} - x_{n} ) / h .
   差分方程就成了離散模型，是我們在“北太天元學(xué)習(xí)18”中得到的。
   基于這個討論，我們就得到了求解常微分方程的一種數(shù)值方法，被稱為歐拉法。我們可以用北太天元寫一個代碼來求微分方程數(shù)值解，實際上就是和我們在北太天元學(xué)習(xí)18中的代碼是一樣的，這里就不多重復(fù)了。
   
   歐拉法是求解常微分方程的一種方法，還有很多其它的方法。這些方法還可以應(yīng)用到
   微分方程組的求解. 在“北太天元學(xué)習(xí)18”中提到的捕食者-被捕食者模型對應(yīng)的連續(xù)模型是一個常微分方程組: 用 x(t) 表示兔子在t時刻的數(shù)量，用 y(t) 表示狼在t時刻的數(shù)量，
   x'(t) = r*x(t) - f*x(t)*y(t),
   y'(t) = w*y(t) + g*x(t)*y(t).
   
   
   在北太天元中提供了很多方法求解常微分方程或者常微分方程組，例如
   ode113?? ?ode15i?? ?ode15ipdinit?? ?ode15ipdupdate?? ?ode15s?? ?ode23?? ?ode23s?? ?ode23t?? ?ode23tb?? ?ode45?? ?ode78?? ?ode89
   這些函數(shù)的用法是類似的，我們以ode45為例子講講如何使用函數(shù)求常微分方程組的數(shù)值解。
   
   % 求解二階ODE方程初值問題
   % 其中未知函數(shù)z(t)的滿足的ODE方程是
   %? z'' = (1-z^2)*z' -z
   % 初值條件是
   %? z(0) = 2
   %? z'(0) = 0
   % ?
   % 我們首先把二階ODE寫成一階ODEs
   % 引入函數(shù) y = [ z; z']
   % 得到 y 的一階ODEs
   %? y' = [ y(2) ; (1-y(1)^2)*y(2)-y(1)]
   % 然后用ode45求解 上面的一階級ODEs，
   % y(1) 就是所求的二階ODE初值問題的解
   clear, clc, close all;
   odefun = @(t,y) [ y(2); (1-y(1)^2)*y(2)-y(1)]; % 必須返回列向量
   tspan = [0 20];
   y0 = [2; 0];? % 這里看起來需要用列向量，但是實際上用行向量也能正確計算 ?
   options = odeset('RelTol',1e-6,'AbsTol',1e-8);
   [t,y] = ode45(odefun,tspan,y0,options); %計算出來的y有兩列，分別對應(yīng)函數(shù)和導(dǎo)數(shù)
   plot(t,y(:,1),'-o',t,y(:,2),'-*');
   legend('二階ODE的解z', 'z''', 'Location', 'SouthEast')
   xlabel('t')
   title("用ode45求解二階常微分方程的初值問題: z'' = (1-z^2)*z' -z? ; z(0)=2; z'(0)=0 ")
   
   詳情看
   【lecture18-1-ode_intro 北太天元上的常微分方程求解工具包的使用】 https://www.bilibili.com/video/BV1FG411A79H/?share_source=copy_web&vd_source=2adc5aa7a702b808eb8b31dbd210f954
   
   還可以參考:
   【情人節(jié)與微分方程組】 https://www.bilibili.com/video/BV1ZT411Q7Uj/?share_source=copy_web&vd_source=2adc5aa7a702b808eb8b31dbd210f954