matlab num2str_MATLAB常微分方程数值求解

最新推荐文章于 2023-09-04 01:00:20 发布

weixin_39943868

最新推荐文章于 2023-09-04 01:00:20 发布

阅读量294

点赞数

文章标签： matlab num2str matlab二元一次方程求解 matlab偏微分方程数值解误差 matlab解微分方程 matlab解微分方程组 num2str函数matlab

点击上方蓝字，关注我们吧！！

今天是白露

白露 09 / 07 农历七月廿十

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

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

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

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

其中，其中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