四、常微分方程数值求解
1、常微分方程数值求解的一般概念
求解常微分方程初值问题就是寻找函数y(t)使之满足如下方程：
所谓其数值解法，就是求y(t)在离散节点tn处的函数近似值
的方法，
。相邻两个节点之间的距离称为步长。一般有单步法和多步法。
2、常微分方程数值求解函数
（1） [t,y]=solver(filename,tspan,y0,option)
其中，t和y分别给出时间向量和相应的数值解。solver为求常微分方程数值解的函数。filename为定义f(t,y)的函数名，该函数必须返回一个列向量。tspan形式为[t0,tf]，表示求解区间。y0是初始状态向量。option是可选参数，用于设置求解属性。常用的属性包括相对误差值RelTol
和绝对误差值AbsTol
（2）常微分方程数值求解函数的统一命名格式：
odennxx
- ode是Ordinary Differential Equation
- nn是数值，代表使用方法的阶数
- xx是字母，用于标注方法的专门特征
例子：求微分方程初值问题，并与精确解
f=@(t,y)(y^2-t-2)/4/(t+1); [t,y]=ode23(f,[0,10],2); y1=sqrt(t+1)+1; plot(t,y,'b:',t,y1,'r');
例子：已知一个二阶线性系统的微分方程为：
取a=2，绘制系统的时间相应曲线和相平面图。
对于高阶的微分方程，一般需要进行降阶处理.
则得到系统的状态方程：
f=@(t,x)[0,-2;1,0]*[x(1);x(2)]; [t,x]=ode45(f,[0,20],[1,0]); subplot(1,2,1); plot(t,x(:,2));%x(:,2)就是x(2),x(2)=x subplot(1,2,2); plot(x(:,2),x(:,1));
3、刚性问题
刚性问题是指一类微分方程，其解的分量有的变化很快，有的变化很慢，且相差悬殊。
刚性问题的数值解算法必须取很小步长才能获得满意的结果。解决刚性问题需要专门的方法，虽然非刚性算法也可以求解，但是需要很长的计算时间。
例子：
- 当lamda=0.01时
lamda=0.01; f=@(t,y)y^2-y^3; tic;%启动计时器 [t,y]=ode45(f,[0,2/lamda],lamda); toc;%结束计时器 disp(['ode45计算的点数' num2str(length(t))]);
历时 0.017177 秒。
ode45计算的点数157
- 当lamda=1e-5时
lamda=1e-5; f=@(t,y)y^2-y^3; tic;%启动计时器 [t,y]=ode45(f,[0,2/lamda],lamda); toc;%结束计时器 disp(['ode45计算的点数' num2str(length(t))]);
历时 0.392539 秒。
ode45计算的点数120565
这时计算时间明显加长，计算的点数剧增，表明此时常微分方程表现为刚性
- 对于刚性问题，我们选择以“s”结尾的函数
lamda=1e-5; f=@(t,y)y^2-y^3; tic;%启动计时器 [t,y]=ode15s(f,[0,2/lamda],lamda); toc;%结束计时器 disp(['ode15s计算的点数' num2str(length(t))]);
历时 0.770228 秒。
ode15s计算的点数98
很明显，计算时间大大缩短，计算的点数大大减少。
五、常微分方程应用举例
1、Lotka-Volterra模型
（1）模型分析
狐狸的增长模型如下：
（2）第一问：兔子数量初始值
，狐狸数量初始值
。
rabbitFox=@(t,x)[x(1)*(2-0.01*x(2));x(2)*(-1+0.01*x(1))]; [t,x]=ode45(rabbitFox,[0,30],[300,150]); subplot(1,2,1); plot(t,x(:,1),'-',t,x(:,2),'-*'); legend('x1(t)','x2(t)'); xlabel('时间');ylabel('物种数量'); grid on subplot(1,2,2); plot(x(:,1),x(:,2)); grid on
结论：兔子和狐狸两者相互制约，当狐狸逐渐增多时，兔子逐渐减少。当狐狸数量增加到一定数量时，由于种群内部竞争激烈，又导致其数量减少，从而兔子天敌减少，导致兔子数量增加。当兔子数量增加到一定数量时，由于种群内部竞争激烈，又导致其数量减少。如此循环，相互制约发展。
（3）第二问：兔子数量初始值
，狐狸数量初始值
。
（4）第三问：兔子数量初始值
=102，狐狸数量初始值
=198。
（5）验证（1/λ，2/λ）是稳定平衡点
- 取λ=0.01，所以稳定平衡点（1/λ，2/λ）即（100，200），以此点作为初值画图。
- 进一步分析初始值变化之后的图像
①当初始值变为(70,150)，向下偏离平衡点比较远，其图像：
可发现兔子和狐狸的数量变化比较剧烈，但是依然成周期性变化。
②当初始值变为(900,1600)，向上偏离平衡点很远，其图像：
可发现，兔子和狐狸的数量变化十分剧烈，但是依然成周期性变化。
2、Lotka-Volterra改进模型
（考虑环境的承受力，限制兔子的数量上限）
（1）第一问：原模型下，狐狸和兔子数量的时间函数曲线
rabbitFox=@(t,x)[x(1)*(2-0.01*x(2));x(2)*(-1+0.01*x(1))]; [t,x]=ode45(rabbitFox,[0,50],[300,150]); plot(t,x(:,1),t,x(:,2)); legend('x1(t)-兔子','x2(t)-狐狸'); xlabel('时间');ylabel('物种数量'); title('原模型下，狐狸和兔子数量的时间函数曲线') grid on
（2）第二问：改进模型下，狐狸和兔子数量的时间函数曲线（令R=400，λ=0.01）
rabbitFox=@(t,x)[2*x(1)*(1-x(1)/400-0.01*x(2));x(2)*(-1+0.01*x(1))]; [t,x]=ode45(rabbitFox,[0,50],[300,150]); plot(t,x(:,1),t,x(:,2)); legend('x1(t)-兔子','x2(t)-狐狸'); xlabel('时间');ylabel('物种数量'); title('改进模型下，狐狸和兔子数量的时间函数曲线') grid on
原模型无论经历多长时间，狐狸和兔子的数量总是在自己的平衡点之间波动。而改进之后的模型，虽然在前一段时间有较大的波动，但随着时间的推移，波动越来越小。在经历足够长的时间后，狐狸和兔子的数量分别达到稳定平衡。这更加接近自然界的实际情况。
（3）第三问：原模型下，狐狸数量相对于兔子数量的函数曲线。
rabbitFox=@(t,x)[x(1)*(2-0.01*x(2));x(2)*(-1+0.01*x(1))]; [t,x]=ode45(rabbitFox,[0,30],[300,150]); plot(x(:,1),x(:,2)); xlabel('时间');ylabel('物种数量'); grid on title('原模型下，狐狸数量相对于兔子数量的函数曲线')
（4）第四问：改进模型下，狐狸数量相对于兔子数量的函数曲线。
rabbitFox=@(t,x)[2*x(1)*(1-x(1)/400-0.01*x(2));x(2)*(-1+0.01*x(1))]; [t,x]=ode45(rabbitFox,[0,50],[300,150]); plot(x(:,1),x(:,2)); grid on xlabel('时间');ylabel('狐狸数量'); title('改进模型下，狐狸数量相对于兔子数量的函数曲线。') grid on
总结
