function
varargout = ode45(ode,tspan,y0,options,varargin)
%ODE45  Solve non-stiff differential equations, medium order method.
%   [TOUT,YOUT]
=
ODE45(ODEFUN,TSPAN,Y0)
with
TSPAN
=
[T0
TFINAL]
integrates
%   the system of differential equations
y' = f(t,y) from time T0 to TFINAL
%   with initial conditions Y0. ODEFUN is a function handle. For a scalar
T
%   and a vector Y, ODEFUN(T,Y) must return a column vector corresponding
%   to f(t,y). Each row in the solution array YOUT corresponds to a time
%   returned in the column vector TOUT.  To obtain solutions at specific
%   times T0,T1,...,TFINAL (all increasing or all decreasing), use TSPAN =
%   [T0 T1 ... TFINAL].
%
%   [TOUT,YOUT] = ODE45(ODEFUN,TSPAN,Y0,OPTIONS) solves as above with
default
%   integration
properties
replaced
by
values
in
OPTIONS,
an
argument
created
%   with
the
ODESET
function.
See
ODESET
for
details.
Commonly
used
options
%   are
scalar
relative
error
tolerance
'RelTol'
(1e-3
by
default)
and
vector
%   of
absolute
error
tolerances
'AbsTol'
(all
components
1e-6
by
default).
%   If certain components of the solution must be non-negative, use
%   ODESET to set the 'NonNegative' property to the indices of these
%   components.
%
%   ODE45 can solve problems M(t,y)*y' = f(t,y) with mass matrix M that is
%   nonsingular.
Use
ODESET
to
set
the
'Mass'
property
to
a
function
handle
%   MASS
if
MASS(T,Y)
returns
the
value
of
the
mass
matrix.
If
the
mass
matrix
%   is constant, the matrix can be used as the value of the 'Mass' option.
If
%   the
mass
matrix
does
not
depend
on
the
state
variable
Y
and
the
function
%   MASS is to be called with one input argument T, set 'MStateDependence'
to
%   'none'.
ODE15S
and
ODE23T
can
solve
problems
with
singular
mass
matrices.
%
%   [TOUT,YOUT,TE,YE,IE]
=
ODE45(ODEFUN,TSPAN,Y0,OPTIONS)
with
the
'Events'
%   property in OPTIONS set to a function handle EVENTS, solves as above
%   while also finding where functions of (T,Y), called event functions,
%   are zero. For each function you specify whether the integration is
%   to terminate at a zero and whether the direction of the zero crossing
%   matters. These are the three column vectors returned by EVENTS:
%   [VALUE,ISTERMINAL,DIRECTION]
=
EVENTS(T,Y).
For
the
I-th
event
function:
%   VALUE(I)
is
the
value
of
the
function,
ISTERMINAL(I)=1
if
the
integration
%   is to terminate at a zero of this event function and 0 otherwise.
%   DIRECTION(I)=0
if
all
zeros
are
to
be
computed
(the
default),
+1
if
only
%   zeros
where
the
event
function
is
increasing,
and
-1
if
only
zeros
where
%   the event function is decreasing. Output TE is a column vector of times
%   at which events occur. Rows of YE are the corresponding solutions, and
%   indices in vector IE specify which event occurred.
%
%   SOL = ODE45(ODEFUN,[T0 TFINAL],Y0...) returns a structure that can be
%   used with DEVAL to evaluate the solution or its first derivative at
%   any point between T0 and TFINAL. The steps chosen by ODE45
are returned
%   in a row vector SOL.x.  For each I, the column SOL.y(:,I) contains
%   the
solution
at
SOL.x(I).
If
events
were
detected,
SOL.xe
is
a
row
vector
%   of points at which events occurred. Columns of SOL.ye are the
corresponding
%   solutions, and indices in vector SOL.ie specify which event occurred.
%
%   Example

目录
写在前面
ode45积分器
带有时变参数的ode45积分
总结
写在后面

写在前面
本人大四狗一名，最近在帮实验室肝项目，毕设用的强化学习暂且放下了一段时间，所以没有更新。在给实验室打工的过程中，遇到了一个需要用到时变参数的微分方程组，解决这种问题利用simulink很好解决，但是项目要求使用m文件进行编程，在以往的学习中，m文件的积分函数一般就是使用ode45，变步长积分器。常用的语法是
[t,y] = ode45(@odefun,tspan,y0)
[t,y] = ode45(@odefun,tspan,y0,options)
[t,y,te,ye,ie] = ode45(@odefun,tspan,y0,options)
sol = ode45(___)

这里面有积分时间域，初值以及积分器的设置，我们要做的就是把微分方程组输入到odefun中然后利用ode45进行积分，积分的过程相对是一个黑盒子，在我以往的学习过程中没有把黑盒子拆开了看明白，在解决时变微分方程组时就遇到了问题，我经过搜索相关问题，没有发现网络上给出很好的解答，于是把我研究的一点小东西写下来，给各位同仁分享，希望大家共同进步，另外，如果这边文章对您的工作生活有所用处，还是请您给个点赞关注评论。
ode45积分器
对matlab有所了解的同志应该知道，matlab是一款除了不能生孩子都可以的软件（希望中国能早日有自己的类似软件。），在解决微分方程组的问题时常用simulink搭积木解决，使用m文件相对就会不直观，不容易抓住各种关节，调试起来比较麻烦，但是理解了ode积分器的工作原理之后，也可以很好的进行debug工作。ode45积分器是一个变步长积分器，对很多问题都能很好的解决，我们的抓手其实就是odefun，编辑好微分函数后，具体的积分我们就管不到了，但是在odefun中我们也可以进行一定程度的debug，比较常用的就是把一个语句后面的分号；去掉，来观察我们要看的语句到底有没有运转，运转结果如何，结合我们的猜想进而更改我们的函数。ode45函数具体的使用方法可以在math work中国看他官方给出的help文件，我在此诸摘录一点，
function dydt = odefcn(t,y)
dydt = zeros(2,1);
dydt(1) = y(2);
dydt(2) = 5*y(1);

这就是我们说的odefun，用一个类似状态观测器的形式给出，我们要把这个函数写好，放在matlab可以找到文件夹中。在调用时可以使用
tspan = [0 5];
y0 = [0 0.01];
[t,y] = ode45(@odefcn, tspan, y0);

带有时变参数的ode45积分
时变参数和普通的参数的区别就是他会跟着时间变化，同时ode45是一个变步长积分器，我们怎么知道积分的每一步时间是多少呢?需要变化的参数怎么传递进去呢？怎么把每一步的时间对应到我们需要用到的参数具体值呢？这就需要在odefcn中找到一个抓手，在mathwork中国网站中，给出了时变参数的相关说明，我在此列下来：
ft = linspace(0,5,25);
f = ft.^2 - ft - 3;

gt = linspace(1,6,25);
g = 3*sin(gt-0.25);

function dydt = myode(t,y,ft,f,gt,g)
f = interp1(ft,f,t); % Interpolate the data set (ft,f) at time t
g = interp1(gt,g,t); % Interpolate the data set (gt,g) at time t
dydt = -f.*y + g; % Evaluate ODE at time t
end

tspan = [1 5];
ic = 1;
opts = odeset('RelTol',1e-2,'AbsTol',1e-4);
[t,y] = ode45(@(t,y) myode(t,y,ft,f,gt,g), tspan, ic, opts);

从上面中的例子中我们可以总结出，怎么把参数带进odefcn，用到的语句就是
[t,y] = ode45(@(t,y) myode(t,y,ft,f,gt,g), tspan, ic, opts);

和一般的用法不同，这里我们@的是（t，y）说明函数的时间、初值为t，y，其他的输入是参数，ode45积分器不会使用，这里的参数使用就和其他函数一样了。

！！！！！！！！！！！！！！！！！！！！！

接下来最重要的一点，就是我们在ode45中的抓手，就是（t，y）中的t，光看官网给出的例子，interp1是差分，在ft，f中间差分到t的值，但是我们还是不知道t究竟是一个值还是一个数组，因为我们输入的时间范围，是一个数组。如果interp1差分出来的是积分所有步骤的参数值也说得通，我在最开始也是这么理解的，但是这里的t其实本身就是时变的，并不是一个大数组（和输出的t数组大不同！）在odefcn中循环调用的t就是一个数，代表着本次积分步骤中时间变量。

有了这个理解我们就可以很好地解决时变参数微分方程组积分的问题了，在这里我给出两个解法
1：利用例子中给出的方法，把要使用的参数的一个时间序列输入到odefcn中，在使用中用matlab自身的插值函数来应对变步长不确定的t，这样的方法可能会导致函数精确度下降，因为本身插值就是有误差的，在积分过程中，误差是会累积的，最终形成一个较大误差。
2.如果我们的时变参数不是查表获得，也就是不需要用到插值，也就是有更好的表达方式来获得更加精确的参数值，这种情况下，我们就可以把需要用到的其他参数输入或者global一下变成全局变量，在odefcn中自我求解，这种方式的准确性较高。
我在项目中用到的就是第二种方法，由于项目本身涉密，函数文件不在此展示。
总结
[dy]=odefcn（t,y,a,b）
.........
end

中，t是一个时变数，大小等于积分步骤正在进行的时间大小，y是等式左边的初值，a，b是可以带进来的参数，dy就是y的一阶导数。如果我们需要解决高阶微分方程组问题，可以把ddy成为dy的一阶导，具体思想和状态空间法类似，也和simulink中搭建高阶微分方程组的思想类似。
在得到这个抓手之后，大家就可以根据自身情况来改变函数的内部结构，不把他当作一个黑盒子，拆开来，走进去。
写在后面
由于网上没有这样的详细说明，大佬们都不屑于解释太多，就只能由我这样的菜鸟谈谈想法，有什么理解不到位的地方欢迎和我多交流，共同进步。

水平有限，请有幸看到的您多多指正。在工作之余，祝您早安，午安，晚安。Have a nice day。

（可能的话点个赞再走吧）

P.S.博主还会不定时更新的，喜欢的话就收藏吧。

上篇博客介绍了Matlab求解常微分方程组解析解的方法：博客地址

微分方程组复杂时，无法求出解析解时，就需要求其数值解，这里来介绍。

以下内容按照Matlab官方文档提供的方程来展开（提议多看官方文档）
介绍一下核心函数ode45()
一般形式：[t,y] = ode45(odefun,tspan,y0)		其中 tspan = [t0 tf]

功能介绍：求微分方程组 y′=f(t,y) 从 t0 到 tf 的积分，初始条件为 y0。解数组 y 中的每一行都与列向量 t 中返回的值相对应。
1. 一阶微分方程求解(简单调用即可)
方程：y’=2*t

代码：
tspan=[1 6]; %定义自变量x的取值空间为1-6
y0=0;%定义因变量的初值，当x=1（x取值空间的第一个数）时，y0=0
[t,y]=ode45(@(t,y) 2*t,tspan,y0); %定义函数y'=2*t,使用ode45求解
plot(t,y,'-o'); %绘制求得的数值曲线

说明：简单的odefun参数就是这个形式，@(x,y) fun	中间是空格，不能用逗号

结果：


2. 二阶微分方程求解(引入函数文件)
方程：范德波尔方程 y1’’-u(1-y1²)*y1’+y1=0;（这里设u=1）

代码：

定义输入的方程，以函数形式定义
function dydt=odefun(t,y)
%二阶方程为y1''-(1-y1²)*y1'+y1=0;
%降阶为两个方程：y1'=y2; 
%               y2'=(1-y1²)*y2-y1;
%t虽然没有使用，但必须要作为参数写入
dydt=[y(2);(1-y(1)^2)*y(2)-y(1)];

求解作图
tspan=[0 20]; %定义自变量x的取值空间为0-20
y0=[2;0];%定义因变量的初值，当x=0时，y1=2,y2=y1'=0;
[t,y]=ode45(@odefun,tspan,y0); %使用ode45求解
%%下面为作图过程，不解释
plot(t,y(:,1),'-o',t,y(:,2),'-o')
title('Solution of van der Pol Equation (\mu = 1) with ODE45');
xlabel('Time t');
ylabel('Solution y');
legend('y_1','y_2')

说明：先进行降阶，使y2=y1’，这时候y2’=y1’’，然后把两个方程写成列向量形式就行了

结果：


3. 求解微分方程组（和2类似）
这里就和求解二阶方程类似的，只不过不需要降阶，仍旧需要一个函数来定义方程组。我们这里不用官方文档的例子，用同学的循坏摆问题来进行演示。

方程：

给定的初值（w接近0，但实际上不能设置为0）：



代码：

定义输入的方程
function dRvw=func(t,Rvw)
%% 函数功能：为ode45提供微分方程
%输入：t:时间序列，就是θ；Rvw:因变量,Rvw(1)代表R,Rvw(2)代表v,Rvw(3)代表w
%输出：dRvw:因变量的一阶微分,dRvw(1)代表dR,dRvw(2)代表dv,dRvw(3)代表dw
%% 初始化因变量的一阶微分，3×1的向量
dRvw=zeros(3,1);
%% 参数初始化
r=0.01;u=0.1;g=9.8;M=10;m=1;
%% 输入微分方程式
dRvw(1)=-Rvw(2)/Rvw(3)-r;
dRvw(2)=(M*g-(m*Rvw(3)^2*Rvw(1)+m*g*sin(t)) *exp(u*(t+pi/2)))/(Rvw(3)*(M+m));
dRvw(3)=Rvw(3)*r/Rvw(1)+2*Rvw(2)/Rvw(1)+g*cos(t)/(Rvw(3)*Rvw(1));

求解作图
%% 参数初始化 start_Theta是θ的初始值 end_Theta是θ的结束值
%R是半径初值；v是线速度初值；w是角速度初值
start_Theta=0;end_Theta=2*pi;R=1;v=0;w=1e-5;
%% 使用ode45方法计算微分方程组func的数值解
%func是带有方程组的函数 
%[start_Theta end_Theta]是自变量范围
%[R;v;w]是方程初值
%T是自变量的数组，Rvw是对应的因变量的数值。Rvw(:,1)=R;Rvw(:,2)=v;Rvw(:,3)=w;
[T,Rvw]=ode45(@func,[start_Theta end_Theta],[R;v;w]);
%% 组合自变量和因变量。TRvw(:,1)=θ;TRvw(:,2)=R;TRvw(:,3)=v;Trvw(:,4);
TRvw=[T,Rvw];
%% 绘制自变量-因变量图像.figure1是R-θ,figure2是v-θ,figure3是w-θ
plot(T,Rvw(:,1),'-o',T,Rvw(:,2),'-o',T,Rvw(:,3),'-o')
title('自变量-因变量图像');
xlabel('Angle θ');
ylabel('Solution');
legend('R','v','w')

说明：注释的应该是比较清楚的，把三个方程写成列向量的形式就行
PS:有些人和我说不能运行，然后我看了他们出错的原因，有点儿哭笑不得。出错的基本上都是运行上面的dRvw=func（t，Rvw）这个函数的。说明一下，这是有参数的函数，不给参数不能直接运行的。下面的求解作图脚本才是需要运行的哈，它调用了函数，才得到的结果。如果大家还发现什么问题，欢迎私戳或评论。
PS+ 有了PS之后，还是很多人问我参数的问题，我在这里直接把文件给大家：cupt.zip

提取码：6k8n。

解压后，有两个文件，大家直接运行cupt.m即可。
结果：


4. 更多形式
讲到这里，大部分我们用到的微分方程形式都可以求解了，Matlab还支持带有时变项和额外参数的微分方程求解，这里不再赘述，大家可以自行参阅官方文档。

ode45解微分equation
[t,y] = ode45(odefun,tspan,y0)example
[t,y] = ode45(odefun,tspan,y0,options)example

首先编写函数文件.m作为ode45的第一个参数odefun
function y = fun(t,x0,flag,k)
%flag is a [] to mark the args 可略去不写y = fun(t,x0,k) %kequ中参数
%x0初值given
%---------------
%some global setting
lamda=1;
%---------------
    y=[x0(2)*k;
    lamda*x0(2)*(1-x0(1)-x0(2))-k*x0(2)];
end

调用：
[t,x]=ode45('fun',tspan,[0,tf],[],params);
%or [t,x]=ode45(@fun,tspan,[0,tf],[],params);

tspan为[0:10] 自变量范围

第三个参数为初值向量

第四个参数空阵组作flag 区别开传入fun的参数和ode45所需要的参数