matlab-m文件常用积分函数-ode45含有时变参数用法/菜鸟理解4

写在前面

本人大四狗一名，最近在帮实验室肝项目，毕设用的强化学习暂且放下了一段时间，所以没有更新。在给实验室打工的过程中，遇到了一个需要用到时变参数的微分方程组，解决这种问题利用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.博主还会不定时更新的，喜欢的话就收藏吧。