注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

一粒浮尘

飘渺虚无

 
 
 

日志

 
 

刚性/非刚性问题(Stiff/Nonstiff)的Matlab解法  

2013-02-24 12:12:47|  分类: Matlab |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
在工程实践中,我们经常遇到一些ODEs,其中某些解变换缓慢,另一些变化很快,且相差悬殊的微分方程,这就是所谓的刚性问题(Stiff),对于所有解的变化相当我们则称为非刚性问题(Nonstiff)。

对于刚性问题一般不适合使用ode45这类函数求解。

由于非刚性问题我们使用的多比较多,我们就不多说,下面主要讲解下Stiff

看一个例子,考虑下面的微分方程

7.gif 

【解】题目中已经帮我们完成了方程组转换,下面我们就直接编程了

(1)编写微分方程函数
  1. odefun=@(t,x)[0.04*(1-x(1))-(1-x(2))*x(1)+0.0001*(1-x(2))^2
  2.     -1e4*x(1)+3000*(1-x(2))^2];
  3. x0=[0 1];
  4. tspan=[0 100];
  5. options=odeset('reltol',1e-6,'abstol',1e-8);
复制代码
(2)对于这个刚性问题,我们先使用ode45函数试试(反正我们是没有信心等待它,太慢了)
  1. tic;[t,y]=ode45(odefun,tspan,x0,options);toc
  2. disp(['ode45计算的点数' num2str(length(t))])
复制代码
(3)换用ode15s函数试试看看
  1. tic;[t2,y2]=ode15s(odefun,tspan,x0,options);toc
  2. disp(['ode15s计算的点数' num2str(length(t2))])
复制代码
(4)绘图看看结果
  1. figure('numbertitle','off','name','Stiff ODEs Demo—by Matlabsky')
  2. subplot(121)
  3. title('Using ode45')
  4. plot(t,y)
  5. subplot(122)
  6. plot(t2,y2)
  7. title('Using ode15s')
复制代码
运行的结果如下,我们可以比较下ode45和ode15s的差距
  1. Elapsed time is 171.005688 seconds.
  2. ode45计算的点数356981
  3. Elapsed time is 0.397287 seconds.
  4. ode15s计算的点数188
复制代码

p2.jpg


  评论这张
 
阅读(3561)| 评论(1)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017