MATLAB之ODE

MATLAB之ODE

一、简介

• 在工程和计算中常常需要求解微分方程，但微分方程的解析解通常很难得到，所以一种退而求其次的方法就是求其数值解。这里是使用Matlab中的ode45求解常微分方程数值解的简单说明。具体可使用"help ode45"命令查看它的帮助文档。
• ode45，常微分方程的数值求解。MATLAB提供了求常微分方程数值解的函数。
• Matlab中求微分方程数值解的函数有七个：ode45，ode23，ode113，ode15s，ode23s，ode23t，ode23tb。
• ode是Matlab专门用于解微分方程的功能函数。该求解器有变步长（variable-step）和定步长（fixed-step）两种类型。不同类型有着不同的求解器，其中ode45求解器属于变步长的一种，采用Runge-Kutta算法；其他采用相同算法的变步长求解器还有ode23。
• ode45表示采用四阶-五阶Runge-Kutta算法，它用4阶方法提供候选解，5阶方法控制误差，是一种自适应步长（变步长）的常微分方程数值解法，其整体截断误差为(Δx)^5。解决的是Nonstiff(非刚性)常微分方程。

二、用法

• [T,Y] = ode45(odefun,tspan,y0)
• odefun 是函数句柄，可以是函数文件名，匿名函数句柄或内联函数名
• tspan 是区间 [t0 tf] 或者一系列散点[t0,t1,…,tf]
• y0 是初始值向量
• T 返回列向量的时间点
• Y 返回对应T的求解列向量

三、举例

假如我们要求解如下方程组，它描述的是无外力作用下的刚体运动：

y_1'= y_2*y_3    y_1(0)=0
y_2'= -y_1*y_3    y_2(0)=1
y_3'= -0.51*y_1*y_2    y_3(0)=1

首先创建一个函数文件，文件名为odefun.m，代码如下：

% odefun.m
function dy = odefun(t,y)
dy = zeros(3,1);    % a column vector
dy(1) = y(2) * y(3);
dy(2) = -y(1) * y(3);
dy(3) = -0.51 * y(1) * y(2);
end

然后创建命令文件，文件名为testode45.m，代码如下：

% testode45.m
% 清空所有变量
clear
% 清空屏幕
clc

% 时间跨度取0-12，可以空格分隔，也可以用逗号分隔
tspan = [0,12];
% 初始值
y0 = [0,1,1];
% 调用语句
[T,Y] = ode45( @odefun, tspan, y0 );
% 绘图
plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.')
legend('x','y','z')

最后运行testode45.m文件，就可以得到所求方程的图像：