pythonnumerical-methodsdifferential-equationssagenumerical-analysis

How to integrate a set of two differential equations numerically in SageMath?

I am trying to solve the following two differential equation (numerically) in SageMath:

1

2

My goal is to obtain the plot of M(r)-r.

I tried the following code:

    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])

But this is returning an error that says: "TypeError: ECL says: Error executing code in Maxima: desolve: can't handle this case."

So I am looking for a numerical solution of the differential equations. But since I am new to SageMath, I don't how to proceed. Can someone suggest me how to proceed for obtaining a numerical solution?

EDIT:

The M(r)-r plot corresponding to the above equations is the following:

3

Solution

  • Following the sage documentation of desolve functions something like the following should work, if you specify the initial conditions and the integration range. The ODE system is presented as a linear system of equations in the derivatives, so use the matrix functionality of sage to solve this linear system to get the symbolic expressions for the explicit first order system.

    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)

    The result is a list of state vectors at the times in times. So with r=times and M=sol[:,0] you should be able to plot M-r against r.