【Mathematica】这个方程哪错了？
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的指数，会发现随着它的增大，整个方程的求解也变得越来越难了，最后终于失败在半途。这个问题的解决方法，前面已经说了，在于初值的精细调节（毕竟我这里的初值是随便给的）。这个工作就我的经验来说是非常艰巨的。对于我这个不了解方程物理背景的人来说那就更困难了。想要比较有效率的求解方程，个人认为最好是从模型层面做一些修改。（比如用别的方法去表示你的非齐次项之类的。）行了，我在这个问题上耗的时间也够多了，就答到这儿。你也可以去我上面所给的那个站（百度知道不允许我直接贴网址，但是我想你应该知道我指的是哪个）碰碰运气。
Mathematica的默认求解精度是10^-16，那么10^-17的数就有可能当做零了。解决的方法是提高求解的精度，有一个选项的，WorkingPrecision -> 30 这样就把求解精度设置为30位了。
