目录

    • 写在前面
  • 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
endtspan = [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-m文件常用积分函数-ode45含有时变参数用法/菜鸟理解4相关推荐

  1. Matlab命令集--常用字符串函数

    Matlab命令集--常用字符串函数 常用函数 eval  :运行字符串表示的表达式 char  :将数组变成字符串 double:将数字字符串变成数字 字符串操作 deblank :去掉字符串末尾的 ...

  2. MATLAB不定积分的运算,matlab中怎么把积分函数 int 得到的不定积分式代入 solve 函数中进行计算?...

    答:syms r x fun=int(r*exp(-2*(r/2)^2),r,0,x); x=solve(fun-0.5) x = 2^(1/2)*log(2)^(1/2) -2^(1/2)*log( ...

  3. matlab 命令文件转成函数文件,科学网—[转载]利用MATLAB将nc文件转成tif - 张乐乐的博文...

    参考链接:https://blog.csdn.net/yangjh1991/article/details/69788778 Lon = ncread(InFile,'lon'); %读取经度数据 L ...

  4. C++常用头文件——常用数学函数头文件

    cmath头文件 cmath是c++语言中的标准库头文件.其中的 "c" 表示其中的函数是来自 C标准库,"math"表示为数学常用库函数.此文件原作为< ...

  5. matlab算定积分时积分函数有变量,MATLAB中计算定积分时可否将一个函数作为积分变量?...

    共回答了17个问题采纳率:82.4% 提供两种解法供参考. 1.解析解法作变量置换t=1/x,则积分上限为1,下限为inf:>> syms x t >> f=(1+1/x)^x ...

  6. MATLAB实现变限积分函数的积分/ 多重积分/ 如何解决求积分显示AB浮点标量报错

    重点是要用arrayfun扩展 求变限积分的积分: fun_inner = @(r) r.^2; fun_integral = @(x) integral(fun_inner, 0,x); fun_o ...

  7. 常用小波基函数以及多尺度多分辨率的理解

    最近在做图像融合相关的工作.需要用到小波分析. 小波分析中用到的小波函数具有不唯一性,即小波函数具有多样性.使用不同的小波基分析同一个问题会产生不同的结果. 以下列出15种小波基函数是matlab支持 ...

  8. 常用小波基函数以及多尺度多分辨率的理解1

    小波分析中用到的小波函数具有不唯一性,即小波函数具有多样性.使用不同的小波基分析同一个问题会产生不同的结果. 以下列出15种小波基函数是matlab支持的15种: 小波函数 Harr Daubechi ...

  9. php文件操作(最后进行文件常用函数封装)

    文件信息相关API $filename="./1-file.php";//filetype($filename):获取文件的类型,返回的是文件的类型echo '文件类型为:',fi ...

  10. matlab中trapz,MATLAB中trapz和cumtrapz函数

    这两个函数都是MATLAB中的内置函数,是基于梯形法则的数值积分公式 例如我们有函数y=x^3-2x-3,为了计算在[0,1]上的积分,可以这么做: 其中x和y分别是自变量和对应的值,trapz其实就 ...

最新文章

  1. c语言int t格式,如何在C中打印int64_t类型
  2. 亲戚再也看不见我一个人食吉野家了
  3. C++cycle sort循环排序的实现算法(附完整源码)
  4. Keyword-Driven Testing
  5. 无法解析此远程名称: 'www.***.com' 解决办法(转)
  6. 通过点击事件监听 setOnClickListener 彻底理解回调-Android
  7. Android FTP Server 1
  8. c++隐式类型转换存在的陷阱
  9. Maven学习(八)-----Maven依赖机制
  10. InfoPath2010表单-IE浏览器2个“微型内嵌工具”的使用和介绍
  11. python中command是什么意思_python 设计模式之命令(Command)模式
  12. 错误代码: 1066 Not unique table/alias: #39;c#39;
  13. ubuntu mysql主从配置_MYSQL 主从数据库的配置 ubuntu 12.04
  14. 三维点云学习(4)5-DBSCNA python 复现-2-kd-_tree加速
  15. [bzoj1566][NOI2009]管道取珠
  16. ffmpeg API变更 2009-03-01—— 2017-05-09变更
  17. 支付宝app登录授权的infoStr授权登录流程
  18. 关于常用(?)字符串处理函数的合集
  19. 用谷歌注册Kaggle没有出现验证码的情况
  20. java 图片相似搜索_JAVA比较两张图片相似度的方法

热门文章

  1. 【SNE-RoadSeg 解读】结合表面法向量的路面分割网络(ECCV2020)
  2. APP测试与WEB测试
  3. DDR4、LPDDR4、LPDDR4x区别及DDR拓展
  4. Windows cmd快捷键
  5. 【Bye-Bye】MMD镜头+动作打包下载.zip
  6. 周志华《机器学习》个人笔记
  7. 计算机为什么不用三十二进制,32位进制导致2038年问题
  8. Hspice 反相器仿真
  9. js刷新页面的几种方式与区别
  10. Powermill数控编程培训,潇洒模具三步教您精通cnc数控编程