ode解微分方程 出错Unable to meet integration tolerances ……

南北左右

请高手不吝赐教，纠结好久了
出错内容如下：
Warning: Failure at t=1.426822e+001.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed
(2.842171e-014) at time t.
> In ode15s at 753

m文件：
function dx=equa(t,x)
global g ruL miu h sga P0 Tb T0 P rug0 k C0 Re CD  Sh Pe CS ;
Re=2*ruL*x(2)*x(1)/miu;
%CD=16*(1+1/(8/Re+(1+3.315/(Re^0.5))/2))/Re;
CD=8.*x(2)*g/3./x(1).^2;%这里的CD是上面那个式子的简化，为了计算方便，这个方程本身也是对的
Pe=2*x(2)*x(1)/k;
%Sh=2/sqrt(pi)*(1-2/3/(1+0.09*Re^(2/3))^(3/4))^(0.5)*(2.5+sqrt(Pe));
Sh=0.65*sqrt(Pe);%这里的Sh是上面那个式子的简化，为了计算方便，这个方程本身也是对的

dx=zeros(4,1);
dx(1)=2*g-0.75.*CD.*x(1).^2./x(2);
dx(2)=x(2)./4./miu.*(x(3).*Tb.*P./rug0./T0-P0-ruL.*g.*(h-x(4))-2.*sga./x(2));
dx(3)=-3.*ruL.*Sh.*k.*(CS-C0)/2./x(2)^2-3.*x(3)/4/miu.*(x(3).*Tb.*P./rug0./T0-P0-ruL*g*(h-x(4))-2*sga./x(2));
dx(4)=x(1);


程序：%ode15s

global g ruL miu h sga P0 Tb T0 P rug0 k CS C0;
%Methane
g=9.8;
ruL=850;  %density
miu=0.25; %kinematic viscosity coefficient
h=10;     %depth of oil
sga=0.03; %surface tension
P0=300000; %pressure above the oil surface
Tb=293;    %bubble temperature
T0=298;    %surrounding temperature
P=101300;  %1atm
rug0=0.7;  %t=25?,1atm
k=1.2*10^(-7);  %D
C0=0;           %concentration in the surrounding
CS=0.0134;      %concentration in the bubble

%x(1) velocity
%x(2)bubble size
%x(3)bubble dencity  
%x(4)rising distance s

t0=0:100;
x0=[0.00001;0.002;0.7;0.1];
[t,x]=ode15s(@equa,t0,x0);
figure;
plot(t,x(:,1));
figure;
plot(t,x(:,2));
figure;
plot(t,x(:,3));
figure;
plot(t,x(:,4));


要求解的方程：

TBE_Legend

有点乱，能贴出要求解的方程吗？

南北左右

TBE_Legend 发表于 2011-9-15 16:23
有点乱，能贴出要求解的方程吗？

你好，方程已经贴出来了，确实比较复杂，能否帮忙看看~~

TBE_Legend

南北左右 发表于 2011-9-16 09:19
你好，方程已经贴出来了，确实比较复杂，能否帮忙看看~~

乱七八糟的，变量到底是4个还是3个，常数是多少给出来啊， 写个word文档吧， 找起来真泪眼。

南北左右

TBE_Legend 发表于 2011-9-16 13:36
乱七八糟的，变量到底是4个还是3个，常数是多少给出来啊， 写个word文档吧， 找起来真泪眼。 ...

你好，是四个方程四个未知数，u，R，rug，s,分别对应着x（1）x（2）x（3）x（4）；我看了下好像这边不可以上传word，把原帖重新标注了下，m文件为蓝色字，程序为紫色那段，（常数的赋值都在这一段），不知道这样有没有清楚些:D~~

TBE_Legend

南北左右 发表于 2011-9-16 16:30
你好，是四个方程四个未知数，u，R，rug，s,分别对应着x（1）x（2）x（3）x（4）；我看了下好像这边不可 ...

你把word文档压缩成 rar 文件就能上传了。

C无穷大 是多少？？
从结果中看到时间长了 rg 会趋于无穷大， 所以 matlab 会给那样的提示。

mmtc 代码：

1. g = 98/10;
   
2. \[Rho]L = 850;
   
3. mu = 25/100;
   
4. h = 10;
   
5. d = 12 10^-8;
   
6. \[Sigma] = 3/100;
   
7. P0 = 300000;
   
8. Tb = 293; T0 = 298;
   
9. \[Rho]g0 = 7/10;
   
10. Cs = 134/10000;
    
11. Cinf = 0;
    
12. 
    
13. re = 2 \[Rho]L R[t] u[t]/mu;
    
14. Cd = 16/re (1 + (8/re + 1/2 (1 + (3315/100)/re^(1/2)))^-1);
    
15. Pb = (\[Rho]g[t] Tb)/(\[Rho]g0 T0) P0;
    
16. Sh = 2/Sqrt[Pi] (1 - 2/(3 (1 + 9/100 re^(2/3))^(3/4)))^(1/
    
17.    2) (25/10 Sqrt[(2 R[t] u[t])/d] + (2 R[t] u[t])/d);
    
18. ode = { u'[t] == 2 g - 3/(4 R[t]) Cd u[t]^2,
    
19.    R'[t] ==
    
20.     R[t]/(4 mu) (Pb - P0 - \[Rho]L g ( h - s[t] ) - (2 \[Sigma])/R[t]),
    
21.    \[Rho]g'[
    
22.      t] == -((3 \[Rho]L Sh d (Cs - Cinf))/(2 R[t]^2)) - (
    
23.       3 \[Rho]g[t])/R[t] R'[t],
    
24.    s'[t] == u[t],
    
25.    u[0] == 1/100000,
    
26.    R[0] == 2/1000,
    
27.    \[Rho]g[0] == 7/10,
    
28.    s[0] == 1/10};
    
29. sol = First[NDSolve[ode, {u, \[Rho]g, s, R}, {t, 0, 0.08}]]
    
30. Table[Plot[f /. sol, {t, 0, 0.08}, PlotRange -> All, Frame -> True,
    
31.   PlotLabel -> f], {f, {u[t], R[t], \[Rho]g[t], s[t]}}]

复制代码

南北左右

TBE_Legend 发表于 2011-9-17 18:58
你把word文档压缩成 rar 文件就能上传了。

C无穷大 是多少？？

真是非常感动，写了那么长的程序代码，太感谢了，你的这个结果比我用matlab的结果要好得多，第一次听说mmtc，看来我得好好学习下~~~~真的十分感激~~

南北左右

TBE_Legend 发表于 2011-9-17 18:58
你把word文档压缩成 rar 文件就能上传了。

C无穷大 是多少？？

你好，根据你的程序，我稍微调整了下参数
g = 98/10;
\[Rho]L = 800;
mu = 25/100;  
h = 20;
d = 12 10^-10;  
\[Sigma] = 3/1000;
P0 = 300000;
Tb = 293;            
T0 = 298;
\[Rho]g0 = 668/1000;
Cs = 296/1000;
Cinf = 125/10000;   

re = 2 \[Rho]L R[t] u[t]/mu;
Cd = 16/re (1 + (8/re + 1/2 (1 + (3315/1000)/re^(1/2)))^-1);
Pb = (\[Rho]g[t] Tb)/(\[Rho]g0 T0) 101300;
Sh = 2/Sqrt[Pi] (1 - 2/(3 (1 + 9/100 re^(2/3))^(3/4)))^(1/2) (25/10 +
     Sqrt[(2 R[t] u[t])/d]);
ode = {u'[t] == 2 g - 3/(4 R[t]) Cd u[t]^2,
   R'[t] ==
    R[t]/(4 mu) (Pb -
       P0 - \[Rho]L g (h - s[t]) - (2 \[Sigma])/R[t]), \[Rho]g'[
     t] == -((3 \[Rho]L Sh d (Cs - Cinf))/(2 R[t]^2)) - (3 \[Rho]g[
          t])/R[t] R'[t], s'[t] == u[t], u[0] == 2/1000,
   R[0] == 2/1000, \[Rho]g[0] == 32/10, s[0] == 1/10};
sol = First[NDSolve[ode, {u, \[Rho]g, s, R}, {t, 0, 1}]]

Table[Plot[f /. sol, {t, 0, 1}, PlotRange -> All, Frame -> True,
  PlotLabel -> f], {f, {u[t], R[t], \[Rho]g[t], s[t]}}]
出现了这样的错误NDSolve::ndsz: At t == 0.5757468138991848`, step size is effectively zero; singularity or stiff system suspected. >>
图片趋势是预想中的，但是如何能够画出完整的图呢，是不是由于参数设置的问题 导致算到t == 0.5757468138991848就进行不下去了，改了一天了 没啥进展，实在没办法，只好还向你请教下，实在不好意思~~·
还是或者因为我的方程是刚性的，这种情况该怎么办呢~

TBE_Legend

matlab 和 mathematica 没啥区别， 你不会看结果图吗？ rhog 都成垂直的了， 你让 mathematica 或 matlab 咋算啊？

调参数没用， 关键是你方程的合理性。rhog 应该是啥样的呢？ 密度都成负值了？？

TBE_Legend

matlab 和 mathematica 没啥区别， 你不会看结果图吗？ rhog 都成垂直的了， 你让 mathematica 或 matlab 咋算啊？

调参数没用， 关键是你方程的合理性。rhog 应该是啥样的呢？ 密度都成负值了？？

南北左右

TBE_Legend 发表于 2011-9-20 18:00
matlab 和 mathematica 没啥区别， 你不会看结果图吗？ rhog 都成垂直的了， 你让 mathematica 或 matlab  ...

  多谢帮助~~~经过修改，结果好多了，另外，还想请教个问题，之前那个图中rhog 出现了y=2.8这条横线，是什么原因？