声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2688|回复: 5

[编程技巧] ode23s解刚性微分方程问题,谢谢各位帮我看看

[复制链接]
发表于 2011-5-22 21:51 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
本帖最后由 komazsl 于 2011-5-22 22:01 编辑

第一个M文件
function dy=rigid(t,y,flag,Ce,Ta,Tb,Tc,Td,Te,Tf,Tg,Th,Ti,Tj,Tk,Tl,Tm,Tn,To,Tp,Tq,Tr);         
dy=zeros(9,1);
dy(1)=-Tb*y(1)*y(4)-Tc*y(3)*y(1)+Td*y(3)*y(4)+Tg*y(3)*y(5)-Th*y(1)*y(4)-Ti*y(1)*y(6)-Tj*y(2)*y(2)*y(5)+2*Tk*y(3)*y(7)+Tl*y(7)*y(4)-Tp*y(8)*y(1);
dy(2)=Tc*y(3)*y(1);
dy(3)=Tb*y(1)*y(4)-Tc*y(3)*y(1)-Td*y(3)*y(4)-Tg*y(3)*y(5)-Tk*y(3)*y(7);
dy(4)=2*Ta*y(5)*Ce-Tb*y(1)*y(4)+Tc*y(3)*y(1)-Td*y(3)*y(4)-Te*y(4)*y(5)-Tf*y(4)*y(6)+Tg*y(3)*y(5)-Th*y(1)*y(4)-Tl*y(7)*y(4)-Tm*y(7)*y(4)-To*y(8)*y(4);
dy(5)=-Ta*y(5)*Ce+Tb*y(1)*y(4)-Te*y(4)*y(5)+2*Tf*y(4)*y(6)-Tg*y(3)*y(5)+Ti*y(1)*y(6)-Tj*y(2)*y(2)*y(5)+Tl*y(7)*y(4)+Tn*y(7)*y(6)+To*y(8)*y(4);
dy(6)=Te*y(4)*y(5)-Tf*y(4)*y(6)-Ti*y(1)*y(6)-Tn*y(7)*y(6);
dy(7)=Th*y(1)*y(4)+Ti*y(1)*y(6)-Tk*y(3)*y(7)-Tl*y(7)*y(4)-Tm*y(7)*y(4)-Tn*y(7)*y(6)+To*y(8)*y(4)+2*Tp*y(8)*y(1)-Tq*y(8)*y(7)+Tr*y(9);
dy(8)=Tm*y(7)*y(4)+Tn*y(7)*y(6)-To*y(8)*y(4)-Tp*y(8)*y(1)-Tq*y(8)*y(7)+Tr*y(9);
dy(9)=Tq*y(8)*y(7)-Tr*y(9);
end
第二个
tspan=[0,6];
y0=[0.001,0.999,0,0,0,0,0,0,0];
Ce=1.00E+13;
Ta=1.29E+15;
Tb=2.77E+9;
Tc=2.04E+13;
Td=1.45E+11;
Te=1.69E+12;
Tf=4.81E+9;
Tg=5.36E+7;
Th=8.43E+11;
Ti=1.09E+10;
Tj=7.09E+9;
Tk=1.81E+12;
Tl=5.84E+12;
Tm=4.38E+7;
Tn=1.97E+9;
To=1.02E+13;
Tp=1.20E+13;
Tq=4.57E+12;
Tr=1.32E+11;
options=odeset('RelTol',1e-6);
t0=cputime;
[t,y]=ode23s('rigid',tspan,y0,options,Ce,Ta,Tb,Tc,Td,Te,Tf,Tg,Th,Ti,Tj,Tk,Tl,Tm,Tn,To,Tp,Tq,Tr);
time=cputime-t0;
figure;
plot(t,y(:,1),'*');
hold on;
plot(t,y(:,2),'o');
set(gca,'Fontsize',12);
ylim([0,0.006]);
xlabel('\itt','Fontsize',16);
set(gcf,'color','w');
算出来出现
Warning: Matrix is close to singular or badly scaled.
         Results may be inaccurate. RCOND = 7.319604e-029.
> In ode23s at 456
  In Untitled11 at 24,
然后能出现图但是y值还是初值,根本没有变化,谢谢各位帮我看看,我毕设遇到难题了
回复
分享到:

使用道具 举报

发表于 2011-5-23 11:49 | 显示全部楼层
从报错来看,可能是你的微分方程病态的原因。
建议你用ode15s来试试!
Matrix is close to singular or badly scaled
  http://www.mathworks.com/matlabc ... /view_thread/170919
  http://groups.google.com/group/c ... ad/f640ae64aa62fbe6
  http://dmpeli.mcmaster.ca/Matlab ... otes/Lecture2_4.htm
from http://forum.vibunion.com/home-spa ... -blog-id-18250.html
 楼主| 发表于 2011-5-23 19:22 | 显示全部楼层
本帖最后由 komazsl 于 2011-5-23 19:23 编辑

回复 2 # meiyongyuandeze 的帖子

为什么会产生病态那,我按照一个论文上的公式编写的,他就是用ode23s方法求解的,而且他解出来了,是不是我第二个M文件写的有问题
发表于 2011-5-23 20:46 | 显示全部楼层
回复 3 # komazsl 的帖子

第二个M文件,我帮你运行了也没发现什么明显的运行错误。
我猜测很可能是某些参数取的不合适,刚性方程本省就是对默写参数敏感依赖的,建议:1. 采用ode15s试试。2.检查参数的取值和初始值取值。

评分

1

查看全部评分

 楼主| 发表于 2011-5-23 22:31 | 显示全部楼层
回复 4 # meiyongyuandeze 的帖子

谢谢我试试,不会再问你,
发表于 2011-5-25 08:58 | 显示全部楼层
我记得ode23s是解含时滞的方程吧!你的方程有时滞项吗
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-29 12:31 , Processed in 0.054557 second(s), 19 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表