Ticket #8950: trac_8950_desolve_numerical.patch

sage/calculus/all.py

@@ -10 +10 @@
 
 from desolvers import (desolve, desolve_laplace, desolve_system,
                        eulers_method, eulers_method_2x2, 
-                       eulers_method_2x2_plot, desolve_rk4, desolve_system_rk4)
+                       eulers_method_2x2_plot, desolve_rk4, desolve_system_rk4, desolve_numerical)
 
 from var import (var, function, clear_vars)
 

sage/calculus/desolvers.py

@@ -1053 +1053 @@
     sol.extend(sol_2)
 
     return sol
+
+def desolve_numerical(des, ics, times, dvars, ivar=None, ComputeJacobian=False, args=(), ml=None, mu=None, rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0):
+        """
+        Solves numerically a system of first-order ordinary differential equations using odeint from scipy.integrate module.
+
+        INPUT:
+                des  -- right hand sides of the system
+                ics  -- initial conditions 
+                times -- a sequence of time points in which the solution must be found
+                dvars -- dependent variables. ATTENTION: the order must be the same as in des, that means: d(dvars[i])/dt=des[i]
+                ivar -- independent variable, default is t
+                ComputeJacobian -- boolean. If True, the Jacobian of des is computed and used during the integration of Stiff Systems. Default value is False.
+
+                Other Parameters (taken from the documentation of odeint function from scipy.integrate module)
+                ----------------------------------------------------------------------------------------------
+                ml, mu : integer
+                    If either of these are not-None or non-negative, then the
+                    Jacobian is assumed to be banded.  These give the number of
+                    lower and upper non-zero diagonals in this banded matrix.
+                    For the banded case, Dfun should return a matrix whose
+                    columns contain the non-zero bands (starting with the
+                    lowest diagonal).  Thus, the return matrix from Dfun should
+                    have shape len(y0) * (ml + mu + 1) when ml >=0 or mu >=0
+                rtol, atol : float
+                    The input parameters rtol and atol determine the error
+                    control performed by the solver.  The solver will control the
+                    vector, e, of estimated local errors in y, according to an
+                    inequality of the form::
+                        max-norm of (e / ewt) <= 1
+                    where ewt is a vector of positive error weights computed as::
+                        ewt = rtol * abs(y) + atol
+                    rtol and atol can be either vectors the same length as y or scalars.
+                tcrit : array
+                    Vector of critical points (e.g. singularities) where integration
+                    care should be taken.
+                h0 : float, (0: solver-determined)
+                    The step size to be attempted on the first step.
+                hmax : float, (0: solver-determined)
+                    The maximum absolute step size allowed.
+                hmin : float, (0: solver-determined)
+                    The minimum absolute step size allowed.
+                ixpr : boolean
+                    Whether to generate extra printing at method switches.
+                mxstep : integer, (0: solver-determined)
+                    Maximum number of (internally defined) steps allowed for each
+                    integration point in t.
+                mxhnil : integer, (0: solver-determined)
+                    Maximum number of messages printed.
+                mxordn : integer, (0: solver-determined)
+                    Maximum order to be allowed for the nonstiff (Adams) method.
+                mxords : integer, (0: solver-determined)
+                    Maximum order to be allowed for the stiff (BDF) method.
+
+
+        OUTPUT:
+                Returns a list with the solution of the system at each time in times.
+
+
+        EXAMPLES:
+        1) Lotka Volterra Equations:
+
+                sage: from sage.calculus.desolvers import desolve_numerical
+                sage: x,y=var('x,y')
+                sage: sol=desolve_numerical([x*(1-y),-y*(1-x)],[0.5,2],srange(0,10,0.1),[x,y])
+                sage: p=line(zip(sol[:,0],sol[:,1]))
+                sage: p.show()
+
+
+        2) Lorenz Equations:
+
+                sage: from sage.calculus.desolvers import desolve_numerical
+                sage: from sage.plot.plot3d.shapes2 import Line
+                sage: x,y,z=var('x,y,z')
+                sage: # Next we define the parameters           
+                sage: sigma=10
+                sage: rho=28
+                sage: beta=8/3
+                sage: # The Lorenz equations
+                sage: lorenz=[sigma*(y-x),x*(rho-z)-y,x*y-beta*z] 
+                sage: # Time and initial conditions
+                sage: times=srange(0,50.05,0.05)
+                sage: ics=[0,1,1]
+                sage: sol=desolve_numerical(lorenz,ics,times,[x,y,z],rtol=1e-13,atol=1e-14)  
+
+
+        3) One-dimensional Stiff system:
+
+                sage: from sage.calculus.desolvers import desolve_numerical
+                sage: y= var('y')
+                sage: epsilon=0.01
+                sage: f=y^2*(1-y)
+                sage: ic=epsilon
+                sage: times=srange(0,2/epsilon,1)
+                sage: sol=desolve_numerical(f,ic,times,y,rtol=1e-13,atol=1e-14,ComputeJacobian=True)
+                sage: p=points(zip(times,sol))
+                sage: p.show()
+
+
+        4) Another Stiff system with some optional parameters with no default value:
+
+                sage: y1,y2,y3=var('y1,y2,y3')
+                sage: f1=77.27*(y2+y1*(1-8.375*1e-6*y1-y2))
+                sage: f2=1/77.27*(y3-(1+y1)*y2)
+                sage: f3=0.16*(y1-y3)
+                sage: f=[f1,f2,f3]
+                sage: ci=[0.2,0.4,0.7]
+                sage: times=srange(0,10,0.01)
+                sage: sol=desolve_numerical(f,ci,times,[y1,y2,y3],ComputeJacobian=True,rtol=1e-3,atol=1e-4,h0=0.1,hmax=1,hmin=1e-4,mxstep=1000,mxords=17)
+
+
+        AUTHOR: Oriol Castejon (05-2010)
+        """
+
+        from scipy.integrate import odeint
+        from sage.ext.fast_eval import fast_float
+        from sage.calculus.functions import jacobian
+
+        if ivar==None:
+                from sage.symbolic.ring import var
+                ivar=var('t')
+
+        # one-dimensional systems:
+        if len(dvars)==1 or len(dvars)==0:
+                func=fast_float(des,dvars,ivar)
+                if ComputeJacobian==False:
+                        Dfun=None
+                else:
+                        J=jacobian(des,dvars)
+                        J=J[0][0]
+                        J=fast_float(J,dvars,ivar)
+                        Dfun=lambda y,t: [J(y,t)]
+
+