如何在SageMath中对一组两个微分方程进行数值积分？

Question

我试图在SageMath中求解以下两个微分方程（数值）：

我的目标是获得M（r）-r的图。

我尝试了以下代码：

    sage: r = var('r')
    sage: M = function('M')(r)
    sage: a = function('a')(r)
    sage: de1 = (M*a*a*diff(M,r) + (M*M*a+6*a)*diff(a,r) + 1/(r*r) == 0)
    sage: de2 = (a*r*diff(M,r) + 7*M*r*diff(a,r) + 2*M*a == 0)
    sage: desolve_system([de1,de2], [M,a])

但这将返回一个错误，该错误表示： “TypeError:ECL表示：执行Maxima:desolve中的代码时出错：无法处理此情况。”

所以我在寻找微分方程的数值解。但是，由于我对数学一无所知，我不知道该如何进行。有人能告诉我如何获得数值解吗

编辑：

与上述方程式对应的M（r）-r图如下所示：

Answer 1

在sage documentation of desolve functions之后 如果您指定了初始条件和积分范围，那么下面这样的方法应该可以工作。ODE系统在导数中表示为线性方程组，因此使用sage的矩阵功能来求解该线性系统，以获得显式一阶系统的符号表达式

A = matrix([[ M*a*a, M*M*a+6*a + 1/(r*r)],
            [ a*r, 7*M*r]])
B= matrix([[ 1/(r*r)], [2*M*a]])

f = -A^-1*B
times=srange(ti,tf,dt)
ics=[M0,a0]
sol=desolve_odeint(f,ics,times,dvars = [M,a], ivar = r)

结果是times中时间的状态向量列表。因此，使用r=times和M=sol[:,0]，您应该能够绘制M-r与r的对比图