【Mathematica】这个方程哪错了?G = 6.67*10^-11;B = 1.23*10^10;NDSolve[{m[r] m'[r]^(2/3) (2 r m''[r] - m'[r]) == 4 B/(3 G) r^(19/3) ,m[2] == 2,m[1] == 1},m[r],{r,1,2}]为什么会encounter 1/0?

【Mathematica】这个方程哪错了?G = 6.67*10^-11;B = 1.23*10^10;NDSolve[{m[r] m'[r]^(2/3) (2 r m''[r] - m'[r]) == 4 B/(3 G) r^(19/3) ,m[2] == 2,m[1] == 1},m[r],{r,1,2}]为什么会encounter 1/0?
【Mathematica】这个方程哪错了?
G = 6.67*10^-11;
B = 1.23*10^10;
NDSolve[{m[r] m'[r]^(2/3) (2 r m''[r] - m'[r]) == 4 B/(3 G) r^(19/3) ,m[2] == 2,m[1] == 1},m[r],{r,1,2}]
为什么会encounter 1/0?

【Mathematica】这个方程哪错了?G = 6.67*10^-11;B = 1.23*10^10;NDSolve[{m[r] m'[r]^(2/3) (2 r m''[r] - m'[r]) == 4 B/(3 G) r^(19/3) ,m[2] == 2,m[1] == 1},m[r],{r,1,2}]为什么会encounter 1/0?
周六到了,那就来答这个题吧.


首先这里需要声明的是,我没能解出你的方程,但我大致知道你的方程为什么不好解,我接下来的答案将会指出这一点.我相信,只要你运用我将要讲述的内容,结合你对这个方程物理背景的了解,你还是有可能把它修改到可以求解的.


@yang_bigarm （好吧我知道百度知道没有@的功能……）认为这个方程的求解是因为G和B的数值过分极端.这的确是一个原因,G和B的极端数值确实使问题变得更严重了,但是这并非最关键的原因.事实上精度调节并不能解决这里的问题.最简单的验证方法,是把G和B的值都改成1,然后你会发现,你的问题依然存在.


你所遇到的问题,其实是在非线性常微分方程（或方程组）的边值问题求解中非常常见的一类问题,我在百度知道和Stackexchange都见过很多次,当然,我也很缺乏微分方程的理论知识,所以在此并不能对这类问题发生的根本原因做出什么解释,但是,我知道,这类问题的求解,通常只有一个有效的（并且常常需要精确调节的）方法,那就是,打靶法（Shooting Method）.


所谓打靶法,简单地说,就是试探性地给出一组某一个边界上的函数值及导数值,然后以此为起点,尝试寻找可以满足另一侧边界条件的解.上面已经提及,这一方法常常需要精确的调节,具体说来,就是用于试探的某一边界上的导数值的极细微的差别,都有可能导致求解的失败.（这也就是所谓的“对初值极度敏感”,我个人猜想这应该和方程的非线性性质有关,没准还能和混沌什么的扯上关系……算了,毕竟相关知识不足,这里就不妄加揣度了.）


你的这一方程就极为敏感,原因很可能在于你方程右端的非齐次项的变化过于激烈,你在此没有给出你的方程的物理背景,我也只能认出那个万有引用常量,不过,我猜,这应该是源于某种假设吧?那么,这里以一个更弱的假设为例展示一下你的方程应该怎么修改：


G = B = 1;(* G和B的极端取值使这个方程的求解变难了,这里全改成1 *)
sol = NDSolve[{m[r] m'[r]^(2/3) (2 r m''[r] - m'[r]) == 4 B/(3 G) r^(1/3)(* 注意我把指数改小了 *),m[2] == 2,m[1] == 1},m[r],{r,1,2}, Method -> {"Shooting", "StartingInitialConditions" -> {m[2] == 2,m'[2] == 1}}]
Plot[m[r] /.sol,{r,1,2}]



你可以试着修改方程右端的那个r的指数,会发现随着它的增大,整个方程的求解也变得越来越难了,最后终于失败在半途.这个问题的解决方法,前面已经说了,在于初值的精细调节（毕竟我这里的初值是随便给的）.这个工作就我的经验来说是非常艰巨的.对于我这个不了解方程物理背景的人来说那就更困难了.想要比较有效率的求解方程,个人认为最好是从模型层面做一些修改.（比如用别的方法去表示你的非齐次项之类的.）行了,我在这个问题上耗的时间也够多了,就答到这儿.你也可以去我上面所给的那个站（百度知道不允许我直接贴网址,但是我想你应该知道我指的是哪个）碰碰运气.

在求解这个方程的时候，遇到奇异点了。你有没有发现，G的值太小了，B的值又太大了，两个一相除，基本就等于零了。
Mathematica的默认求解精度是10^-16，那么10^-17的数就有可能当做零了。解决的方法是提高求解的精度，有一个选项的，WorkingPrecision -> 30 这样就把求解精度设置为30位了。...

全部展开

在求解这个方程的时候，遇到奇异点了。你有没有发现，G的值太小了，B的值又太大了，两个一相除，基本就等于零了。
Mathematica的默认求解精度是10^-16，那么10^-17的数就有可能当做零了。解决的方法是提高求解的精度，有一个选项的，WorkingPrecision -> 30 这样就把求解精度设置为30位了。

收起