matlab num2str_MATLAB常微分方程数值求解

aea2d43f1e305b8273997fde3c658de4.png 点击上方蓝字,关注我们吧!! d8bf099a255c4c5cd0a0ba7e1a47d0bd.gif

今天是白露138071e97b2d0a88c1f8446d94b92f67.png138071e97b2d0a88c1f8446d94b92f67.png

白露 09 / 07 农历七月廿十 9c7e35cf1d161417ed412cc620dcfd14.png

本节将介绍常微分方程初值问题的数值求解,主要内容分为三个部分:常微分方程数值求解的概念、求解函数及刚性问题

一、常微分方程数值求解的一般概念

首先,凡含有参数,未知函数和未知函数导数 (或微分) 的方程,称为微分方程,有时简称为方程,未知函数是一元函数的微分方程称作常微分方程,未知函数是多元函数的微分方程称作偏微分方程。微分方程中出现的未知函数最高阶导数的阶数,称为微分方程的阶。

常微分方程数值求解,指研究求解初值问题各类数值方法的构造、理论分析与数值实现问题。研究的主要对象为一阶方程组初值问题:

3b6cf263db0b7d4400d1ec2ab291e952.png

其中,其中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 ),表示函数适用于刚性方程。

下表列出了求常微分方程数值解的函数:

求解函数采用方法适用场合
ode232 阶或3阶 Runge- - Kutta 算法,低精度非刚性
ode454 阶或5阶 Runge- - Kutta 算法,中精度非刚性
ode113Adams 算法,精度可到10的-3次方至10的-6次方非刚性,计算时间比 ode45
ode23t梯形算法适度刚性
ode15sGear’s 反向数值微分算法,中精度刚性
ode23s2阶 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)

31ce160670054cf815f22dfb2d6a95ec.png

三、刚性问题

有一类常微分方程,其解的分量有的变化很快,有的变化很慢,且相差悬殊,这就是所谓的刚性问题 (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

389efdd523a8abd8d97d162ca8273d94.png 389efdd523a8abd8d97d162ca8273d94.png a15b24829f0cb772f5b53fcde3e26c2c.png扫码关注我们吧!出品:Asoul水云天课堂工作室 550acf46c30e5dce8d9a8b36454a48b0.png在看点一下
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值