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

如何解决如何在SageMath中将两个微分方程组进行数值积分？

我正在尝试在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：无法处理这种情况。”

所以我正在寻找微分方程的数值解。但是由于我是SageMath的新手，所以我不怎么进行。有人可以建议我如何进行数值求解吗？

编辑：

与上述等式对应的M（r）-r曲线如下：

解决方法

跟随sage documentation of desolve functions 如果您指定初始条件和积分范围，则类似以下内容的方法应该起作用。 ODE系统是作为导数方程组的线性系统表示的，因此使用鼠尾草的矩阵功能来求解该线性系统，以获得显式一阶系统的符号表达式。

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。