【求求】DAE求解三元两次微分方程组如何避开0处的奇点?
mathematica吧
全部回复
仅看楼主
level 1
CheyneZ 楼主
Mathematica 12.0.0.0,Mac版
一个三元微分方程组用DAE求解,Vt、v 和 Lambda 是待求变量 [开心]
三个初值条件 Vt[0]==0,v[0]==383.431,\[Lambda][0]==0.0291924 确定是没有错的。
但让他求解的时候总是报错[怒],说找不到能满足残差函数的初始条件。
我怀疑是 t=0 处为 Vt 的一阶倒数的奇点(因此手算的时候也求不出 Vt ’ [0] 等于多少 ),想问问怎么让它在计算的时候避开这个奇点呢?
同时,这个方程还有2个终值条件 Vt [Inf] == 1.8 * 10^(-7),v [Inf] == 0;
但软件说 NDSolve只能求解初值问题??
调试的时候参考了官方教学文档,然鹅还是没有解决:
https://reference.wolfram.com/language/tutorial/NDSolveDAE.html
最后想问的是求解这种方程的Method用什么好?之前试过 "Collocation","Shooting","QR",都没有解决问题(我其实还不明白这几种method的原理 TaT [委屈]
==============源代码如下==================
Clear[Vt,v,\[Lambda],t,dae];
Clear[Derivative];
Clear["Global`*"];
d= 250*10^(-6);
Pm = 2*10^5;
M = 0.029;
\[Rho] = 2.33;
R = 8.31;
L = 1*10^(-2);
T = 293;
\[Nu] = 7.73 * 10^(-6);
DAE={
((\[Rho]*R*T/(M*Pm))-1)*v[t]*(Vt'[t])^2+v[t]*Vt''[t]*Vt[t]==v'[t]*Vt'[t]*Vt[t],
Vt'[t]==((\[Pi]*d^2)/4)*(Pm*v[t]/(Pm-0.5*\[Lambda][t]*\[Rho]*L*(v[t]^2)/d)),(\[Lambda][t])^(-0.5)==Log[10,(d^2*(v[t])^2*\[Lambda][t])/(2.51^2*\[Nu]^2)]
};
a=NDSolve[{DAE,Vt[0]==0,v[0]==383.431,\[Lambda][0]==0.0291924},{Vt,v,\[Lambda]},{t,0,10}]
2019年11月28日 10点11分 1
吧务
level 15
“同时,这个方程还有2个终值条件 Vt [Inf] == 1.8 * 10^(-7),v [Inf] == 0;”如果这两个条件也没错,那你的方程组很可能错了。
要让Mathematica解出这个并不难,基本思想是把方程形式尽量向接近常微分方程组标准形式的形式转化。(具体原理参Stackexchange问题《What's behind Method -> {“EquationSimplification” -> “Residual”}》,问题编号136114,具体链接就不贴了,怕吞。)
首先通过观察可以发现,利用第3个方程可以消掉v[t]:
rule = Solve[{(\[Lambda][t])^(-1/2) ==
Log[10, (d^2*(v[t])^2*\[Lambda][t])/((251/100)^2*\[Nu]^2)]}, v[t]][[1]]
newsys = Most@dae /. rule /. D[rule, t];
到这一步NDSolve依旧有困难,那么我们进一步消掉Vt的导数:
finalsys = {newsys[[2]], newsys[[1]] /. Rule @@ newsys[[2]] /. D[Rule @@ newsys[[2]], t]};
至此,NDSolve依旧解不了这个问题,那么我们适当地错开初值条件:
test = With[{eps = 0},
NDSolveValue[{finalsys,
Vt[eps] == -10^-3, \[Lambda][eps] == 0.0291924}, {Vt, \[Lambda]}, {t, eps, 10}]]
所得图像为
ListPlot[test, PlotRange -> All]
显然Vt在无穷远处不符合你说的条件。
把-10^-3改成10^-3所得的结果类似,因此有理由认为方程本身错了。
2019年12月07日 12点12分 2
1