ode解微分方程 出错Unable to meet integration tolerances ……
本帖最后由 南北左右 于 2011-9-20 14:46 编辑http://cache.soso.com/img/img/e100.gif 请高手不吝赐教，纠结好久了http://cache.soso.com/img/img/e100.gif
出错内容如下：
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 CDSh 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=;
=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 发表于 2011-9-15 16:23 static/image/common/back.gif
有点乱，能贴出要求解的方程吗？
你好，方程已经贴出来了，确实比较复杂，能否帮忙看看~~ 南北左右 发表于 2011-9-16 09:19 static/image/common/back.gif
你好，方程已经贴出来了，确实比较复杂，能否帮忙看看~~
乱七八糟的，变量到底是4个还是3个，常数是多少给出来啊， 写个word文档吧， 找起来真泪眼。 TBE_Legend 发表于 2011-9-16 13:36 static/image/common/back.gif
乱七八糟的，变量到底是4个还是3个，常数是多少给出来啊， 写个word文档吧， 找起来真泪眼。 ...
你好，是四个方程四个未知数，u，R，rug，s,分别对应着x（1）x（2）x（3）x（4）；我看了下好像这边不可以上传word，把原帖重新标注了下，m文件为蓝色字，程序为紫色那段，（常数的赋值都在这一段），不知道这样有没有清楚些:D~~ 本帖最后由 TBE_Legend 于 2011-9-17 19:00 编辑
南北左右 发表于 2011-9-16 16:30 static/image/common/back.gif
你好，是四个方程四个未知数，u，R，rug，s,分别对应着x（1）x（2）x（3）x（4）；我看了下好像这边不可 ...
你把word文档压缩成 rar 文件就能上传了。
C无穷大 是多少？？
从结果中看到时间长了 rg 会趋于无穷大， 所以 matlab 会给那样的提示。
mmtc 代码：
g = 98/10;
\L = 850;
mu = 25/100;
h = 10;
d = 12 10^-8;
\ = 3/100;
P0 = 300000;
Tb = 293; T0 = 298;
\g0 = 7/10;
Cs = 134/10000;
Cinf = 0;
re = 2 \L R u/mu;
Cd = 16/re (1 + (8/re + 1/2 (1 + (3315/100)/re^(1/2)))^-1);
Pb = (\g Tb)/(\g0 T0) P0;
Sh = 2/Sqrt (1 - 2/(3 (1 + 9/100 re^(2/3))^(3/4)))^(1/
2) (25/10 Sqrt[(2 R u)/d] + (2 R u)/d);
ode = { u' == 2 g - 3/(4 R) Cd u^2,
R' ==
R/(4 mu) (Pb - P0 - \L g ( h - s ) - (2 \)/R),
\g'[
t] == -((3 \L Sh d (Cs - Cinf))/(2 R^2)) - (
3 \g)/R R',
s' == u,
u == 1/100000,
R == 2/1000,
\g == 7/10,
s == 1/10};
sol = Firstg, s, R}, {t, 0, 0.08}]]
Table[Plot[f /. sol, {t, 0, 0.08}, PlotRange -> All, Frame -> True,
PlotLabel -> f], {f, {u, R, \g, s}}]
本帖最后由 南北左右 于 2011-9-18 19:14 编辑
TBE_Legend 发表于 2011-9-17 18:58 http://forum.simwe.com/static/image/common/back.gif
你把word文档压缩成 rar 文件就能上传了。
C无穷大 是多少？？
真是非常感动，写了那么长的程序代码，太感谢了，你的这个结果比我用matlab的结果要好得多，第一次听说mmtc，看来我得好好学习下~~~~真的十分感激~~ 本帖最后由 南北左右 于 2011-9-20 15:37 编辑
TBE_Legend 发表于 2011-9-17 18:58 http://forum.simwe.com/static/image/common/back.gif
你把word文档压缩成 rar 文件就能上传了。
C无穷大 是多少？？
你好，根据你的程序，我稍微调整了下参数
g = 98/10;
\L = 800;
mu = 25/100;
h = 20;
d = 12 10^-10;
\ = 3/1000;
P0 = 300000;
Tb = 293;
T0 = 298;
\g0 = 668/1000;
Cs = 296/1000;
Cinf = 125/10000;
re = 2 \L R u/mu;
Cd = 16/re (1 + (8/re + 1/2 (1 + (3315/1000)/re^(1/2)))^-1);
Pb = (\g Tb)/(\g0 T0) 101300;
Sh = 2/Sqrt (1 - 2/(3 (1 + 9/100 re^(2/3))^(3/4)))^(1/2) (25/10 +
Sqrt[(2 R u)/d]);
ode = {u' == 2 g - 3/(4 R) Cd u^2,
R' ==
R/(4 mu) (Pb -
P0 - \L g (h - s) - (2 \)/R), \g'[
t] == -((3 \L Sh d (Cs - Cinf))/(2 R^2)) - (3 \g[
t])/R R', s' == u, u == 2/1000,
R == 2/1000, \g == 32/10, s == 1/10};
sol = Firstg, s, R}, {t, 0, 1}]]
Table[Plot[f /. sol, {t, 0, 1}, PlotRange -> All, Frame -> True,
PlotLabel -> f], {f, {u, R, \g, s}}]
出现了这样的错误NDSolve::ndsz: At t == 0.5757468138991848`, step size is effectively zero; singularity or stiff system suspected. >>
图片趋势是预想中的，但是如何能够画出完整的图呢，是不是由于参数设置的问题 导致算到t == 0.5757468138991848就进行不下去了，改了一天了 没啥进展，实在没办法，只好还向你请教下，实在不好意思~~·
还是或者因为我的方程是刚性的，这种情况该怎么办呢~
matlab 和 mathematica 没啥区别， 你不会看结果图吗？ rhog 都成垂直的了， 你让 mathematica 或 matlab 咋算啊？
调参数没用， 关键是你方程的合理性。rhog 应该是啥样的呢？ 密度都成负值了？？
matlab 和 mathematica 没啥区别， 你不会看结果图吗？ rhog 都成垂直的了， 你让 mathematica 或 matlab 咋算啊？
调参数没用， 关键是你方程的合理性。rhog 应该是啥样的呢？ 密度都成负值了？？
本帖最后由 南北左右 于 2011-9-20 19:15 编辑
TBE_Legend 发表于 2011-9-20 18:00 http://forum.simwe.com/static/image/common/back.gif
matlab 和 mathematica 没啥区别， 你不会看结果图吗？ rhog 都成垂直的了， 你让 mathematica 或 matlab...
多谢帮助~~~经过修改，结果好多了，另外，还想请教个问题，之前那个图中rhog 出现了y=2.8这条横线，是什么原因？
页:
[1]