0

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

Richard
  • 155
  • 6
  • The image is very, very unclear tho – Shivam Jha Sep 07 '20 at 14:40
  • 1
    If you rearrange the equations to get `dM/dr = ...` and `da/dr = ...`, you should be able to use `desolve_system_rk4` https://doc.sagemath.org/html/en/reference/calculus/sage/calculus/desolvers.html#sage.calculus.desolvers.desolve_system_rk4 or `desolve_odeint` https://doc.sagemath.org/html/en/reference/calculus/sage/calculus/desolvers.html#sage.calculus.desolvers.desolve_odeint – John Palmieri Sep 07 '20 at 16:54

1 Answers1

1

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.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
  • I had used the code you had written in the answer with r=times with the range of r being (3,1000). The initial condition was chosen to be the point of intersection of the two curves in the figure in the question (I had edited the question), but I am getting an error. One thing that I need to mention is that the point of intersection of the two curves is a critical point where dM/dr is of 0/0 form and one needs to evaluate dM/dr at that point using the L'Hospital rule. – Richard Sep 07 '20 at 20:20