komazsl 发表于 2011-5-22 21:51

ode23s解刚性微分方程问题，谢谢各位帮我看看

本帖最后由 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=;
y0=;
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;
=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();
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值还是初值，根本没有变化，谢谢各位帮我看看，我毕设遇到难题了

meiyongyuandeze 发表于 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

komazsl 发表于 2011-5-23 19:22

本帖最后由 komazsl 于 2011-5-23 19:23 编辑

回复 2 # meiyongyuandeze 的帖子

为什么会产生病态那，我按照一个论文上的公式编写的，他就是用ode23s方法求解的，而且他解出来了，是不是我第二个M文件写的有问题

meiyongyuandeze 发表于 2011-5-23 20:46

回复 3 # komazsl 的帖子

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

komazsl 发表于 2011-5-23 22:31

回复 4 # meiyongyuandeze 的帖子

谢谢我试试，不会再问你，{:4_74:}

kangarooli 发表于 2011-5-25 08:58

我记得ode23s是解含时滞的方程吧！你的方程有时滞项吗
页: [1]
查看完整版本: ode23s解刚性微分方程问题，谢谢各位帮我看看