Ticket #8132: trac_8132.patch

File trac_8132.patch, 82.0 KB (added by Robert Marik, 13 years ago)
  • doc/en/constructions/calculus.rst

    # HG changeset patch
    # User Robert Marik <marik@mendelu.cz>
    # Date 1264898007 -3600
    # Node ID ea16778cc7ae89ebf12c73c371e2fdc0cbab046e
    # Parent  b89cf68ccc1ca1cfc9f29e86e9b61770b41d7df7
    Ticket #8132 - fix documentation related to ODE solvers
    
    diff -r b89cf68ccc1c -r ea16778cc7ae doc/en/constructions/calculus.rst
    a b  
    2929    sage: latex(f.diff(x))
    3030    k x^{3} e^{\left(k x\right)} \sin\left(w x\right) + w x^{3} e^{\left(k x\right)} \cos\left(w x\right) + 3 \, x^{2} e^{\left(k x\right)} \sin\left(w x\right)
    3131
    32 If you type ``view(f.diff('x'))`` another window will open up
     32If you type ``view(f.diff(x))`` another window will open up
    3333displaying the compiled output. In the notebook, you can enter
    3434
    3535::
    3636
    37     f = maxima('x^3 * %e^(k*x) * sin(w*x)')
    38     show(f)
    39     show(f.diff('x'))
     37    var('x k w')
     38    f = x^3 * e^(k*x) * sin(w*x)
     39    show(f)       
     40    show(f.diff(x))
    4041
    4142into a cell and press ``shift-enter`` for a similar result. You can
    42 also call Maxima indirectly using the commands
     43also differentiate and integrate using the commands
    4344
    4445::
    4546
     
    263264Ordinary differential equations
    264265===============================
    265266
    266 Symbolically solving ODEs can be done using 's interface with
    267 Maxima. Numerical solution of ODEs can be done using 's interface
     267Symbolically solving ODEs can be done using Sage interface with
     268Maxima. See
     269
     270::
     271
     272    sage:desolvers?
     273
     274for available commands. Numerical solution of ODEs can be done using Sage interface
    268275with Octave (an experimental package), or routines in the GSL (Gnu
    269276Scientific Library).
    270277
    271 You can solve ODE's symbolically in Sage using the Maxima interface
     278An example, how to solve ODE's symbolically in Sage using the Maxima interface
    272279(do not type the ``...``):
    273280
    274281::
    275282
    276     sage: maxima.de_solve('diff(y,x,2) + 3*x = y', ['x','y'], [1,1,1])
    277     y=3*x-2*%e^(x-1)
    278     sage: maxima.de_solve('diff(y,x,2) + 3*x = y', ['x','y'])
    279     y=%k1*%e^x+%k2*%e^-x+3*x
    280     sage: maxima.de_solve('diff(y,x) + 3*x = y', ['x','y'])
    281     y=(%c-3*(-x-1)*%e^-x)*%e^x
    282     sage: maxima.de_solve('diff(y,x) + 3*x = y', ['x','y'],[1,1])
    283     y=-%e^-1*(5*%e^x-3*%e*x-3*%e)
     283    sage: y=function('y',x); desolve(diff(y,x,2) + 3*x == y, dvar = y, ics = [1,1,1])
     284    3*x - 2*e^(x - 1)
     285    sage: desolve(diff(y,x,2) + 3*x == y, dvar = y)
     286    k1*e^x + k2*e^(-x) + 3*x
     287    sage: desolve(diff(y,x) + 3*x == y, dvar = y)
     288    (3*(x + 1)*e^(-x) + c)*e^x
     289    sage: desolve(diff(y,x) + 3*x == y, dvar = y, ics = [1,1]).expand()
     290    3*x - 5*e^(x - 1) + 3
    284291   
    285     sage: maxima.clear('x'); maxima.clear('f')
    286     sage: maxima.de_solve_laplace("diff(f(x),x,2) = 2*diff(f(x),x)-f(x)",\
    287     ...   ["x","f"], [0,1,2])
    288     f(x)=x*%e^x+%e^x
     292    sage: f=function('f',x); desolve_laplace(diff(f,x,2) == 2*diff(f,x)-f, dvar = f, ics = [0,1,2])
     293    x*e^x + e^x
    289294               
    290     sage: maxima.clear('x'); maxima.clear('f')           
    291     sage: f = maxima.de_solve_laplace("diff(f(x),x,2) = 2*diff(f(x),x)-f(x)",\
    292     ...   ["x","f"])
    293     sage: f
    294     f(x)=x*%e^x*('at('diff(f(x),x,1),x=0))-f(0)*x*%e^x+f(0)*%e^x
    295     sage: f.display2d()
    296                                                    !
    297                                        x  d        !                  x          x
    298                             f(x) = x %e  (-- (f(x))!     ) - f(0) x %e  + f(0) %e
    299                                           dx       !
    300                                                    !x = 0
     295    sage: desolve_laplace(diff(f,x,2) == 2*diff(f,x)-f, dvar = f)
     296    -x*e^x*f(0) + x*e^x*D[0](f)(0) + e^x*f(0)
    301297
    302298.. index:
    303299   pair: differential equations; plot
     
    314310.. math::
    315311    x' = x+y, x(0) = 1; y' = x-y, y(0) = -1,
    316312
    317 for :math:`0 <= t <= 2`. Another way this system can be solved is to use the
    318 command ``desolve_system`` in ``calculus/examples``.
     313for :math:`0 <= t <= 2`. The same result can be obtained by using ``desolve_system_rk4``::
     314
     315    sage: var('x y t')
     316    sage: P=desolve_system_rk4([x+y, x-y], [x,y], ics=[0,1,-1], ivar=t, end_points=2)
     317    sage: p1 = list_plot([[i,j] for i,j,k in P], plotjoined=True)
     318    sage: p2 = list_plot([[i,k] for i,j,k in P], plotjoined=True, color='red')
     319    sage: p1+p2
     320
     321Another way this system can be solved is to use the command ``desolve_system``.
    319322
    320323.. skip
    321324
    322325::
    323326
    324     sage: attach os.environ['SAGE_ROOT'] + '/examples/calculus/desolvers.sage'
    325     sage: des = ["'diff(x(t),t)=x(t)+y(t)","'diff(y(t),t)=x(t)-y(t)"]
    326     sage: vars = ["t","x","y"]
    327     sage: ics = [0,1,-1]
    328     sage: desolve_system(des,vars,ics)
    329     [x(t)=cosh(2^(1/2)*t),y(t)=2*sinh(2^(1/2)*t)/sqrt(2)-cosh(2^(1/2)*t)]
     327    sage: t=var('t'); x=function('x',t); y=function('y',t)
     328    sage: des = [diff(x,t) == x+y, diff(y,t) == x-y]
     329    sage: desolve_system(des, [x,y], ics = [0, 1, -1])
     330    [x(t) == cosh(sqrt(2)*t), y(t) == sqrt(2)*sinh(sqrt(2)*t) - cosh(sqrt(2)*t)]
    330331
    331332The output of this command is *not* a pair of functions.
    332333
  • sage/calculus/calculus.py

    diff -r b89cf68ccc1c -r ea16778cc7ae sage/calculus/calculus.py
    a b  
    514514    If ``self`` has only one variable, then it returns the
    515515    integral with respect to that variable.
    516516
    517      INPUT:
    518  
     517    If definite integration fails, it could be still possbile to
     518    evaluate definite integral using indefinite integral and Newton -
     519    Leibniz theorem (however, the user has to ensure that the
     520    indefinite integral is continuous on compact interval `[a,b]` and
     521    this theorem can be applied).
     522
     523    INPUT:
    519524 
    520525    - ``v`` - (optional) a variable or variable name.  This can also
    521        be a tuple of the variable (optional) and endpoints (i.e.,
    522        ``(x,0,1)`` or ``(0,1)``).
     526      be a tuple of the variable (optional) and endpoints (i.e.,
     527      ``(x,0,1)`` or ``(0,1)``).
    523528 
    524529    - ``a`` - (optional) lower endpoint of definite integral
    525530 
     
    533538               
    534539       - 'mathematica_free' - use http://integrals.wolfram.com/
    535540 
    536      EXAMPLES::
     541    EXAMPLES::
    537542
    538543        sage: x = var('x')
    539544        sage: h = sin(x)/(cos(x))^2
  • sage/calculus/desolvers.py

    diff -r b89cf68ccc1c -r ea16778cc7ae sage/calculus/desolvers.py
    a b  
    11"""
    22This file contains functions useful for solving differential equations
    3 which occur commonly in a 1st semester differential equations course.
     3which occur commonly in a 1st semester differential equations
     4course. For another numerical solver see :meth:`ode_solver` function
     5and optional package Octave.
    46
    5 * desolve -- Computes the "general solution" to a 1st or 2nd order
    6              ODE via Maxima.
     7Commands:
    78
    8 * desolve_laplace -- Solves an ODE using laplace transforms via Maxima.
    9                      Initials conditions are optional.
     9- ``desolve`` - Computes the "general solution" to a 1st or 2nd order
     10  ODE via Maxima.
    1011
    11 * desolve_system -- Solves any size system of 1st order odes using Maxima.
    12                     Initials conditions are optional.
     12- ``desolve_laplace`` - Solves an ODE using laplace transforms via
     13  Maxima. Initials conditions are optional.
    1314
    14 * desolve_rk4 -- Solves numerically IVP for one first order equation,
    15                  returns list of points or plot
     15- ``desolve_system`` - Solves any size system of 1st order odes using
     16  Maxima. Initials conditions are optional.
    1617
    17 * desolve_system_rk4 -- Solves numerically IVP for system of first order
    18                         equations, returns list of points
     18- ``desolve_rk4`` - Solves numerically IVP for one first order
     19  equation, returns list of points or plot
    1920
    20 * eulers_method -- Approximate solution to a 1st order DE, presented as a table.
     21- ``desolve_system_rk4`` - Solves numerically IVP for system of first
     22  order equations, returns list of points
    2123
    22 * eulers_method_2x2 -- Approximate solution to a 1st order system of DEs,
    23                        presented as a table.
     24- ``eulers_method`` - Approximate solution to a 1st order DE,
     25  presented as a table.
    2426
    25 * eulers_method_2x2_plot -- Plots the sequence of points obtained from
    26                     Euler's method.
     27- ``eulers_method_2x2`` - Approximate solution to a 1st order system
     28  of DEs, presented as a table.
    2729
    28 AUTHORS: David Joyner (3-2006)     -- Initial version of functions
    29          Marshall Hampton (7-2007) -- Creation of Python module and testing
    30          Robert Bradshaw (10-2008) -- Some interface cleanup.
    31          Robert Marik (10-2009) -- Some bugfixes and enhancements
     30- ``eulers_method_2x2_plot`` - Plots the sequence of points obtained
     31  from Euler's method.
     32
     33AUTHORS:
     34
     35- David Joyner (3-2006) - Initial version of functions
     36
     37- Marshall Hampton (7-2007) - Creation of Python module and testing
     38
     39- Robert Bradshaw (10-2008) - Some interface cleanup.
     40
     41- Robert Marik (10-2009) - Some bugfixes and enhancements
    3242
    3343"""
    3444
     
    5262def desolve(de, dvar, ics=None, ivar=None, show_method=False, contrib_ode=False):
    5363    """
    5464    Solves a 1st or 2nd order linear ODE via maxima. Including IVP and BVP.
     65   
     66    *Use* ``desolve? <tab>`` *if the output in truncated in notebook.*
    5567
    5668    INPUT:
    57         de    -- an expression or equation representing the ODE
    58         dvar  -- the dependent variable (hereafter called y)
    59         ics   -- (optional) the initial or boundary conditions
    60                     for a first-order equation, specify the initial x and y
    61                     for a second-order equation, specify the initial x, y, and dy/dx
    62                     for a second-order boundary solution, specify initial and final x and y
    63                     initial conditions gives error if the solution is not SymbolisEquation
    64                     (as happens for example for clairot equation)
    65         ivar  -- (optional) the independent variable  (hereafter called x),
     69   
     70    - ``de`` - an expression or equation representing the ODE
     71   
     72    - ``dvar`` - the dependent variable (hereafter called ``y``)
     73
     74    - ``ics`` - (optional) the initial or boundary conditions
     75
     76      - for a first-order equation, specify the initial ``x`` and ``y``
     77
     78      - for a second-order equation, specify the initial ``x``, ``y``,
     79        and ``dy/dx``
     80     
     81      - for a second-order boundary solution, specify initial and
     82        final ``x`` and ``y`` initial conditions gives error if the
     83        solution is not SymbolisEquation (as happens for example for
     84        clairot equation)
     85       
     86    - ``ivar`` - (optional) the independent variable  (hereafter called x),
    6687                    which must be specified if there is more than one
    6788                    independent variable in the equation.
    68         show_method  -- (optional) if true, then Sage returns pair [solution, method], where method
    69                            is the string describing method which has been used to get solution
    70                            (Maxima uses the following order for first order equations: linear, separable,
    71                            exact (including exact with integrating factor), homogeneous, bernoulli,
    72                            generalized homogeneous) - use carefully in class, see below for the example
    73                            of the equation which is separable but this property is not recognized
    74                            by Maxima and equation is solved as exact.
    75         contrib_ode  -- (optional) if true, desolve allows to solve clairot, lagrange, riccati and
    76                            some other equations. May take a long time and thus turned off by default.
    77                            Initial conditions can be used only if the result is one SymbolicEquation
    78                            (does not contain singular solution, for example)
     89                   
     90    - ``show_method`` - (optional) if true, then Sage returns pair
     91      ``[solution, method]``, where method is the string describing
     92      method which has been used to get solution (Maxima uses the
     93      following order for first order equations: linear, separable,
     94      exact (including exact with integrating factor), homogeneous,
     95      bernoulli, generalized homogeneous) - use carefully in class,
     96      see below for the example of the equation which is separable but
     97      this property is not recognized by Maxima and equation is solved
     98      as exact.
     99
     100    - ``contrib_ode`` - (optional) if true, desolve allows to solve
     101      clairot, lagrange, riccati and some other equations. May take
     102      a long time and thus turned off by default.  Initial conditions
     103      can be used only if the result is one SymbolicEquation (does not
     104      contain singular solution, for example)
    79105
    80106
    81107    OUTPUT:
    82         In most cases returns SymbolicEquation which defines the solution implicitly.
    83         If the result is in the form y(x)=... (happens for linear eqs.),
    84         returns the right-hand side only.
    85         The possible constant solutions of separable ODE's are omitted.
     108   
     109    In most cases returns SymbolicEquation which defines the solution
     110    implicitly.  If the result is in the form y(x)=... (happens for
     111    linear eqs.), returns the right-hand side only.  The possible
     112    constant solutions of separable ODE's are omitted.
    86113
    87114
    88     EXAMPLES:
     115    EXAMPLES::
     116
    89117        sage: x = var('x')
    90118        sage: y = function('y', x)
    91119        sage: desolve(diff(y,x) + y - 1, y)
    92120        (c + e^x)*e^(-x)
     121   
     122    ::
     123
    93124        sage: f = desolve(diff(y,x) + y - 1, y, ics=[10,2]); f
    94125        (e^10 + e^x)*e^(-x)
     126
     127    ::
     128
    95129        sage: plot(f)
    96130       
    97     We can also solve second-order differential equations.
     131    We can also solve second-order differential equations.::
    98132   
    99133        sage: x = var('x')
    100134        sage: y = function('y', x)
    101135        sage: de = diff(y,x,2) - y == x
    102136        sage: desolve(de, y)
    103137        k1*e^x + k2*e^(-x) - x
     138       
     139
     140    ::
     141
    104142        sage: f = desolve(de, y, [10,2,1]); f
    105143        -x + 5*e^(-x + 10) + 7*e^(x - 10)
    106144        sage: f(x=10)
     
    108146        sage: diff(f,x)(x=10)
    109147        1
    110148
     149    ::
     150
    111151        sage: de = diff(y,x,2) + y == 0
    112152        sage: desolve(de, y)
    113153        k1*sin(x) + k2*cos(x)
    114154        sage: desolve(de, y, [0,1,pi/2,4])
    115155        4*sin(x) + cos(x)
    116 
     156       
    117157        sage: desolve(y*diff(y,x)+sin(x)==0,y)
    118158        -1/2*y(x)^2 == c - cos(x)
    119159
    120     Clairot equation: general and sungular solutions
     160    Clairot equation: general and singular solutions::
    121161
    122162        sage: desolve(diff(y,x)^2+x*diff(y,x)-y==0,y,contrib_ode=True,show_method=True)
    123163        [[y(x) == c^2 + c*x, y(x) == -1/4*x^2], 'clairault']
    124  
    125     For equations involving more variables we specify independent variable
     164    
     165    For equations involving more variables we specify independent variable::
    126166
    127167        sage: a,b,c,n=var('a b c n')
    128168        sage: desolve(x^2*diff(y,x)==a+b*x^n+c*x^2*y^2,y,ivar=x,contrib_ode=True)
     
    131171        [[[y(x) == 0, (b*x^(n - 2) + a/x^2)*c^2*u == 0]], 'riccati']
    132172
    133173
    134     Higher orded, not involving independent variable
     174    Higher orded, not involving independent variable::
    135175
    136176        sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y).expand()
    137177        1/6*y(x)^3 + k1*y(x) == k2 + x
     
    140180        sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,1,3],show_method=True)
    141181        [1/6*y(x)^3 - 5/3*y(x) == x - 3/2, 'freeofx']
    142182
    143     Separable equations - Sage returns solution in implicit form
     183    Separable equations - Sage returns solution in implicit form::
    144184
    145185        sage: desolve(diff(y,x)*sin(y) == cos(x),y)
    146186        -cos(y(x)) == c + sin(x)
     
    149189        sage: desolve(diff(y,x)*sin(y) == cos(x),y,[pi/2,1])
    150190        -cos(y(x)) == sin(x) - cos(1) - 1
    151191
    152     Linear equation - Sage returns the expression on the right hand side only
     192    Linear equation - Sage returns the expression on the right hand side only::
    153193
    154194        sage: desolve(diff(y,x)+(y) == cos(x),y)
    155195        1/2*((sin(x) + cos(x))*e^x + 2*c)*e^(-x)
     
    158198        sage: desolve(diff(y,x)+(y) == cos(x),y,[0,1])
    159199        1/2*(e^x*sin(x) + e^x*cos(x) + 1)*e^(-x)
    160200
    161     This ODE with separated variables is solved as exact
    162     Explanation - factor does not split exp(x-y) in maxima into exp(x)*exp(-y)
     201    This ODE with separated variables is solved as
     202    exact. Explanation - factor does not split `e^{x-y}` in Maxima
     203    into `e^{x}e^{y}`::
    163204
    164205        sage: desolve(diff(y,x)==exp(x-y),y,show_method=True)
    165206        [-e^x + e^y(x) == c, 'exact']
     
    167208    You can solve Bessel equations. You can also use initial
    168209    conditions, but you cannot put (sometimes desired) initial
    169210    condition at x=0, since this point is singlar point of the
    170     eruation. Anyway, if the solution should be bounded at x=0, then
    171     k2=0.
     211    equation. Anyway, if the solution should be bounded at x=0, then
     212    k2=0.::
    172213           
    173214        sage: desolve(x^2*diff(y,x,x)+x*diff(y,x)+(x^2-4)*y==0,y)
    174215        k1*bessel_j(2, x) + k2*bessel_y(2, x)
    175216   
    176     This example asks (silly) question, whether sqrt(3) is an integer.
    177     has been reported to Maxima developers.
    178         # sage: desolve(x^2*diff(y,x,x)+x*diff(y,x)+(x^2-3)*y==0,y)
     217    Difficult ODE produces error::
    179218
     219        sage: desolve(sqrt(y)*diff(y,x)+e^(y)+cos(x)-sin(x+y)==0,y) # not tested
     220        Traceback (click to the left for traceback)
     221        ...
     222        NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True."
     223       
     224    Difficult ODE produces error - moreover, takes a long time ::
    180225
    181     Produces error (difficult ODE)
    182         #sage: desolve(sqrt(y)*diff(y,x)+e^(y)+cos(x)-sin(x+y)==0,y)
    183         #Traceback (click to the left for traceback)
    184         #...
    185         #NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True."
     226        sage: desolve(sqrt(y)*diff(y,x)+e^(y)+cos(x)-sin(x+y)==0,y,contrib_ode=True) # not tested
    186227
    187         #Produces error (difficult ODE) - moreover, takes a long time
    188         #sage: desolve(sqrt(y)*diff(y,x)+e^(y)+cos(x)-sin(x+y)==0,y,contrib_ode=True)
    189 
    190     Some more types od ODE's
     228    Some more types od ODE's::
    191229
    192230        sage: desolve(x*diff(y,x)^2-(1+x*y)*diff(y,x)+y==0,y,contrib_ode=True,show_method=True)
    193231        [[y(x) == c + log(x), y(x) == c*e^x], 'factor']
    194232
     233    ::
     234   
    195235        sage: desolve(diff(y,x)==(x+y)^2,y,contrib_ode=True,show_method=True)
    196236        [[[x == c - arctan(sqrt(t)), y(x) == -x - sqrt(t)], [x == c + arctan(sqrt(t)), y(x) == -x + sqrt(t)]], 'lagrange']
    197237
    198     These two examples produce error (as expected, Maxima 5.18 cannot solve equations from initial conditions)
    199     Current Maxima 5.18 returns false answer in this case!
     238    These two examples produce error (as expected, Maxima 5.18 cannot
     239    solve equations from initial conditions). Current Maxima 5.18
     240    returns false answer in this case!::
    200241
    201         #sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,2]).expand()
    202         #Traceback (click to the left for traceback)
    203         #...
    204         #NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True."
    205         #sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,2],show_method=True)
    206         #Traceback (click to the left for traceback)
    207         #...
    208         #NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True."
     242        sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,2]).expand() # not tested
     243        Traceback (click to the left for traceback)
     244        ...
     245        NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True."
     246       
     247    ::
    209248
     249        sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,2],show_method=True) # not tested
     250        Traceback (click to the left for traceback)
     251        ...
     252        NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True."
    210253
    211     Second order linear ODE
     254    Second order linear ODE::
    212255
    213256        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y)
    214257        (k2*x + k1)*e^(-x) + 1/2*sin(x)
     
    222265        3*((e^(1/2*pi) - 2)*x/pi + 1)*e^(-x) + 1/2*sin(x)
    223266        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,[0,3,pi/2,2],show_method=True)
    224267        [3*((e^(1/2*pi) - 2)*x/pi + 1)*e^(-x) + 1/2*sin(x), 'variationofparameters']
    225 
     268       
    226269        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y)
    227270        (k2*x + k1)*e^(-x)
    228271        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,show_method=True)
     
    237280        [(2*(2*e^(1/2*pi) - 3)*x/pi + 3)*e^(-x), 'constcoeff']
    238281
    239282       
    240     This equation can be solved within Maxima but not within Sage
    241     needs assumptions assume(x>0,y>0) and works in maxima, but not in Sage
    242         # sage: assume(x>0)
    243         # sage: assume(y>0)
    244         # sage: desolve(x*diff(y,x)-x*sqrt(y^2+x^2)-y,y,show_method=True)   
     283    This equation can be solved within Maxima but not within Sage. It
     284    needs assumptions assume(x>0,y>0) and works in Maxima, but not in Sage::
    245285
    246     This equation can be solved within Maxima but not within Sage
    247     Perhaps upgrading Maxima help
    248         # sage: desolve(x*(x+1)*diff(y,x,2)+(x+5)*diff(y,x,1)+(-4)*y,y,contrib_ode=True)
     286        sage: assume(x>0) # not tested
     287        sage: assume(y>0) # not tested
     288        sage: desolve(x*diff(y,x)-x*sqrt(y^2+x^2)-y,y,show_method=True) # not tested
    249289
    250290       
    251291    TESTS:
    252292   
    253     Trac #6479 fixed
     293    Trac #6479 fixed::
    254294
    255295        sage: x = var('x')
    256296        sage: y = function('y', x)
    257297        sage: desolve( diff(y,x,x) == 0, y, [0,0,1])
    258298        x
     299
     300    ::
     301   
    259302        sage: desolve( diff(y,x,x) == 0, y, [0,1,1])
    260303        x + 1
    261304
    262305
    263     AUTHOR: David Joyner (1-2006)
    264             Robert Bradshaw (10-2008)
    265             Robert Marik (10-2009)
     306    AUTHORS:
     307 
     308    - David Joyner (1-2006)
     309
     310    - Robert Bradshaw (10-2008)
     311
     312    - Robert Marik (10-2009)
    266313
    267314    """
    268315    if is_SymbolicEquation(de):
     
    411458    Solves an ODE using laplace transforms. Initials conditions are optional.
    412459
    413460    INPUT:
    414         de    -- a lambda expression representing the ODE
    415                  (eg, de = diff(y,x,2) == diff(y,x)+sin(x))
    416         dvar  -- the dependent variable (eg y)
    417         ivar  -- (optional) the independent variable  (hereafter called x),
    418                  which must be specified if there is more than one
    419                  independent variable in the equation.
    420         ics   -- a list of numbers representing initial conditions,
    421                  (eg, f(0)=1, f'(0)=2 is ics = [0,1,2])
     461   
     462    - ``de`` - a lambda expression representing the ODE (eg, de =
     463      diff(y,x,2) == diff(y,x)+sin(x))
     464
     465    - ``dvar`` - the dependent variable (eg y)
     466
     467    - ``ivar`` - (optional) the independent variable (hereafter called
     468      x), which must be specified if there is more than one
     469      independent variable in the equation.
     470
     471    - ``ics`` - a list of numbers representing initial conditions, (eg,
     472      f(0)=1, f'(0)=2 is ics = [0,1,2])
    422473
    423474    OUTPUT:
    424         Solution of the ODE as symbolic expression
    425475
    426     EXAMPLES:
     476    Solution of the ODE as symbolic expression
     477
     478    EXAMPLES::
     479   
    427480        sage: u=function('u',x)
    428481        sage: eq = diff(u,x) - exp(-x) - u == 0
    429482        sage: desolve_laplace(eq,u)
    430483        1/2*(2*u(0) + 1)*e^x - 1/2*e^(-x)
    431484   
    432     We can use initial conditions
     485    We can use initial conditions::
    433486
    434487        sage: desolve_laplace(eq,u,ics=[0,3])
    435488        -1/2*e^(-x) + 7/2*e^x
    436489
    437490    The initial conditions do not persist in the system (as they persisted
    438     in previous versions)
     491    in previous versions)::
    439492
    440493        sage: desolve_laplace(eq,u)
    441494        1/2*(2*u(0) + 1)*e^x - 1/2*e^(-x)
    442495
     496    ::
     497
    443498        sage: f=function('f', x)
    444499        sage: eq = diff(f,x) + f == 0
    445500        sage: desolve_laplace(eq,f,[0,1])
    446501        e^(-x)
    447 
     502   
     503    ::
     504   
    448505        sage: x = var('x')
    449506        sage: f = function('f', x)
    450507        sage: de = diff(f,x,x) - 2*diff(f,x) + f
     
    454511        x*e^x + e^x
    455512
    456513    TESTS:
    457         #4839 fixed
     514
     515    Trac #4839 fixed::
    458516
    459517        sage: t=var('t')
    460518        sage: x=function('x', t)
     
    463521        e^(-t) + 1
    464522        sage: soln(t=3)
    465523        e^(-3) + 1
     524   
     525    AUTHORS:
    466526
    467     AUTHOR: David Joyner (1-2006,8-2007)
    468             Robert Marik (10-2009)
     527    - David Joyner (1-2006,8-2007)
     528   
     529    - Robert Marik (10-2009)
    469530    """
    470531    #This is the original code from David Joyner (inputs and outputs strings)
    471532    #maxima("de:"+de._repr_()+"=0;")
     
    514575    Solves any size system of 1st order ODE's. Initials conditions are optional.
    515576
    516577    INPUT:
    517         des   -- list of ODEs
    518         vars  -- list of dependent variables
    519         ics   -- (optional) list of initial values for ivar and vars
    520         ivar  -- (optional) the independent variable, which must be specified
    521                     if there is more than one independent variable in the equation.
     578   
     579    - ``des`` - list of ODEs
     580   
     581    - ``vars`` - list of dependent variables
     582   
     583    - ``ics`` - (optional) list of initial values for ivar and vars
     584   
     585    - ``ivar`` - (optional) the independent variable, which must be
     586      specified if there is more than one independent variable in the
     587      equation.
    522588
    523     EXAMPLES:
     589    EXAMPLES::
     590
    524591        sage: t = var('t')
    525592        sage: x = function('x', t)
    526593        sage: y = function('y', t)
     
    530597        [x(t) == (x(0) - 1)*cos(t) - (y(0) - 1)*sin(t) + 1,
    531598         y(t) == (x(0) - 1)*sin(t) + (y(0) - 1)*cos(t) + 1]
    532599         
    533     Now we give some initial conditions:
     600    Now we give some initial conditions::
     601   
    534602        sage: sol = desolve_system([de1, de2], [x,y], ics=[0,1,2]); sol
    535603        [x(t) == -sin(t) + 1, y(t) == cos(t) + 1]
    536604        sage: solnx, solny = sol[0].rhs(), sol[1].rhs()
    537605        sage: plot([solnx,solny],(0,1))
    538606        sage: parametric_plot((solnx,solny),(0,1))
    539607
    540     AUTHOR: Robert Bradshaw (10-2008)
     608    AUTHORS:
     609
     610    - Robert Bradshaw (10-2008)
    541611    """
    542612    ivars = set([])
    543613    for i, de in enumerate(des):
     
    611681        Now type show(P1), show(P2) to view these.
    612682                     
    613683
    614     AUTHOR: David Joyner (3-2006, 8-2007)
     684    AUTHORS: David Joyner (3-2006, 8-2007)
    615685    """
    616686    d = len(des)
    617687    dess = [de._maxima_init_() + "=0" for de in des]
     
    633703
    634704def eulers_method(f,x0,y0,h,x1,method="table"):
    635705    """
    636     This implements Euler's method for finding numerically the solution of the 1st order
    637     ODE y' = f(x,y), y(a)=c. The "x" column of the table increments from
    638     x0 to x1 by h (so (x1-x0)/h must be an integer). In the "y" column,
    639     the new y-value equals the old y-value plus the corresponding entry in the
     706    This implements Euler's method for finding numerically the
     707    solution of the 1st order ODE ``y' = f(x,y)``, ``y(a)=c``. The "x"
     708    column of the table increments from ``x0`` to ``x1`` by ``h`` (so
     709    ``(x1-x0)/h`` must be an integer). In the "y" column, the new
     710    y-value equals the old y-value plus the corresponding entry in the
    640711    last column.
    641712
    642713    *For pedagogical purposes only.*
    643714   
    644     EXAMPLES:
     715    EXAMPLES::
     716   
    645717        sage: from sage.calculus.desolvers import eulers_method
    646718        sage: x,y = PolynomialRing(QQ,2,"xy").gens()
    647719        sage: eulers_method(5*x+y-5,0,1,1/2,1)
     
    649721             0                    1                   -2
    650722           1/2                   -1                 -7/4
    651723             1                -11/4                -11/8
     724
     725    ::
     726
    652727        sage: x,y = PolynomialRing(QQ,2,"xy").gens()
    653728        sage: eulers_method(5*x+y-5,0,1,1/2,1,method="none")
    654729        [[0, 1], [1/2, -1], [1, -11/4], [3/2, -33/8]]
     730
     731    ::
     732
    655733        sage: RR = RealField(sci_not=0, prec=4, rnd='RNDU')
    656734        sage: x,y = PolynomialRing(RR,2,"xy").gens()
    657735        sage: eulers_method(5*x+y-5,0,1,1/2,1,method="None")
    658736        [[0, 1], [1/2, -1.0], [1, -2.7], [3/2, -4.0]]
     737
     738    ::
     739
    659740        sage: RR = RealField(sci_not=0, prec=4, rnd='RNDU')
    660741        sage: x,y=PolynomialRing(RR,2,"xy").gens()
    661742        sage: eulers_method(5*x+y-5,0,1,1/2,1)
     
    663744             0                    1                 -2.0
    664745           1/2                 -1.0                 -1.7
    665746             1                 -2.7                 -1.3
     747
     748    ::
     749
    666750        sage: x,y=PolynomialRing(QQ,2,"xy").gens()
    667751        sage: eulers_method(5*x+y-5,1,1,1/3,2)
    668752                 x                    y                  h*f(x,y)
     
    670754               4/3                  4/3                    1
    671755               5/3                  7/3                 17/9
    672756                 2                 38/9                83/27
     757
     758    ::
     759
    673760        sage: eulers_method(5*x+y-5,0,1,1/2,1,method="none")
    674761        [[0, 1], [1/2, -1], [1, -11/4], [3/2, -33/8]]
     762
     763    ::
     764
    675765        sage: pts = eulers_method(5*x+y-5,0,1,1/2,1,method="none")
    676766        sage: P1 = list_plot(pts)
    677767        sage: P2 = line(pts)
    678768        sage: (P1+P2).show()
    679769
    680     AUTHOR: David Joyner
     770    AUTHORS:
     771
     772    - David Joyner
    681773    """
    682774    if method=="table":
    683775        print "%10s %20s %25s"%("x","y","h*f(x,y)")
     
    695787
    696788def eulers_method_2x2(f,g, t0, x0, y0, h, t1,method="table"):
    697789    """
    698     This implements Euler's method for finding numerically the solution of the 1st order
    699     system of two ODEs
    700         x' = f(t, x, y), x(t0)=x0.
    701         y' = g(t, x, y), y(t0)=y0.
    702     The "t" column of the table increments from t0 to t1 by
    703     h (so (t1-t0)/h must be an integer). In the "x" column, the new x-value equals the
    704     old x-value plus the corresponding entry in the next (third) column.
    705     In the "y" column, the new y-value equals the old y-value plus the
    706     corresponding entry in the next (last) column.
     790    This implements Euler's method for finding numerically the
     791    solution of the 1st order system of two ODEs
     792
     793    ``x' = f(t, x, y), x(t0)=x0.``
     794
     795    ``y' = g(t, x, y), y(t0)=y0.``
     796
     797    The "t" column of the table increments from `t_0` to `t_1` by `h`
     798    (so `\\frac{t_1-t_0}{h}` must be an integer). In the "x" column,
     799    the new x-value equals the old x-value plus the corresponding
     800    entry in the next (third) column.  In the "y" column, the new
     801    y-value equals the old y-value plus the corresponding entry in the
     802    next (last) column.
    707803
    708804    *For pedagogical purposes only.*
    709805   
    710     EXAMPLES:
     806    EXAMPLES::
     807
    711808        sage: from sage.calculus.desolvers import eulers_method_2x2
    712809        sage: t, x, y = PolynomialRing(QQ,3,"txy").gens()
    713810        sage: f = x+y+t; g = x-y
    714811        sage: eulers_method_2x2(f,g, 0, 0, 0, 1/3, 1,method="none")
    715812        [[0, 0, 0], [1/3, 0, 0], [2/3, 1/9, 0], [1, 10/27, 1/27], [4/3, 68/81, 4/27]]
     813
     814    ::
     815
    716816        sage: eulers_method_2x2(f,g, 0, 0, 0, 1/3, 1)
    717817             t                    x                h*f(t,x,y)                    y           h*g(t,x,y)
    718818             0                    0                         0                    0                    0
    719819           1/3                    0                       1/9                    0                    0
    720820           2/3                  1/9                      7/27                    0                 1/27
    721821             1                10/27                     38/81                 1/27                  1/9
     822
     823    ::
     824
    722825        sage: RR = RealField(sci_not=0, prec=4, rnd='RNDU')
    723826        sage: t,x,y=PolynomialRing(RR,3,"txy").gens()
    724827        sage: f = x+y+t; g = x-y
     
    729832           2/3                 0.13                      0.29                 0.00                0.043
    730833             1                 0.41                      0.57                0.043                 0.15
    731834
    732     To numerically approximate y(1), where (1+t^2)y''+y'-y=0, y(0)=1,y'(0)=-1,
    733     using 4 steps of Euler's method, first convert to a system:
    734     y1' = y2, y1(0)=1; y2' = (y1-y2)/(1+t^2), y2(0)=-1.
     835    To numerically approximate `y(1)`, where `(1+t^2)y''+y'-y=0`,
     836    `y(0)=1`, `y'(0)=-1`, using 4 steps of Euler's method, first
     837    convert to a system: `y_1' = y_2`, `y_1(0)=1`; `y_2' =
     838    \\frac{y_1-y_2}{1+t^2}`, `y_2(0)=-1`.::
    735839
    736840         sage: RR = RealField(sci_not=0, prec=4, rnd='RNDU')
    737841         sage: t, x, y=PolynomialRing(RR,3,"txy").gens()
     
    744848           3/4                 0.63                   -0.0078               -0.031                 0.11
    745849             1                 0.63                     0.020                0.079                0.071
    746850
    747     To numerically approximate y(1), where y''+ty'+y=0, y(0)=1,y'(0)=0:
     851    To numerically approximate y(1), where `y''+ty'+y=0`, `y(0)=1`, `y'(0)=0`::
    748852
    749853        sage: t,x,y=PolynomialRing(RR,3,"txy").gens()
    750854        sage: f = y; g = -x-y*t
     
    756860           3/4                 0.88                     -0.15                -0.62                -0.10
    757861             1                 0.75                     -0.17                -0.68               -0.015
    758862
    759     AUTHOR: David Joyner
     863    AUTHORS:
     864
     865    - David Joyner
    760866    """
    761867    if method=="table":
    762868        print "%10s %20s %25s %20s %20s"%("t", "x","h*f(t,x,y)","y", "h*g(t,x,y)")
     
    776882
    777883def eulers_method_2x2_plot(f,g, t0, x0, y0, h, t1):
    778884    """
    779     This plots the soln in the rectangle
    780     (xrange[0],xrange[1]) x (yrange[0],yrange[1])
    781     and plots using Euler's method the numerical solution of
    782     the 1st order ODEs x' = f(t,x,y), x(a)=x0, y' = g(t,x,y), y(a) = y0.
     885    This plots the soln in the rectangle ``(xrange[0],xrange[1])
     886    x (yrange[0],yrange[1])`` and plots using Euler's method the
     887    numerical solution of the 1st order ODEs `x' = f(t,x,y)`,
     888    `x(a)=x_0`, `y' = g(t,x,y)`, `y(a) = y_0`.
    783889
    784890    *For pedagogical purposes only.*
    785891   
    786     ===
    787     The following example plots the solution to theta''+sin(theta)=0,
    788     theta(0)=3/4, theta'(0) = 0.  Type P[0].show() to plot the solution,
    789     (P[0]+P[1]).show() to plot (t,theta(t)) and (t,theta'(t)).
    790     ===
    791     EXAMPLES:
     892    EXAMPLES::
     893
     894    The following example plots the solution to
     895    `\\theta''+\\sin(\\theta)=0`, `\\theta(0)=\\frac 34`, `\\theta'(0) =
     896    0`.  Type ``P[0].show()`` to plot the solution,
     897    ``(P[0]+P[1]).show()`` to plot `(t,\\theta(t))` and
     898    `(t,\\theta'(t))`::
     899
    792900        sage: from sage.calculus.desolvers import eulers_method_2x2_plot
    793901        sage: f = lambda z : z[2]; g = lambda z : -sin(z[1])
    794902        sage: P = eulers_method_2x2_plot(f,g, 0.0, 0.75, 0.0, 0.1, 1.0)
     
    810918    """
    811919    Used to determine bounds for numerical integration.
    812920   
    813     If end_points is None, the interval for integration is from ics[0]
    814     to ics[0]+10
     921    - If end_points is None, the interval for integration is from ics[0]
     922      to ics[0]+10
    815923
    816     If end_points is a or [a], the interval for integration is from min(ics[0],a)
    817     to max(ics[0],a)
     924    - If end_points is a or [a], the interval for integration is from min(ics[0],a)
     925      to max(ics[0],a)
    818926
    819     If end_points is [a,b], the interval for integration is from min(ics[0],a)
    820     to max(ics[0],b)
     927    - If end_points is [a,b], the interval for integration is from min(ics[0],a)
     928      to max(ics[0],b)
    821929
    822    EXAMPLES:
     930    EXAMPLES::
     931
    823932        sage: from sage.calculus.desolvers import desolve_rk4_determine_bounds
    824933        sage: desolve_rk4_determine_bounds([0,2],1)
    825934        (0, 1)
     935   
     936    ::
     937
    826938        sage: desolve_rk4_determine_bounds([0,2]) 
    827939        (0, 10)
     940   
     941    ::
     942
    828943        sage: desolve_rk4_determine_bounds([0,2],[-2])
    829944        (-2, 0)
     945   
     946    ::
     947
    830948        sage: desolve_rk4_determine_bounds([0,2],[-2,4])
    831949        (-2, 4)
    832950   
     
    843961
    844962def desolve_rk4(de, dvar, ics=None, ivar=None, end_points=None, step=0.1, output='list', **kwds):
    845963    """
    846     Solves numerically one first-order ordinary differential equation.
     964    Solves numerically one first-order ordinary differential
     965    equation. See also ``ode_solver``.
    847966
    848967    INPUT:
    849         input is similar to desolve command. The differential equation can be
    850         written in a form close to the plot_slope_field or desolve command
     968   
     969    input is similar to ``desolve`` command. The differential equation can be
     970    written in a form close to the plot_slope_field or desolve command
    851971
    852         Variant 1 (function in two variables)
    853         de   -- right hand side, i.e. the function f(x,y) from ODE y'=f(x,y)
    854         dvar -- dependent variable (symbolic variable declared by var)
     972    - Variant 1 (function in two variables)
     973     
     974      - ``de`` - right hand side, i.e. the function `f(x,y)` from ODE `y'=f(x,y)`
     975     
     976      - ``dvar`` - dependent variable (symbolic variable declared by var)
    855977
    856         Variant 2 (symbolic equation)
    857         de   -- equation, including term with diff(y,x)
    858         dvar -- dependent variable (declared as funciton of independent variable)
     978    - Variant 2 (symbolic equation)
    859979
    860         Other parameters
    861         ivar  -- should be specified, if there are more variables or if the equation is autonomous
    862         ics   -- initial conditions in the form [x0,y0]
    863         end_points  -- the end points of the interval
    864            if end_points is a or [a], we integrate on between min(ics[0],a) and max(ics[0],a)
    865            if end_points is None, we use end_points=ics[0]+10
    866            if end_points is [a,b] we integrate on between min(ics[0],a) and max(ics[0],b)
    867         step   -- the length of the step (positive number)
    868         output -- 'list', 'plot', 'slope_field' (graph of the solution with slope field)
     980      - ``de`` - equation, including term with ``diff(y,x)``
    869981
     982      - ``dvar``` - dependent variable (declared as funciton of independent variable)
    870983
    871     OUTPUT:
    872         Returns a list of points, or plot produced by list_plot, optionally with slope field.
     984    - Other parameters
    873985
     986      - ``ivar`` - should be specified, if there are more variables or if the equation is autonomous
    874987
    875     EXAMPLES:
     988      - ``ics`` - initial conditions in the form [x0,y0]
     989     
     990      - ``end_points`` - the end points of the interval
     991
     992        - if end_points is a or [a], we integrate on between min(ics[0],a) and max(ics[0],a)
     993        - if end_points is None, we use end_points=ics[0]+10
     994
     995        - if end_points is [a,b] we integrate on between min(ics[0],a) and max(ics[0],b)
     996
     997      - ``step`` - (optional, default:0.1) the length of the step (positive number)
     998 
     999      - ``output`` - (optional, default: 'list') one of 'list',
     1000        'plot', 'slope_field' (graph of the solution with slope field)
     1001
     1002    OUTPUT:
     1003
     1004    Returns a list of points, or plot produced by list_plot,
     1005    optionally with slope field.
     1006
     1007
     1008    EXAMPLES::
     1009
    8761010        sage: from sage.calculus.desolvers import desolve_rk4
    8771011
    878     Variant 2 for input - more common in numerics
     1012    Variant 2 for input - more common in numerics::
    8791013
    8801014        sage: x,y=var('x y')
    8811015        sage: desolve_rk4(x*y*(2-y),y,ics=[0,1],end_points=1,step=0.5)
     
    8831017
    8841018    Variant 1 for input - we can pass ODE in the form used by
    8851019    desolve function In this example we integrate bakwards, since
    886     end_points < ics[0]
     1020    ``end_points < ics[0]``::
    8871021
    8881022        sage: y=function('y',x)
    8891023        sage: desolve_rk4(diff(y,x)+y*(y-1) == x-2,y,ics=[1,1],step=0.5, end_points=0)
    8901024        [[0.0, 8.90425710896], [0.5, 1.90932794536], [1, 1]]
    8911025
    892     Here we show how to plot simple pictures. For more advanced aplications use
    893     list_plot instead. To see the resulting picture use show(P) in Sage notebook.
     1026    Here we show how to plot simple pictures. For more advanced
     1027    aplications use list_plot instead. To see the resulting picture
     1028    use ``show(P)`` in Sage notebook. ::
    8941029
    8951030        sage: x,y=var('x y')
    8961031        sage: P=desolve_rk4(y*(2-y),y,ics=[0,.1],ivar=x,output='slope_field',end_points=[-4,6],thickness=3)
    8971032               
    8981033    ALGORITHM:
    899         4th order Runge-Kutta method. Wrapper for command rk in Maxima's dynamics package.
    900         Perhaps could be faster by using fast_float instead.
    9011034
    902     AUTHOR:
    903         Robert Marik (10-2009)
     1035    4th order Runge-Kutta method. Wrapper for command ``rk`` in
     1036    Maxima's dynamics package.  Perhaps could be faster by using
     1037    fast_float instead.
     1038
     1039    AUTHORS:
     1040
     1041    - Robert Marik (10-2009)
    9041042    """
    9051043    if ics is None:
    9061044        raise ValueError, "No initial conditions, specify with ics=[x0,y0]."
     
    9731111    """
    9741112    Solves numerically system of first-order ordinary differential
    9751113    equations using the 4th order Runge-Kutta method. Wrapper for
    976     Maxima command rk.
     1114    Maxima command ``rk``. See also ``ode_solver``.
    9771115
    9781116    INPUT:
    979         input is similar to desolve_system and desolve_rk4 commands
    980         des  -- right hand sides of the system
    981         vars -- dependent variables
    982         ivar -- should be specified, if there are more variables or if the equation is autonomous
    983                and the independent variable is missing
    984         ics  -- initial conditions in the form [x0,y01,y02,y03,....]
    985         end_points  -- the end points of the interval
    986            if end_points is a or [a], we integrate on between min(ics[0],a) and max(ics[0],a)
    987            if end_points is None, we use end_points=ics[0]+10
    988            if end_points is [a,b] we integrate on between min(ics[0],a) and max(ics[0],b)
    989         step -- the length of the step
     1117   
     1118    input is similar to desolve_system and desolve_rk4 commands
     1119   
     1120    - ``des`` - right hand sides of the system
    9901121
     1122    - ``vars`` - dependent variables
     1123
     1124    - ``ivar`` - (optional) should be specified, if there are more variables or
     1125      if the equation is autonomous and the independent variable is
     1126      missing
     1127
     1128    - ``ics`` - initial conditions in the form [x0,y01,y02,y03,....]
     1129   
     1130    - ``end_points`` - the end points of the interval
     1131   
     1132      - if end_points is a or [a], we integrate on between min(ics[0],a) and max(ics[0],a)
     1133      - if end_points is None, we use end_points=ics[0]+10
     1134
     1135      - if end_points is [a,b] we integrate on between min(ics[0],a) and max(ics[0],b)
     1136
     1137    - ``step`` -- (optional, default: 0.1) the length of the step
    9911138
    9921139    OUTPUT:
    993         Returns a list of points.
    9941140
     1141    Returns a list of points.
    9951142
    996     EXAMPLES:
    997     Lotka Volterra system
     1143    EXAMPLES::
     1144
     1145    Lotka Volterra system::
    9981146
    9991147        sage: from sage.calculus.desolvers import desolve_system_rk4
    10001148        sage: x,y,t=var('x y t')
     
    10061154        sage: LP=list_plot(Q)
    10071155   
    10081156    ALGORITHM:
    1009         4th order Runge-Kutta method. Wrapper for command rk in Maxima's dynamics package.
    1010         Perhaps could be faster by using fast_float instead.
    10111157
    1012     AUTHOR: Robert Marik (10-2009)
     1158    4th order Runge-Kutta method. Wrapper for command ``rk`` in Maxima's
     1159    dynamics package.  Perhaps could be faster by using ``fast_float``
     1160    instead.
     1161
     1162    AUTHOR:
     1163
     1164    - Robert Marik (10-2009)
    10131165    """
    10141166
    10151167    if ics is None:
  • sage/gsl/ode.pyx

    diff -r b89cf68ccc1c -r ea16778cc7ae sage/gsl/ode.pyx
    a b  
    1111include 'gsl.pxi'
    1212
    1313
    14 import sage.gsl.interpolation 
     14import sage.gsl.interpolation
    1515
    1616cdef class PyFunctionWrapper:
    17    cdef object the_function
    18    cdef object the_jacobian
    19    cdef object the_parameters
    20    cdef int y_n
    21    
     17    cdef object the_function
     18    cdef object the_jacobian
     19    cdef object the_parameters
     20    cdef int y_n
    2221
    23    
    24    cdef set_yn(self,x):
    25       self.y_n = x
     22
     23
     24    cdef set_yn(self,x):
     25        self.y_n = x
    2626
    2727
    2828cdef class ode_system:
    29    cdef int  c_j(self,double t, double *y, double *dfdy,double *dfdt): #void *params):
    30       return 0
     29    cdef int  c_j(self,double t, double *y, double *dfdy,double *dfdt): #void *params):
     30        return 0
    3131
    32    cdef int c_f(self,double t, double* y, double* dydt):  #void *params):
    33       return 0
     32    cdef int c_f(self,double t, double* y, double* dydt):  #void *params):
     33        return 0
    3434
    3535
    3636
    3737
    3838cdef int c_jac_compiled(double t, double *y, double *dfdy,double *dfdt, void * params):
    39    cdef int status
    40    cdef ode_system wrapper
    41    wrapper = <ode_system> params
    42    status = wrapper.c_j(t,y,dfdy,dfdt)    #Could add parameters
    43    return status
     39    cdef int status
     40    cdef ode_system wrapper
     41    wrapper = <ode_system> params
     42    status = wrapper.c_j(t,y,dfdy,dfdt)    #Could add parameters
     43    return status
    4444cdef int c_f_compiled(double t, double *y, double *dydt, void *params):
    45    cdef int status
    46    cdef ode_system wrapper
    47    wrapper = <ode_system> params
    48    status =  wrapper.c_f(t,y,dydt)  #Could add parameters
    49    return status
     45    cdef int status
     46    cdef ode_system wrapper
     47    wrapper = <ode_system> params
     48    status =  wrapper.c_f(t,y,dydt)  #Could add parameters
     49    return status
    5050
    5151cdef int c_jac(double t,double *y,double *dfdy,double *dfdt,void *params):
    52  
    53    cdef int i
    54    cdef int j
    55    cdef int y_n
    56    cdef int param_n
    57    cdef PyFunctionWrapper wrapper
    58    wrapper = <PyFunctionWrapper > params
    59    y_n=wrapper.y_n
    60    y_list=[]
    61    for i from 0<=i<y_n:
    62       y_list.append(y[i])
    63    try:
    64       if len(wrapper.the_parameters)==0:
    65          jac_list=wrapper.the_jacobian(t,y_list)
    66       else:
    67          jac_list=wrapper.the_jacobian(t,y_list,wrapper.the_parameters)
    68       for i from 0<=i<y_n:
    69          for j from 0<=j<y_n:
    70             dfdy[i*y_n+j]=jac_list[i][j]
    7152
    72       for i from 0 <=i<y_n:
    73          dfdt[i]=jac_list[y_n][i]
     53    cdef int i
     54    cdef int j
     55    cdef int y_n
     56    cdef int param_n
     57    cdef PyFunctionWrapper wrapper
     58    wrapper = <PyFunctionWrapper > params
     59    y_n=wrapper.y_n
     60    y_list=[]
     61    for i from 0<=i<y_n:
     62        y_list.append(y[i])
     63    try:
     64        if len(wrapper.the_parameters)==0:
     65            jac_list=wrapper.the_jacobian(t,y_list)
     66        else:
     67            jac_list=wrapper.the_jacobian(t,y_list,wrapper.the_parameters)
     68        for i from 0<=i<y_n:
     69            for j from 0<=j<y_n:
     70                dfdy[i*y_n+j]=jac_list[i][j]
    7471
    75       return GSL_SUCCESS
    76    except:
    77       return -1
     72        for i from 0 <=i<y_n:
     73            dfdt[i]=jac_list[y_n][i]
     74
     75        return GSL_SUCCESS
     76    except:
     77        return -1
    7878
    7979cdef int c_f(double t,double* y, double* dydt,void *params):
    80    cdef int i
    81    cdef int y_n
    82    cdef int param_n
     80    cdef int i
     81    cdef int y_n
     82    cdef int param_n
    8383
    84    cdef PyFunctionWrapper wrapper
    85    wrapper = <PyFunctionWrapper> params
    86    y_n= wrapper.y_n
    87    y_list=[]
    88    for i from 0<=i<y_n:
    89        y_list.append(y[i])
    90    try:
    91       if len(wrapper.the_parameters)!=0:
    92          dydt_list=wrapper.the_function(t,y_list,wrapper.the_parameters)
    93       else:
    94          dydt_list=wrapper.the_function(t,y_list)
    95       for i from 0<=i<y_n:
    96           dydt[i]=dydt_list[i]
    97       return GSL_SUCCESS
    98    except:
    99       return -1
     84    cdef PyFunctionWrapper wrapper
     85    wrapper = <PyFunctionWrapper> params
     86    y_n= wrapper.y_n
     87    y_list=[]
     88    for i from 0<=i<y_n:
     89        y_list.append(y[i])
     90    try:
     91        if len(wrapper.the_parameters)!=0:
     92            dydt_list=wrapper.the_function(t,y_list,wrapper.the_parameters)
     93        else:
     94            dydt_list=wrapper.the_function(t,y_list)
     95        for i from 0<=i<y_n:
     96            dydt[i]=dydt_list[i]
     97        return GSL_SUCCESS
     98    except:
     99        return -1
    100100
    101101
    102102
    103103
    104104class ode_solver(object):
    105    """ode_solver is a class that wraps the GSL libraries ode solver routines
    106          To use it instantiate a class,
    107          sage: T=ode_solver()
     105    r"""
     106    :meth:`ode_solver` is a class that wraps the GSL libraries ode
     107    solver routines To use it instantiate a class,::
    108108
    109          To solve a system of the form dy_i/dt=f_i(t,y), you must supply a
    110          vector or tuple/list valued function f representing f_i.
    111          The functions f and the jacobian should have the form foo(t,y) or foo(t,y,params).
    112          params which is optional allows for your function to depend on one or a tuple of parameters.
    113          Note if you use it, params must be a tuple even if it only has one component.
    114          For example if you wanted to solve y''+y=0. You need to write it as a first order system
    115          y_0' = y_1
    116          y1_1 = -y_0
     109        sage: T=ode_solver()
    117110
    118          In code,
    119          
    120             sage: f = lambda t,y:[y[1],-y[0]]
    121             sage: T.function=f
     111    To solve a system of the form ``dy_i/dt=f_i(t,y)``, you must
     112    supply a vector or tuple/list valued function ``f`` representing
     113    ``f_i``.  The functions ``f`` and the jacobian should have the
     114    form ``foo(t,y)`` or ``foo(t,y,params)``.  ``params`` which is
     115    optional allows for your function to depend on one or a tuple of
     116    parameters.  Note if you use it, ``params`` must be a tuple even
     117    if it only has one component.  For example if you wanted to solve
     118    `y''+y=0`. You need to write it as a first order system::
     119   
     120        y_0' = y_1
     121        y_1' = -y_0
     122   
     123    In code::
     124   
     125        sage: f = lambda t,y:[y[1],-y[0]]
     126        sage: T.function=f
    122127
    123          For some algorithms the jacobian must be supplied as well,
    124          the form of this should be a function return a list of lists
    125          of the form
     128    For some algorithms the jacobian must be supplied as well, the
     129    form of this should be a function return a list of lists of the
     130    form ``[ [df_1/dy_1,...,df_1/dy_n], ...,
     131    [df_n/dy_1,...,df_n,dy_n], [df_1/dt,...,df_n/dt] ]``.
    126132
    127          [ [df_1/dy_1,...,df_1/dy_n], ..., [df_n/dy_1,...,df_n,dy_n],
    128          [df_1/dt,...,df_n/dt] ].
     133    There are examples below, if your jacobian was the function
     134    ``my_jacobian`` you would do::
    129135
    130          There are examples below, if your jacobian was the function
    131          my_jacobian you would do.
     136        sage: T.jacobian = my_jacobian     # not tested, since it doesn't make sense to test this
    132137
    133              sage: T.jacobian = my_jacobian     # not tested, since it doesn't make sense to test this
     138    There are a variety of algorithms available for different types of systems. Possible algorithms are
    134139
    135          
    136          There are a variety of algorithms available for different types of systems. Possible algorithms are
    137          "rkf45" - runga-kutta-felhberg (4,5)
    138          "rk2" - embedded runga-kutta (2,3)
    139          "rk4" - 4th order classical runga-kutta
    140          "rk8pd" - runga-kutta prince-dormand (8,9)
    141          "rk2imp" - implicit 2nd order runga-kutta at gaussian points
    142          "rk4imp" - implicit 4th order runga-kutta at gaussian points
    143          "bsimp" - implicit burlisch-stoer (requires jacobian)
    144          "gear1" - M=1 implicit gear
    145          "gear2" - M=2 implicit gear
     140    - ``rkf45`` - runga-kutta-felhberg (4,5)
    146141
    147          The default algorithm if rkf45. If you instead wanted to use bsimp you would do
     142    - ``rk2`` - embedded runga-kutta (2,3)
    148143
    149             sage: T.algorithm="bsimp"
     144    - ``rk4`` - 4th order classical runga-kutta
    150145
    151          The user should supply initial conditions in y_0. For example
    152          if your initial conditions are y_0=1,y_1=1, do
     146    - ``rk8pd`` - runga-kutta prince-dormand (8,9)
    153147
    154             sage: T.y_0=[1,1]
     148    - ``rk2imp`` - implicit 2nd order runga-kutta at gaussian points
    155149
    156          The actual solver is invoked by the method ode_solve.
    157          It has arguments t_span, y_0,num_points, params.
    158          y_0 must be supplied either as an argument or above by assignment.
    159          params which is optional and only necessary if your system uses params can be supplied to ode_solve
    160          or by assignment.
    161          
    162          t_span is the time interval on which to solve the ode.
    163          If t_span is a tuple with just two time
    164          values then the user must specify num_points,
    165          and the system will be evaluated at num_points equally spaced points between t_span[0]
    166          and t_span[1]. If t_span is a tuple with more than two values than the
    167          values of y_i at points in time listed in t_span will be returned.
     150    - ``rk4imp`` - implicit 4th order runga-kutta at gaussian points
    168151
    169          Error is estimated via the expression
     152    - ``bsimp`` - implicit burlisch-stoer (requires jacobian)
    170153
    171          D_i = error_abs*s_i+error_rel*(a|y_i|+a_dydt*h*|y_i'|)
     154    - ``gear1`` - M=1 implicit gear
    172155
    173          The user can specify  error_abs (1e-10 by default),  error_rel (1e-10 by default)
    174          a (1 by default),  a_(dydt) (0 by default) and s_i (as scaling_abs which should be a
    175          tuple and is 1 in all components by default).
    176          If you specify one of a or a_dydt you must specify the other.
    177          You may specify a and a_dydt without
    178          scaling_abs (which will be taken =1 be default)
     156    - ``gear2`` - M=2 implicit gear
    179157
    180          h is the initial step size which is (1e-2) by default.
     158    The default algorithm if ``rkf45``. If you instead wanted to use
     159    ``bsimp`` you would do::
    181160
    182          ode_solve solves the solution as a list of tuples of the form,
    183          [ (t_0,[y_1,...,y_n]),(t_1,[y_1,...,y_n]),...,(t_n,[y_1,...,y_n])].
     161        sage: T.algorithm="bsimp"
    184162
    185          This data is stored in the variable solutions
    186          
    187              sage: T.solution               # not tested
     163    The user should supply initial conditions in y_0. For example if
     164    your initial conditions are y_0=1,y_1=1, do::
    188165
    189          Examples:
    190          Consider solving the Van der pol oscillator x''(t) + ux'(t)(x(t)^2-1)+x(t)=0 between t=0 and t= 100.
    191          As a first order system it is
    192          x'=y
    193          y'=-x+uy(1-x^2)
    194          Let us take u=10 and use initial conditions (x,y)=(1,0) and use the runga-kutta prince-dormand
    195          algorithm.
     166        sage: T.y_0=[1,1]
    196167
    197          sage: def f_1(t,y,params):
    198          ...      return[y[1],-y[0]-params[0]*y[1]*(y[0]**2-1.0)]
     168    The actual solver is invoked by the method :meth:`ode_solve`.  It
     169    has arguments ``t_span``, ``y_0``, ``num_points``, ``params``.
     170    ``y_0`` must be supplied either as an argument or above by
     171    assignment.  Params which are optional and only necessary if your
     172    system uses params can be supplied to ``ode_solve`` or by
     173    assignment.
    199174
    200          sage: def j_1(t,y,params):
    201          ...      return [ [0.0, 1.0],[-2.0*params[0]*y[0]*y[1]-1.0,-params[0]*(y[0]*y[0]-1.0)], [0.0,0.0] ]
     175    ``t_span`` is the time interval on which to solve the ode.  If
     176    ``t_span`` is a tuple with just two time values then the user must
     177    specify num_points, and the system will be evaluated at num_points
     178    equally spaced points between ``t_span[0]`` and ``t_span[1]``. If
     179    ``t_span`` is a tuple with more than two values than the values of
     180    y_i at points in time listed in ``t_span`` will be returned.
    202181
    203          sage: T=ode_solver()
    204          sage: T.algorithm="rk8pd"
    205          sage: T.function=f_1
    206          sage: T.jacobian=j_1
    207          sage: T.ode_solve(y_0=[1,0],t_span=[0,100],params=[10.0],num_points=1000)
    208          sage: outfile = SAGE_TMP + 'sage.png'
    209          sage: T.plot_solution(filename=outfile)
     182    Error is estimated via the expression ``D_i =
     183    error_abs*s_i+error_rel*(a|y_i|+a_dydt*h*|y_i'|)``.  The user can
     184    specify ``error_abs`` (1e-10 by default), ``error_rel`` (1e-10 by
     185    default) ``a`` (1 by default), ``a_(dydt)`` (0 by default) and
     186    ``s_i`` (as scaling_abs which should be a tuple and is 1 in all
     187    components by default).  If you specify one of ``a`` or ``a_dydt``
     188    you must specify the other.  You may specify ``a`` and ``a_dydt``
     189    without ``scaling_abs`` (which will be taken =1 be default).
     190    ``h`` is the initial step size which is (1e-2) by default.
    210191
    211          The solver line is equivalent to
    212          sage: T.ode_solve(y_0=[1,0],t_span=[x/10.0 for x in range(1000)],params = [10.0])
     192    ``ode_solve`` solves the solution as a list of tuples of the form,
     193    ``[ (t_0,[y_1,...,y_n]),(t_1,[y_1,...,y_n]),...,(t_n,[y_1,...,y_n])]``.
    213194
     195    This data is stored in the variable solutions::
     196   
     197        sage: T.solution               # not tested
    214198
    215          Lets try a system
     199    EXAMPLES::
    216200
    217          y_0'=y_1*y_2
    218          y_1'=-y_0*y_2
    219          y_2'=-.51*y_0*y_1
     201    Consider solving the Van der Pol oscillator `x''(t) +
     202    ux'(t)(x(t)^2-1)+x(t)=0` between `t=0` and `t= 100`.  As a first
     203    order system it is `x'=y`, `y'=-x+uy(1-x^2)`. Let us take `u=10`
     204    and use initial conditions `(x,y)=(1,0)` and use the runga-kutta
     205    prince-dormand algorithm.::
    220206
    221          We will not use the jacobian this time and will change the error tolerances.
     207        sage: def f_1(t,y,params):
     208        ...      return[y[1],-y[0]-params[0]*y[1]*(y[0]**2-1.0)]
    222209
    223          sage: g_1= lambda t,y: [y[1]*y[2],-y[0]*y[2],-0.51*y[0]*y[1]]
    224          sage: T.function=g_1
    225          sage: T.y_0=[0,1,1]
    226          sage: T.scale_abs=[1e-4,1e-4,1e-5]
    227          sage: T.error_rel=1e-4
    228          sage: T.ode_solve(t_span=[0,12],num_points=100)
     210        sage: def j_1(t,y,params):
     211        ...      return [ [0.0, 1.0],[-2.0*params[0]*y[0]*y[1]-1.0,-params[0]*(y[0]*y[0]-1.0)], [0.0,0.0] ]
    229212
    230          By default T.plot_solution() plots the y_0, to plot general y_i use
    231          sage: T.plot_solution(i=0, filename=outfile)
    232          sage: T.plot_solution(i=1, filename=outfile)
    233          sage: T.plot_solution(i=2, filename=outfile)
     213        sage: T=ode_solver()
     214        sage: T.algorithm="rk8pd"
     215        sage: T.function=f_1
     216        sage: T.jacobian=j_1
     217        sage: T.ode_solve(y_0=[1,0],t_span=[0,100],params=[10.0],num_points=1000)
     218        sage: outfile = SAGE_TMP + 'sage.png'
     219        sage: T.plot_solution(filename=outfile)
     220       
     221    The solver line is equivalent to::
    234222
    235          The method interpolate_solution will return a spline interpolation through the points found by the solver. By default y_0 is interpolated.
    236          You can interpolate y_i through the keyword argument i.
     223        sage: T.ode_solve(y_0=[1,0],t_span=[x/10.0 for x in range(1000)],params = [10.0])
    237224
    238          sage: f = T.interpolate_solution()
    239          sage: plot(f,0,12).show()
    240          sage: f = T.interpolate_solution(i=1)
    241          sage: plot(f,0,12).show()
    242          sage: f = T.interpolate_solution(i=2)
    243          sage: plot(f,0,12).show()
    244          sage: f = T.interpolate_solution()
    245          sage: f(pi)
    246          0.5379...
     225    Lets try a system::
    247226
    248          The solver attributes may also be set up using arguments to ode_solver.  The previous example can be rewritten as
     227        y_0'=y_1*y_2
     228        y_1'=-y_0*y_2
     229        y_2'=-.51*y_0*y_1
    249230
    250          sage: T = ode_solver(g_1,y_0=[0,1,1],scale_abs=[1e-4,1e-4,1e-5],error_rel=1e-4, algorithm="rk8pd")
    251          sage: T.ode_solve(t_span=[0,12],num_points=100)
    252          sage: f = T.interpolate_solution()
    253          sage: f(pi)
    254          0.5379...       
     231    We will not use the jacobian this time and will change the
     232    error tolerances.::
    255233
    256          Unfortunately because python functions are used, this solver is slow on system that require many function evaluations.
    257          It is possible to pass a compiled function by deriving from the class ode_sysem and overloading c_f and c_j with C functions that
    258          specify the system. The following will work in notebook
     234        sage: g_1= lambda t,y: [y[1]*y[2],-y[0]*y[2],-0.51*y[0]*y[1]]
     235        sage: T.function=g_1
     236        sage: T.y_0=[0,1,1]
     237        sage: T.scale_abs=[1e-4,1e-4,1e-5]
     238        sage: T.error_rel=1e-4
     239        sage: T.ode_solve(t_span=[0,12],num_points=100)
    259240
    260          %cython
    261          cimport sage.gsl.ode
    262          import sage.gsl.ode
    263          include 'gsl.pxi'
     241    By default T.plot_solution() plots the y_0, to plot general y_i use::
    264242
    265          cdef class van_der_pol(sage.gsl.ode.ode_system):
    266              cdef int c_f(self,double t, double *y,double *dydt):
    267                  dydt[0]=y[1]
    268                  dydt[1]=-y[0]-1000*y[1]*(y[0]*y[0]-1)
    269                  return GSL_SUCCESS
    270              cdef int c_j(self, double t,double *y,double *dfdy,double *dfdt):
    271                  dfdy[0]=0
    272                  dfdy[1]=1.0
    273                  dfdy[2]=-2.0*1000*y[0]*y[1]-1.0
    274                  dfdy[3]=-1000*(y[0]*y[0]-1.0)
    275                  dfdt[0]=0
    276                  dfdt[1]=0
    277                  return GSL_SUCCESS
     243        sage: T.plot_solution(i=0, filename=outfile)
     244        sage: T.plot_solution(i=1, filename=outfile)
     245        sage: T.plot_solution(i=2, filename=outfile)
    278246
    279          After executing the above block of code you can do the
    280          following (WARNING: the following is *not* automatically
    281          doctested):
     247    The method interpolate_solution will return a spline interpolation
     248    through the points found by the solver. By default y_0 is
     249    interpolated.  You can interpolate y_i through the keyword
     250    argument i.::
    282251
    283          sage: T = ode_solver()                     # not tested
    284          sage: T.algorithm = "bsimp"                # not tested
    285          sage: vander = van_der_pol()               # not tested
    286          sage: T.function=vander                    # not tested
    287          sage: T.ode_solve(y_0 = [1,0], t_span=[0,2000], num_points=1000)   # not tested
    288          sage: T.plot_solution(i=0, filename=SAGE_TMP + '/test.png')        # not tested
     252        sage: f = T.interpolate_solution()
     253        sage: plot(f,0,12).show()
     254        sage: f = T.interpolate_solution(i=1)
     255        sage: plot(f,0,12).show()
     256        sage: f = T.interpolate_solution(i=2)
     257        sage: plot(f,0,12).show()
     258        sage: f = T.interpolate_solution()
     259        sage: f(pi)
     260        0.5379...
    289261
    290    
    291    """
    292    def __init__(self,function=None,jacobian=None,h = 1e-2,error_abs=1e-10,error_rel=1e-10, a=False,a_dydt=False,scale_abs=False,algorithm="rkf45",y_0=None,t_span=None,params = []):
    293       self.function = function
    294       self.jacobian = jacobian
    295       self.h = h
    296       self.error_abs = error_abs
    297       self.error_rel = error_rel
    298       self.a = a
    299       self.a_dydt = a_dydt
    300       self.scale_abs = scale_abs
    301       self.algorithm = algorithm
    302       self.y_0 = y_0
    303       self.t_span = t_span
    304       self.params = params
    305       self.solution = []
     262    The solver attributes may also be set up using arguments to
     263    ode_solver.  The previous example can be rewritten as::
    306264
    307    def __setattr__(self,name,value):
    308       if(hasattr(self,'solution')):
     265        sage: T = ode_solver(g_1,y_0=[0,1,1],scale_abs=[1e-4,1e-4,1e-5],error_rel=1e-4, algorithm="rk8pd")
     266        sage: T.ode_solve(t_span=[0,12],num_points=100)
     267        sage: f = T.interpolate_solution()
     268        sage: f(pi)
     269        0.5379...
     270
     271    Unfortunately because python functions are used, this solver
     272    is slow on system that require many function evaluations.  It
     273    is possible to pass a compiled function by deriving from the
     274    class ode_sysem and overloading c_f and c_j with C functions
     275    that specify the system. The following will work in notebook::
     276
     277          %cython
     278          cimport sage.gsl.ode
     279          import sage.gsl.ode
     280          include 'gsl.pxi'
     281
     282          cdef class van_der_pol(sage.gsl.ode.ode_system):
     283              cdef int c_f(self,double t, double *y,double *dydt):
     284                  dydt[0]=y[1]
     285                  dydt[1]=-y[0]-1000*y[1]*(y[0]*y[0]-1)
     286                  return GSL_SUCCESS
     287              cdef int c_j(self, double t,double *y,double *dfdy,double *dfdt):
     288                  dfdy[0]=0
     289                  dfdy[1]=1.0
     290                  dfdy[2]=-2.0*1000*y[0]*y[1]-1.0
     291                  dfdy[3]=-1000*(y[0]*y[0]-1.0)
     292                  dfdt[0]=0
     293                  dfdt[1]=0
     294                  return GSL_SUCCESS
     295
     296    After executing the above block of code you can do the
     297    following (WARNING: the following is *not* automatically
     298    doctested)::
     299
     300        sage: T = ode_solver()                     # not tested
     301        sage: T.algorithm = "bsimp"                # not tested
     302        sage: vander = van_der_pol()               # not tested
     303        sage: T.function=vander                    # not tested
     304        sage: T.ode_solve(y_0 = [1,0], t_span=[0,2000], num_points=1000)   # not tested
     305        sage: T.plot_solution(i=0, filename=SAGE_TMP + '/test.png')        # not tested
     306
     307
     308    """
     309    def __init__(self,function=None,jacobian=None,h = 1e-2,error_abs=1e-10,error_rel=1e-10, a=False,a_dydt=False,scale_abs=False,algorithm="rkf45",y_0=None,t_span=None,params = []):
     310        self.function = function
     311        self.jacobian = jacobian
     312        self.h = h
     313        self.error_abs = error_abs
     314        self.error_rel = error_rel
     315        self.a = a
     316        self.a_dydt = a_dydt
     317        self.scale_abs = scale_abs
     318        self.algorithm = algorithm
     319        self.y_0 = y_0
     320        self.t_span = t_span
     321        self.params = params
     322        self.solution = []
     323
     324    def __setattr__(self,name,value):
     325        if(hasattr(self,'solution')):
    309326            object.__setattr__(self,'solution',[])
    310       object.__setattr__(self,name,value)
     327        object.__setattr__(self,name,value)
    311328
    312    def interpolate_solution(self,i=0):
    313       l=eval('[ (x[0],x[1][i]) for x in solution]',{'solution':self.solution,'i':i})
    314       return sage.gsl.interpolation.spline(l)
    315      
     329    def interpolate_solution(self,i=0):
     330        l=eval('[ (x[0],x[1][i]) for x in solution]',{'solution':self.solution,'i':i})
     331        return sage.gsl.interpolation.spline(l)
    316332
    317    def plot_solution(self, i=0, filename=None, interpolate=False):
    318       from sage.plot.all import plot, point
    319       points=[]
    320       for x in self.solution:
    321          points.append(point((x[0],x[1][i])))
    322       t = plot(points)
    323       if filename is None:
    324          t.show()
    325       else:
    326          t.save(filename=filename)
    327    
    328    def ode_solve(self,t_span=False,y_0=False,num_points=False,params=[]):
    329       import inspect
    330       cdef double h # step size
    331       h=self.h
    332       cdef int i
    333       cdef int j
    334       cdef int type
    335       cdef int dim
    336       cdef PyFunctionWrapper wrapper #struct to pass information into GSL C function
    337       self.params=params
    338      
    339       if t_span != False:
    340          self.t_span = t_span     
    341       if y_0 != False:
    342          self.y_0 = y_0
    343333
    344       dim = len(self.y_0)
    345       type = isinstance(self.function,ode_system)
    346       if type == 0:
    347          wrapper = PyFunctionWrapper()
    348          if self.function!=None:
    349             wrapper.the_function = self.function
    350          else:
    351             raise ValueError, "ODE system not yet defined"
    352          if self.jacobian is None:
    353             wrapper.the_jacobian = None
    354          else:   
    355             wrapper.the_jacobian = self.jacobian
    356          if self.params==[] and len(inspect.getargspec(wrapper.the_function)[0])==2:
    357             wrapper.the_parameters=[]
    358          elif self.params==[] and len(inspect.getargspec(wrapper.the_function)[0])>2:
    359             raise ValueError, "ODE system has a parameter but no parameters specified"
    360          elif self.params!=[]:
    361             wrapper.the_parameters = self.params
    362          wrapper.y_n = dim
     334    def plot_solution(self, i=0, filename=None, interpolate=False):
     335        from sage.plot.all import plot, point
     336        points=[]
     337        for x in self.solution:
     338            points.append(point((x[0],x[1][i])))
     339        t = plot(points)
     340        if filename is None:
     341            t.show()
     342        else:
     343            t.save(filename=filename)
    363344
     345    def ode_solve(self,t_span=False,y_0=False,num_points=False,params=[]):
     346        import inspect
     347        cdef double h # step size
     348        h=self.h
     349        cdef int i
     350        cdef int j
     351        cdef int type
     352        cdef int dim
     353        cdef PyFunctionWrapper wrapper #struct to pass information into GSL C function
     354        self.params=params
    364355
    365       cdef double t
    366       cdef double t_end
    367       cdef double *y
    368       cdef double * scale_abs_array
    369       scale_abs_array=NULL   
    370      
    371       y= <double*> malloc(sizeof(double)*(dim))
    372       if y==NULL:
    373          raise MemoryError,"error allocating memory"
    374       result=[]
    375       v=[0]*dim   
    376       cdef gsl_odeiv_step_type * T
     356        if t_span != False:
     357            self.t_span = t_span
     358        if y_0 != False:
     359            self.y_0 = y_0
    377360
    378       for i from 0 <=i< dim: #copy initial conditions into C array
    379           y[i]=self.y_0[i]
     361        dim = len(self.y_0)
     362        type = isinstance(self.function,ode_system)
     363        if type == 0:
     364            wrapper = PyFunctionWrapper()
     365            if self.function!=None:
     366                wrapper.the_function = self.function
     367            else:
     368                raise ValueError, "ODE system not yet defined"
     369            if self.jacobian is None:
     370                wrapper.the_jacobian = None
     371            else:
     372                wrapper.the_jacobian = self.jacobian
     373            if self.params==[] and len(inspect.getargspec(wrapper.the_function)[0])==2:
     374                wrapper.the_parameters=[]
     375            elif self.params==[] and len(inspect.getargspec(wrapper.the_function)[0])>2:
     376                raise ValueError, "ODE system has a parameter but no parameters specified"
     377            elif self.params!=[]:
     378                wrapper.the_parameters = self.params
     379            wrapper.y_n = dim
    380380
    381       if self.algorithm == "rkf45":
    382          T=gsl_odeiv_step_rkf45
    383       elif self.algorithm == "rk2":
    384          T=gsl_odeiv_step_rk2
    385       elif self.algorithm == "rk4":
    386          T=gsl_odeiv_step_rk4
    387       elif self.algorithm == "rkck":
    388          T=gsl_odeiv_step_rkck
    389       elif self.algorithm == "rk8pd":
    390          T=gsl_odeiv_step_rk8pd
    391       elif self.algorithm == "rk2imp":
    392          T= gsl_odeiv_step_rk2imp
    393       elif self.algorithm == "rk4imp":
    394          T= gsl_odeiv_step_rk4imp
    395       elif self.algorithm == "bsimp":
    396          T = gsl_odeiv_step_bsimp
    397          if not type and self.jacobian==None:
    398             raise TypeError,"The jacobian must be provided for the implicit Burlisch-Stoer method"
    399       elif self.algorithm == "gear1":
    400          T = gsl_odeiv_step_gear1
    401       elif self.algorithm == "gear2":
    402          T = gsl_odeiv_step_gear2
    403       else:
    404          raise TypeError,"algorithm not valid"
    405381
     382        cdef double t
     383        cdef double t_end
     384        cdef double *y
     385        cdef double * scale_abs_array
     386        scale_abs_array=NULL
    406387
    407       cdef gsl_odeiv_step * s 
    408       s  = gsl_odeiv_step_alloc (T, dim)
    409       if s==NULL:
    410          free(y)
    411          raise MemoryError, "error setting up solver"
     388        y= <double*> malloc(sizeof(double)*(dim))
     389        if y==NULL:
     390            raise MemoryError,"error allocating memory"
     391        result=[]
     392        v=[0]*dim
     393        cdef gsl_odeiv_step_type * T
    412394
     395        for i from 0 <=i< dim: #copy initial conditions into C array
     396            y[i]=self.y_0[i]
    413397
    414       cdef gsl_odeiv_control * c
     398        if self.algorithm == "rkf45":
     399            T=gsl_odeiv_step_rkf45
     400        elif self.algorithm == "rk2":
     401            T=gsl_odeiv_step_rk2
     402        elif self.algorithm == "rk4":
     403            T=gsl_odeiv_step_rk4
     404        elif self.algorithm == "rkck":
     405            T=gsl_odeiv_step_rkck
     406        elif self.algorithm == "rk8pd":
     407            T=gsl_odeiv_step_rk8pd
     408        elif self.algorithm == "rk2imp":
     409            T= gsl_odeiv_step_rk2imp
     410        elif self.algorithm == "rk4imp":
     411            T= gsl_odeiv_step_rk4imp
     412        elif self.algorithm == "bsimp":
     413            T = gsl_odeiv_step_bsimp
     414            if not type and self.jacobian==None:
     415                raise TypeError,"The jacobian must be provided for the implicit Burlisch-Stoer method"
     416        elif self.algorithm == "gear1":
     417            T = gsl_odeiv_step_gear1
     418        elif self.algorithm == "gear2":
     419            T = gsl_odeiv_step_gear2
     420        else:
     421            raise TypeError,"algorithm not valid"
    415422
    416       if self.a == False and self.a_dydt==False:
    417          c  = gsl_odeiv_control_y_new (self.error_abs, self.error_rel)
    418       elif self.a !=False and self.a_dydt != False:
    419          if self.scale_abs==False:
    420             c = gsl_odeiv_control_standard_new(self.error_abs,self.error_rel,self.a,self.a_dydt)
    421          elif hasattr(self.scale_abs,'__len__'):
    422             if len(self.scale_abs)==dim:
    423                scale_abs_array =<double *> malloc(dim*sizeof(double))
    424                for i from 0 <=i<dim:
    425                   scale_abs_array[i]=self.scale_abs[i]
    426                c = gsl_odeiv_control_scaled_new(self.error_abs,self.error_rel,self.a,self.a_dydt,scale_abs_array,dim)
    427423
    428       if c == NULL:
    429          gsl_odeiv_control_free (c)
    430          gsl_odeiv_step_free (s)
    431          free(y)
    432          if scale_abs_array!=NULL:
    433             free(scale_abs_array)
    434          raise MemoryError, "error setting up solver"
     424        cdef gsl_odeiv_step * s
     425        s  = gsl_odeiv_step_alloc (T, dim)
     426        if s==NULL:
     427            free(y)
     428            raise MemoryError, "error setting up solver"
    435429
    436430
    437       cdef gsl_odeiv_evolve * e
    438       e  = gsl_odeiv_evolve_alloc(dim)
     431        cdef gsl_odeiv_control * c
    439432
    440       if e == NULL:
    441          gsl_odeiv_control_free (c)
    442          gsl_odeiv_step_free (s)
    443          free(y)
    444          if scale_abs_array!=NULL:
    445             free(scale_abs_array)
    446          raise MemoryError, "error setting up solver"
     433        if self.a == False and self.a_dydt==False:
     434            c  = gsl_odeiv_control_y_new (self.error_abs, self.error_rel)
     435        elif self.a !=False and self.a_dydt != False:
     436            if self.scale_abs==False:
     437                c = gsl_odeiv_control_standard_new(self.error_abs,self.error_rel,self.a,self.a_dydt)
     438            elif hasattr(self.scale_abs,'__len__'):
     439                if len(self.scale_abs)==dim:
     440                    scale_abs_array =<double *> malloc(dim*sizeof(double))
     441                    for i from 0 <=i<dim:
     442                        scale_abs_array[i]=self.scale_abs[i]
     443                    c = gsl_odeiv_control_scaled_new(self.error_abs,self.error_rel,self.a,self.a_dydt,scale_abs_array,dim)
    447444
    448 
    449       cdef gsl_odeiv_system sys
    450       if type:               # The user has passed a class with a compiled function, use that for the system
    451          sys.function = c_f_compiled             
    452          sys.jacobian = c_jac_compiled
    453 #         (<ode_system>self.function).the_parameters = self.params
    454          sys.params = <void *> self.function   
    455       else:                  # The user passed a python function.
    456          sys.function = c_f
    457          sys.jacobian = c_jac
    458          sys.params = <void *> wrapper
    459       sys.dimension = dim
    460 
    461 
    462       cdef int status
    463       import copy
    464       cdef int n
    465 
    466       if len(self.t_span)==2 and num_points!=False:
    467          try:
    468             n = num_points
    469          except TypeError:
    470             gsl_odeiv_evolve_free (e)
     445        if c == NULL:
    471446            gsl_odeiv_control_free (c)
    472447            gsl_odeiv_step_free (s)
    473448            free(y)
    474449            if scale_abs_array!=NULL:
    475                free(scale_abs_array)
    476             raise TypeError,"numpoints must be integer"
    477          result.append( (self.t_span[0],self.y_0))
    478          delta = (self.t_span[1]-self.t_span[0])/(1.0*num_points)
    479          t =self.t_span[0]
    480          t_end=self.t_span[0]+delta
    481          for i from 0<i<=n:
    482             while (t < t_end):
    483                try:
    484                    _sig_on
    485                    status = gsl_odeiv_evolve_apply (e, c, s, &sys, &t, t_end, &h, y)
    486                    _sig_off
    487                except:
    488                   gsl_odeiv_evolve_free (e)
    489                   gsl_odeiv_control_free (c)
    490                   gsl_odeiv_step_free (s)
    491                   free(y)
    492                   if scale_abs_array!=NULL:
    493                      free(scale_abs_array)               
    494                   raise ValueError,"error solving"
     450                free(scale_abs_array)
     451            raise MemoryError, "error setting up solver"
    495452
    496                  
    497                if (status != GSL_SUCCESS):               
    498                   gsl_odeiv_evolve_free (e)
    499                   gsl_odeiv_control_free (c)
    500                   gsl_odeiv_step_free (s)
    501                   free(y)
    502                   if scale_abs_array!=NULL:
    503                      free(scale_abs_array)               
    504                   raise ValueError,"error solving"
    505453
    506             for j  from 0<=j<dim:
    507                   v[j]=<double> y[j]
    508             result.append( (t,copy.copy(v)) )
    509             t = t_end
    510             t_end= t+delta
     454        cdef gsl_odeiv_evolve * e
     455        e  = gsl_odeiv_evolve_alloc(dim)
    511456
    512       elif len(self.t_span) > 2:
    513          n = len(self.t_span)
    514          result.append((self.t_span[0],self.y_0))
    515          t=self.t_span[0]
    516          t_end=self.t_span[1]
    517          for i from 0<i<n-1:
    518             while (t < t_end):
    519                try:
    520                    _sig_on
    521                    status = gsl_odeiv_evolve_apply (e, c, s, &sys, &t, t_end, &h, y)
    522                    _sig_off
    523                except:
    524                   gsl_odeiv_evolve_free (e)
    525                   gsl_odeiv_control_free (c)
    526                   gsl_odeiv_step_free (s)
    527                   free(y)
    528                   if scale_abs_array!=NULL:
    529                      free(scale_abs_array)               
    530                   raise ValueError,"error solving"
     457        if e == NULL:
     458            gsl_odeiv_control_free (c)
     459            gsl_odeiv_step_free (s)
     460            free(y)
     461            if scale_abs_array!=NULL:
     462                free(scale_abs_array)
     463            raise MemoryError, "error setting up solver"
    531464
    532                  
    533                if (status != GSL_SUCCESS):
    534                   gsl_odeiv_evolve_free (e)
    535                   gsl_odeiv_control_free (c)
    536                   gsl_odeiv_step_free (s)
    537                   free(y)
    538                   if scale_abs_array!=NULL:
    539                      free(scale_abs_array)               
    540                   raise ValueError,"error solving"
    541465
     466        cdef gsl_odeiv_system sys
     467        if type:               # The user has passed a class with a compiled function, use that for the system
     468            sys.function = c_f_compiled
     469            sys.jacobian = c_jac_compiled
     470#         (<ode_system>self.function).the_parameters = self.params
     471            sys.params = <void *> self.function
     472        else:                  # The user passed a python function.
     473            sys.function = c_f
     474            sys.jacobian = c_jac
     475            sys.params = <void *> wrapper
     476        sys.dimension = dim
    542477
    543                
    544             for j from 0<=j<dim:
    545                v[j]=<double> y[j]
    546             result.append( (t,copy.copy(v)) )
    547478
    548             t=t_span[i]
    549             t_end=t_span[i+1]
     479        cdef int status
     480        import copy
     481        cdef int n
    550482
     483        if len(self.t_span)==2 and num_points!=False:
     484            try:
     485                n = num_points
     486            except TypeError:
     487                gsl_odeiv_evolve_free (e)
     488                gsl_odeiv_control_free (c)
     489                gsl_odeiv_step_free (s)
     490                free(y)
     491                if scale_abs_array!=NULL:
     492                    free(scale_abs_array)
     493                raise TypeError,"numpoints must be integer"
     494            result.append( (self.t_span[0],self.y_0))
     495            delta = (self.t_span[1]-self.t_span[0])/(1.0*num_points)
     496            t =self.t_span[0]
     497            t_end=self.t_span[0]+delta
     498            for i from 0<i<=n:
     499                while (t < t_end):
     500                    try:
     501                        _sig_on
     502                        status = gsl_odeiv_evolve_apply (e, c, s, &sys, &t, t_end, &h, y)
     503                        _sig_off
     504                    except:
     505                        gsl_odeiv_evolve_free (e)
     506                        gsl_odeiv_control_free (c)
     507                        gsl_odeiv_step_free (s)
     508                        free(y)
     509                        if scale_abs_array!=NULL:
     510                            free(scale_abs_array)
     511                        raise ValueError,"error solving"
    551512
    552       gsl_odeiv_evolve_free (e)
    553       gsl_odeiv_control_free (c)
    554       gsl_odeiv_step_free (s)
    555       free(y)
    556       if scale_abs_array!=NULL:
    557          free(scale_abs_array)
    558       self.solution = result
    559513
     514                    if (status != GSL_SUCCESS):
     515                        gsl_odeiv_evolve_free (e)
     516                        gsl_odeiv_control_free (c)
     517                        gsl_odeiv_step_free (s)
     518                        free(y)
     519                        if scale_abs_array!=NULL:
     520                            free(scale_abs_array)
     521                        raise ValueError,"error solving"
    560522
     523                for j  from 0<=j<dim:
     524                    v[j]=<double> y[j]
     525                result.append( (t,copy.copy(v)) )
     526                t = t_end
     527                t_end= t+delta
    561528
     529        elif len(self.t_span) > 2:
     530            n = len(self.t_span)
     531            result.append((self.t_span[0],self.y_0))
     532            t=self.t_span[0]
     533            t_end=self.t_span[1]
     534            for i from 0<i<n-1:
     535                while (t < t_end):
     536                    try:
     537                        _sig_on
     538                        status = gsl_odeiv_evolve_apply (e, c, s, &sys, &t, t_end, &h, y)
     539                        _sig_off
     540                    except:
     541                        gsl_odeiv_evolve_free (e)
     542                        gsl_odeiv_control_free (c)
     543                        gsl_odeiv_step_free (s)
     544                        free(y)
     545                        if scale_abs_array!=NULL:
     546                            free(scale_abs_array)
     547                        raise ValueError,"error solving"
    562548
     549
     550                    if (status != GSL_SUCCESS):
     551                        gsl_odeiv_evolve_free (e)
     552                        gsl_odeiv_control_free (c)
     553                        gsl_odeiv_step_free (s)
     554                        free(y)
     555                        if scale_abs_array!=NULL:
     556                            free(scale_abs_array)
     557                        raise ValueError,"error solving"
     558
     559
     560
     561                for j from 0<=j<dim:
     562                    v[j]=<double> y[j]
     563                result.append( (t,copy.copy(v)) )
     564
     565                t=t_span[i]
     566                t_end=t_span[i+1]
     567
     568
     569        gsl_odeiv_evolve_free (e)
     570        gsl_odeiv_control_free (c)
     571        gsl_odeiv_step_free (s)
     572        free(y)
     573        if scale_abs_array!=NULL:
     574            free(scale_abs_array)
     575        self.solution = result
  • sage/misc/functional.py

    diff -r b89cf68ccc1c -r ea16778cc7ae sage/misc/functional.py
    a b  
    663663
    664664def integral(x, *args, **kwds):
    665665    """
    666     Returns an indefinite integral of an object x.
     666    Returns an indefinite or definite integral of an object x.
    667667   
    668668    First call x.integrate() and if that fails make an object and
    669669    integrate it using Maxima, maple, etc, as specified by algorithm.
     670
     671    For symbolic expression calls
     672    ``sage.calculus.calculus.integral`` - see this function for
     673    available options.
    670674   
    671675    EXAMPLES::
    672676   
    673677        sage: f = cyclotomic_polynomial(10)
    674678        sage: integral(f)
    675679        1/5*x^5 - 1/4*x^4 + 1/3*x^3 - 1/2*x^2 + x
     680
     681    ::
     682
    676683        sage: integral(sin(x),x)
    677684        -cos(x)
     685
     686    ::
     687
    678688        sage: y = var('y')
    679689        sage: integral(sin(x),y)
    680690        y*sin(x)
     691
     692    ::
     693
    681694        sage: integral(sin(x), x, 0, pi/2)
    682695        1
    683696        sage: sin(x).integral(x, 0,pi/2)