本节将介绍常微分方程初值问题的数值求解,主要内容分为三个部分:常微分方程数值求解的概念、求解函数及刚性问题。
一、常微分方程数值求解的一般概念
首先,凡含有参数,未知函数和未知函数导数 (或微分) 的方程,称为微分方程,有时简称为方程,未知函数是一元函数的微分方程称作常微分方程,未知函数是多元函数的微分方程称作偏微分方程。微分方程中出现的未知函数最高阶导数的阶数,称为微分方程的阶。
常微分方程数值求解,指研究求解初值问题各类数值方法的构造、理论分析与数值实现问题。研究的主要对象为一阶方程组初值问题:

其中,其中y=y(x)是未知函数,y(x0)=y0是初值条件,而f(x,y)是给定的二元函数。
所谓其数值解法,就是求y(x)在离散结点xn处的函数近似值yn 的方法 ,yn≈y(xn)。这些近似值称为常微分方程初值问题的数值解。相邻两个结点之间的距离称为步长。
这里主要介绍两种:单步法和多步法
单步法:在计算y(n+1)时只用到前一步的y(n),因此在有了初值之后就可以逐步往下计算,其代表是龙格- - 库塔( Runge- - Kutta ) 法
多步法:在计算y(n+1)时,除了用到前一步的值y(n)之外, , 还要用到y(n-p)( p=1,2 , … ,k,k>0)的值, , 即前面的k步。其代表就是亚当斯 (Adams) 法
更多介绍可参见这个链接:
https://wenku.baidu.com/view/cafe161b9a6648d7c1c708a1284ac850ad0204fc.html
二、常微分方程求解函数
MATLAB 提供了多个求常微分方程初值问题数值解的函数,一般调用格式为:
[t,y]=solver(filename,tspan,y0,option)
其中,t和y分别给出时间向量和相应的数值解。solver为求常微分方程数值解的函数。filename 是定义 f(t ,y) 的函数名,该函数必须返回一个列向量。
tspan 形式为 [t0,tf],表示求解区间。 y0 是初始状态向量。Option 是可选参数,用于设置求解属性,常用的属性包括相对误差值 RelTol(默认值是10的-3次方)和绝对误差值 AbsTol( ( 默认值是10的-6次方)。
常微分方程数值求解函数的统一命名格式:
odennx
ode是Ordinary Differential Equation 的缩写,是常微分方程的意思。
nn 是数字,代表所用方法的阶数。例如,ode23采用2阶龙格- - 库塔( Runge- - Kutta )算法 ,用3阶公式做误差估计来调节步长,具有低等精度。ode45 采用4阶龙格- - 库塔算 法 ,用5阶公式做误差估计来调节步长,具有中等精度。
x是字母,用于标注方法的专门特征。例如, ode15s 、ode23s 中的“s”代表( Stiff ),表示函数适用于刚性方程。
下表列出了求常微分方程数值解的函数:

求解函数 采用方法 适用场合
ode23 2 阶或3阶 Runge- - Kutta 算法,低精度 非刚性
ode45 4 阶或5阶 Runge- - Kutta 算法,中精度 非刚性
ode113 Adams 算法,精度可到10的-3次方至10的-6次方 非刚性,计算时间比 ode45
ode23t 梯形算法 适度刚性
ode15s Gear’s 反向数值微分算法,中精度 刚性
ode23s 2阶 Rosebrock 算法,低精度 刚性,当精度较低时,计算时间比 ode15s
ode23tb 梯形算法,低精度 刚性,当精度较低时,计算时间比 ode15s

先看一个简单例子,y‘=y+3x/x^2,初值y(0)=-2,求解区间为[1,4]。

>> odefun=@(x,y) (y+3*x)/x^2;
tspan=[1 4];
y0=-2;
[x y]=ode45(odefun,tspan,y0)
plot(x,y)

二、刚性问题
有一类常微分方程,其解的分量有的变化很快,有的变化很慢,且相差悬殊,这就是所谓的刚性问题 (Stiff) 。
对于刚性问题,数值解算法必须取很小步长才能获得满意的结果,导致计算量会大大增加。解决刚性问题需要有专门方法。非刚性算法可以求解刚性问题,只不过需要很长的计算时间。
例、假如点燃一个火柴,火焰球迅速增大直至某个临界体积,然后,维持这一体积不变,原因是火焰球内部燃烧耗费的氧气和从球表面所获氧气达到平衡。其简化模型如下:
y’=y2-y3,y(0)=L,0<=t<=2/L
其中, y(t) 代表火焰球半径,初始半径是λ ,它很小。分析 λ 的大小对方程求解过程的影响。

>>  L=0.01;
f=@(t,y) y^2-y^3;
tic;[ t,y ]=ode45(f,[0,2/L],L); toc
disp (['ode45 计算的点数' num2str(length(t))]);
时间已过 0.003704 秒。
ode45 计算的点数157
>>  L=1e-5;
f=@(t,y) y^2-y^3;
tic;[ t,y ]=ode45(f,[0,2/L],L); toc
disp (['ode45 计算的点数' num2str(length(t))]);
时间已过 1.010016 秒。
ode45 计算的点数120565
>> L=1e-5;
f=@(t,y) y^2-y^3;
tic;[ t,y ]=ode15s(f,[0,2/L],L); toc
disp (['ode45 计算的点数' num2str(length(t))]);
时间已过 0.189977 秒。
ode45 计算的点数98

注:tic 和 toc 函数 用来记录 微分方程求解 命令执行的时间 ,使用 tic 函数启动计时器,使用 toc 函数显示从计时器启动到当前所经历的时间。
上述采用了三种不同方法,可以发现,第一种的运行结果表明这时常微分方程不算很刚性;第二种这时计算时间明显加长,计算的点数剧增,表明这时常微分方程表现为刚性;因此对于刚性问题,我们需要改变求解算法,我们选择以“s”结尾的函数,例如第三种方法他们专门用于求解刚性方程。计算时间大大缩短,计算的点数大大减少。
常微分方程数值解法的研究已发展得相当成熟,理论上也颇为完善,小编由于能力有限,只能了解个大概,本文讲解也就写的是比较基础的一些方面,如果大家有更多需求,可以自己查阅相关资料。

关于MATLAB的学习:

大家可以关注我们的知乎专栏——数据可视化和数据分析中matlab的使用:
https://zhuanlan.zhihu.com/c_1131568134137692160

欢迎大家加入我们的MATLAB学习交流群:
953314432
扫码关注我们
发现更多精彩

matlab常微分方程数值求解相关推荐

  1. [MATLAB]常微分方程数值求解(ode45 ode15s ode23 solver)

    本次讨论取材于中南大学<科学计算与matlab语言>内容为常微分方程数值解法 常微分方程数值求解的一般概念 常微分方程数值求解函数 刚性问题 常微分方程数值求解的一般概念 求解常微分方程初 ...

  2. 常微分方程数值解matlab欧拉,matlab 常微分方程数值解法 源程序代码

    matlab 常微分方程数值解法 源程序代码 所属分类:其他 开发工具:matlab 文件大小:16KB 下载次数:41 上传日期:2019-02-13 11:03:29 上 传 者:XWLYF 说明 ...

  3. matlab数学实验报告西安交通大学微分方程模型高为16米,数学实验第二次作业——常微分方程数值求解...

    实验4常微分方程数值解 实验目的: 1. 练习数值积分的计算: 2. 掌握用MATLAB软件求微分方程初值问题数值解的方法: 3. 通过实例学习用微分方程模型解决简化的实际问题: 4. 了解欧拉方法和 ...

  4. 常微分方程数值求解【python】

    简述 这里只考虑最为简单的一种常微分方程 dydx=f(x,y)\frac{dy}{dx} = f(x,y)dxdy​=f(x,y) 然后这里的实例都是以下面这个方程来做展示的. dydx=y∗(1− ...

  5. 微分方程matlab源代码,matlab 常微分方程数值解法 源程序代码

    11.1 Euler方法 380 11.1.1 Euler公式的推导 380 11.1.2 Euler方法的改进 383 11.2 Runge-Kutta方法 385 11.2.1 二阶Runge-K ...

  6. matlab的数值求解实验报告,matlab计算方法实验报告5(数值积分)

    计算方法实验报告(5) 学生姓名杨贤邦学号指导教师吴明芬实验时间2014.4.16地点综合实验大楼203 实验题目数值积分方法 实验目的●利用复化梯形.辛普森公式和龙贝格数值积分公式计算定积分的 近似 ...

  7. matlab之常微分方程(ODE)求解

    问题描述:已知理想全混釜的初始进料条件,并知道各物料的动力学,求解足够长时间后全混釜中各物质的浓度 常微分方程:只包含一个自变量的微分方程是常微分方程(Ordinary differential eq ...

  8. matlab定积分上界求解,定积分问题的数值求解及Matlab实现.pdf

    定积分问题的数值求解及Matlab实现 第28卷第5期 哈 尔滨 商 业 大 学 学报 (自然科学版) Vo1.28No.5 2012年 10月 JournalofHarbinUniversityof ...

  9. matlab求解f非线性微分方程数值解,非线性﹑微分方程数值求解.PPT

    非线性﹑微分方程数值求解 数 学 系 University of Science and Technology of China DEPARTMENT OF MATHEMATICS 计算方法(B) 主 ...

最新文章

  1. Jmeter工具中参数化、正则表达式提取器、响应断言的实现
  2. 网站被降权后该如何操作?
  3. csrf 攻击和防御
  4. 文本编辑器添加文本编辑区
  5. 浅谈如何更好的打开和关闭ADO.NET连接池
  6. 更改span标签样式_CSS 内嵌样式
  7. python opencv2_Python + OpenCV2 系列:2 - 图片操作
  8. 《新一代人工智能发展白皮书(2017年)》重磅发布(100页完整版PPT)
  9. ThihkPHP开发聚合支付系统源码 兼容所有易支付程序
  10. 【英语学习】【Level 07】U06 First Time L3 Subway everyday
  11. 悲催!谷歌员工中位数年薪达 170 万元,却仍买不起房!
  12. Linux线程屏障,线程屏障(基于linuxthreads-2.3)
  13. android手机测试用例,Android手机测试用例-从事手机测试必备
  14. 云效研发效能度量体系,如何展示和解读交付效能数据
  15. Chinese Whisper 人脸聚类算法实现
  16. 输入年份月份实现日历打印,C到C++过渡。
  17. 倾斜摄影房屋轮廓线提取思路
  18. (详细易懂)一篇文章让你读懂到底什么是Ajax
  19. 史蒂夫·乔布斯(简介)
  20. 那些我们对2019技术世界趋势的预测都说准了吗?

热门文章

  1. li、ul如何去除小黑点样式
  2. VB.NET 创建打印机选择列表
  3. 关于UAP无法运行server的解决办法
  4. 微信jsapi获取用户地理位置接口开发(第八课)
  5. MTTR、MTBF、MTTF、MTTD
  6. 解决CentOS 7,ATI显卡,屏幕亮度调节问题
  7. 将安卓系统迁移到鸿蒙OS,华为或很快推出鸿蒙系统手机 将安卓系统迁移到鸿蒙OS只需1-2天即可实现...
  8. html5语音框架,使用Html5多媒体实现微信语音功能
  9. C++入门基础(下)
  10. GMGC数娱节前瞻,好玩好看有逼格