连续型模型-北太天元学习30
在"北太天元学习18"中，我们以离散的时间步长对现象进行了建模。今天我们学习连续型模型，它在连续的时间尺度上跟踪正在发生的事情。连续模型通常采用微分方程的形式，而不是差分方程。由于这些微分方程中的许多都没有解析解，因此也有必要理解数值方法。
我们举一个物理上的例子，一个质点在下落的过程中收到重力和阻力，其中阻力大小和速度成正比， 方向与速度方向相反。 我们使用牛顿第二定律建立数学模型
设质点的质量是m, 在0时刻速度是 v0, 我们可以得到下面的方程
v'(t) = g - k v(t)
其中g 是重力常数，k 是阻力系数。
这是一个常微分方程， 要求的未知量(unknow)是函数 v(t)。
这是一个连续型的模型，不再像我们在“北太天元学习18”中那样仅仅给出了速度在两个不同时刻的关系，而是使用了导数，给出速度在每个时刻的变化率（对时间的导数) 和其它的量之间的关系。
我们以前学习的 3x+2= 5 这样的一元一次方程，这是一个代数方程，它的未知数不是一个函数, 而仅仅是一个数。 再看一个方程: 3x(t)+2 = 5*cos(t), 这个方程的未知数x(t)是一个函数， 但是我们求解这个方程和求解 3x+2=5 相比 一点都不困难：例如t=0, 我们轻松地求解 x(0), 在求解的过程中，会看到 x(0) 的求解 和 t=1,t=1.5 等t在别的地方的值是没有关系的，因此这还是一个代数方程。
但是再看我们刚才说的常微分方程: v'(t) = g- k v(t) 要求出 v(t) 显然不仅仅依赖于 t 在一点的值，也就是说t的所有点信息是耦合在一起的。
简单地说, 常微分方程的定义包含了三个要素，
它是一个方程,
它的unknown(未知量)是一个一元函数（只有一个自变量），
它涉及到这个函数的导数或者高阶导数。
在"北太天元学习18"，我们用离散模型研究人口问题，我们给出了两个时刻的人口之间的关系 x_{n+1} = x_n + k *x_n
这被称为离散马尔萨斯人口模型(discrete Malthusian population model)，表示两个时刻的人口变化量和上一个时刻的人口成正比，比例系数是k. 我们可以假设人口增长率和当前时刻的人口成正比，这里的人口增长率是指人口对时间的导数, 这样就给出了马尔萨斯人口模型， 是一个常微分方程
x'(t) = r * x(t),
其中 r 是增长常数。
这个方程是可以写出解析解的，但是很多常微分方程是无法解析求解，需要进行数值求解。我们可以把 x'(t) = r*x(t) 离散成差分方程来求解，例如
假设我们对 x(t) , 0<=t <= T 感兴趣，那么我们可以把[0,T]分成N 个小区间,其中小区间的端点x_n = n h, h=T/N 是时间步长，n=0,1,...N. 我们寻找 x(t) 在t_n 上近似值 x_n, 因为我们是从 x'(t) = r*x(t) 出发找到了 x_n 满足的一个近似的方程
( x_{n+1} - x_{n} ) / h = r* x_n
这个被称为差分方程，因为导数又被称为微商（微分之商), 我们刚刚得到的这个差分方程和微分方程的区别是把微商 x'(t_n)换成了差商 ( x_{n+1} - x_{n} ) / h .
差分方程就成了离散模型，是我们在“北太天元学习18”中得到的。
基于这个讨论，我们就得到了求解常微分方程的一种数值方法，被称为欧拉法。我们可以用北太天元写一个代码来求微分方程数值解，实际上就是和我们在北太天元学习18中的代码是一样的，这里就不多重复了。
欧拉法是求解常微分方程的一种方法，还有很多其它的方法。这些方法还可以应用到
微分方程组的求解. 在“北太天元学习18”中提到的捕食者-被捕食者模型对应的连续模型是一个常微分方程组: 用 x(t) 表示兔子在t时刻的数量，用 y(t) 表示狼在t时刻的数量，
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
这些函数的用法是类似的，我们以ode45为例子讲讲如何使用函数求常微分方程组的数值解。
% 求解二阶ODE方程初值问题
% 其中未知函数z(t)的满足的ODE方程是
% z'' = (1-z^2)*z' -z
% 初值条件是
% z(0) = 2
% z'(0) = 0
%
% 我们首先把二阶ODE写成一阶ODEs
% 引入函数 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有两列，分别对应函数和导数
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
还可以参考:
【情人节与微分方程组】 https://www.bilibili.com/video/BV1ZT411Q7Uj/?share_source=copy_web&vd_source=2adc5aa7a702b808eb8b31dbd210f954
