# Having trouble in solving two differential equations using desolve_system

I am trying to solve the following two differential equations simultaneously: $$Ma^2\frac{dM}{dr}+(M^2a+6a)\frac{da}{dr}+\frac{1}{r^2}=0$$ $$ar\frac{dM}{dr}+7Mr\frac{da}{dr}+2Ma=0$$ where $M=M(r)$ and $a=a(r)$ are the variables.

I had written the following code in Sage:

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


After writing the above code in Sage, I am getting the following error:

TypeError                                 Traceback (most recent call last)
<ipython-input-26-6bce8159491d> in <module>()
----> 1 desolve_system([de1, de2], [M,a])

/usr/lib/sagemath/local/lib/python2.7/site-packages/sage/calculus/desolvers.py in desolve_system(des, vars, ics, ivar)
822         for dvar, ic in zip(dvars, ics[1:]):
823             dvar.atvalue(ivar==ivar_ic, ic)
--> 824     soln = dvars.parent().desolve(des, dvars)
825     if str(soln).strip() == 'false':
826         raise NotImplementedError("Maxima was unable to solve this system.")

/usr/lib/sagemath/local/lib/python2.7/site-packages/sage/interfaces/interface.py in __call__(self, *args, **kwds)
606
607     def __call__(self, *args, **kwds):
--> 608         return self._parent.function_call(self._name, list(args), kwds)
609
610     def _sage_doc_(self):

/usr/lib/sagemath/local/lib/python2.7/site-packages/sage/interfaces/interface.py in function_call(self, function, args, kwds)
532                                        [s.name() for s in args],
533                                        ['%s=%s'%(key,value.name()) for key, value in kwds.items()])
--> 534         return self.new(s)
535
536     def _function_call_string(self, function, args, kwds):

/usr/lib/sagemath/local/lib/python2.7/site-packages/sage/interfaces/interface.py in new(self, code)
307
308     def new(self, code):
--> 309         return self(code)
310
311     ###################################################################

/usr/lib/sagemath/local/lib/python2.7/site-packages/sage/interfaces/interface.py in __call__(self, x, name)
242
243         if isinstance(x, six.string_types):
--> 244             return cls(self, x, name=name)
245         try:
246             return self._coerce_from_special_method(x)

/usr/lib/sagemath/local/lib/python2.7/site-packages/sage/interfaces/interface.py in __init__(self, parent, value, is_name, name)
670                 self._name = parent._create(value, name=name)
671             except (TypeError, RuntimeError, ValueError) as x:
--> 672                 raise TypeError(x)
673
674     def _latex_(self):

TypeError: ECL says: Error executing code in Maxima: desolve: can't handle this case.


Can someone help me with the problem. I am new to Sage and so I could not interpret the error.

edit retag close merge delete

1

Are you looking for a symbolic or a numerical solution? I don't think Maxima can handle such a complex differential equation symbolically. Is there even a solution?

The solution for differential equations, i.e., the plot of M(r)-r can be found in the following link: plot

It would be fine if I could get a numerical solution. Actually I require to obtain the M(r)-r plot.

Sort by » oldest newest most voted

Your system of equations has no "easy" analytical (closed-form) solution It might be possible to get one, but I won't delve into it here.

Your question lacks precision on the meaning of your variables and initial numerical values of your functions. Ihave made a number of (likely incorrect) assumptions :

Data:

a,M=function("a,M")
r=var("r")
Eq1=M(r)*a(r)^2*diff(M(r),r)+(M(r)^2*a(r)+6*a(r))*diff(a(r),r)+1/r^2==0
Eq2=a(r)*M(r)*diff(M(r),r)+7*M(r)*r*diff(a(r),r)+2*M(r)*a(r)==0


desolve_system_rk4 is a function integrating systems of differential equations. It has a set of peculiar requirements :

• Equations are expressions of the first-order derivatives of the functions of interest :

.

var(("da,dM"))
Sys1=[u.subs({diff(a(r),r):da, diff(M(r),r):dM}) for u in (Eq1,Eq2)]
Sys2=solve(Sys1,[da,dM])
Sys3=[u.subs({da:diff(a(r),r), dM:diff(M(r),r)}) for u in Sys2]

• It uses representation of functions by variables (deprecated in Sage...)

.

var("va,vM")
Sys4=[u.rhs().subs({a(r):va, M(r):vM}) for u in Sys3]


Now we can solve the system, conditional to our assumptions on the initial values :

NSol=desolve_system_rk4(Sys4,[va,vM],ics=[1,6,0],ivar=r,end_points=exp(3).n())


The result is a list of (independent variable, value1, value2,...) "points". We can plot those points (and functions thereof), use them to create interpolation functions, etc...

line([[log(r),M-r] for r,a,M in NSol])


HTH,

more

Maybe desolve_odeint can be used the same way?

Yup. As well as using ode_solver (more Fortranish than Pythonic...).

ISTR that there is a also a boatload of Python libraries floating on the 'Net for DE solving. But this is the kind of problems easy to solve but fiendishly difficult to solve well.

Numerical analysis is a bitch...