Ticket #10682: trac_10682-reviewer.patch

File trac_10682-reviewer.patch, 55.8 KB (added by Jean-Pierre Flori, 11 years ago)

Reviewer patch, formatting issues

  • sage/calculus/calculus.py

    # HG changeset patch
    # User Jean-Pierre Flori <jean-pierre.flor@ssi.gouv.fr>
    # Date 1331030079 -3600
    # Node ID 4ebd0aae52ae9b89870a211578c815127fd98b5b
    # Parent  0a84cd9c2c97b45d3bf8e2e8a10d337f7139b98f
    #10682: Reviewer patch
    
    diff --git a/sage/calculus/calculus.py b/sage/calculus/calculus.py
    a b  
    247247    sage: CC(f)
    248248    1.12762596520638 + 1.17520119364380*I
    249249    sage: ComplexField(200)(f)
    250     1.1276259652063807852262251614026720125478471180986674836290 + 1.1752011936438014568823818505956008151557179813340958702296*I
     250    1.1276259652063807852262251614026720125478471180986674836290
     251    + 1.1752011936438014568823818505956008151557179813340958702296*I
    251252    sage: ComplexField(100)(f)
    252253    1.1276259652063807852262251614 + 1.1752011936438014568823818506*I
    253254
     
    256257
    257258    sage: f = sum(1/var('n%s'%i)^i for i in range(10))
    258259    sage: f
    259     1/n1 + 1/n2^2 + 1/n3^3 + 1/n4^4 + 1/n5^5 + 1/n6^6 + 1/n7^7 + 1/n8^8 + 1/n9^9 + 1
     260    1/n1 + 1/n2^2 + 1/n3^3 + 1/n4^4 + 1/n5^5 + 1/n6^6 + 1/n7^7 + 1/n8^8 + 1/n9^9
     261    + 1
    260262
    261263Note that after calling var, the variables are immediately
    262264available for use::
     
    267269We can, of course, substitute::
    268270
    269271    sage: f(n9=9,n7=n6)
    270     1/n1 + 1/n2^2 + 1/n3^3 + 1/n4^4 + 1/n5^5 + 1/n6^6 + 1/n6^7 + 1/n8^8 + 387420490/387420489
     272    1/n1 + 1/n2^2 + 1/n3^3 + 1/n4^4 + 1/n5^5 + 1/n6^6 + 1/n6^7 + 1/n8^8
     273    + 387420490/387420489
    271274
    272275TESTS:
    273276
     
    632635      - ``5`` - integral is probably divergent or slowly
    633636        convergent
    634637
    635       - ``6`` - the input is invalid; this includes the case of desired_relative_error
    636         being too small to be achieved
     638      - ``6`` - the input is invalid; this includes the case of
     639                desired_relative_error being too small to be achieved
    637640
    638641    ALIAS: nintegrate is the same as nintegral
    639642
     
    716719        else:
    717720            raise TypeError, err
    718721
    719     #This is just a work around until there is a response to
    720     #http://www.math.utexas.edu/pipermail/maxima/2008/012975.html
     722    # Maxima returns an unevaluated expression when the underlying SLATEC
     723    # library fails. See:
     724    # http://www.math.utexas.edu/pipermail/maxima/2008/012975.html
     725    # and Trac ticket #10682 for more information.
    721726    if 'quad_qags' in str(v):
    722727        raise ValueError, "Maxima (via quadpack) cannot compute the integral"
    723728
     
    806811        sage: sin(pi/7).minpoly()
    807812        x^6 - 7/4*x^4 + 7/8*x^2 - 7/64
    808813        sage: minpoly(exp(I*pi/17))
    809         x^16 - x^15 + x^14 - x^13 + x^12 - x^11 + x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1
     814        x^16 - x^15 + x^14 - x^13 + x^12 - x^11 + x^10 - x^9 + x^8 - x^7 + x^6
     815        - x^5 + x^4 - x^3 + x^2 - x + 1
    810816
    811817    Here we verify it gives the same result as the abstract number
    812818    field.
     
    884890    algorithm='numerical'::
    885891
    886892        sage: cos(pi/33).minpoly(algorithm='algebraic')
    887         x^10 + 1/2*x^9 - 5/2*x^8 - 5/4*x^7 + 17/8*x^6 + 17/16*x^5 - 43/64*x^4 - 43/128*x^3 + 3/64*x^2 + 3/128*x + 1/1024
     893        x^10 + 1/2*x^9 - 5/2*x^8 - 5/4*x^7 + 17/8*x^6 + 17/16*x^5 - 43/64*x^4
     894        - 43/128*x^3 + 3/64*x^2 + 3/128*x + 1/1024
    888895        sage: cos(pi/33).minpoly(algorithm='numerical')
    889896        Traceback (most recent call last):
    890897        ...
     
    13061313
    13071314    .. math::
    13081315
    1309                       F(s) = \frac{1}{2\pi i} \int_{\gamma-i\infty}^{\gamma + i\infty} e^{st} F(s) dt,
     1316        F(s) = \frac{1}{2\pi i} \int_{\gamma-i\infty}^{\gamma + i\infty} e^{st} F(s) dt,
    13101317
    13111318    where `\gamma` is chosen so that the contour path of
    13121319    integration is in the region of convergence of `F(s)`.
     
    13271334        sage: inverse_laplace(L, s, t)
    13281335        t |--> t*cos(t)
    13291336        sage: inverse_laplace(1/(s^3+1), s, t)
    1330         1/3*(sqrt(3)*sin(1/2*sqrt(3)*t) - cos(1/2*sqrt(3)*t))*e^(1/2*t) + 1/3*e^(-t)
     1337        1/3*(sqrt(3)*sin(1/2*sqrt(3)*t) - cos(1/2*sqrt(3)*t))*e^(1/2*t)
     1338        + 1/3*e^(-t)
    13311339
    13321340    No explicit inverse Laplace transform, so one is returned formally
    13331341    as a function ``ilt``::
     
    13691377        sage: diff(u(x+h), x)
    13701378        D[0](u)(h + x)
    13711379        sage: taylor(u(x+h),h,0,4)
    1372         1/24*h^4*D[0, 0, 0, 0](u)(x) + 1/6*h^3*D[0, 0, 0](u)(x) + 1/2*h^2*D[0, 0](u)(x) + h*D[0](u)(x) + u(x)
     1380        1/24*h^4*D[0, 0, 0, 0](u)(x) + 1/6*h^3*D[0, 0, 0](u)(x)
     1381        + 1/2*h^2*D[0, 0](u)(x) + h*D[0](u)(x) + u(x)
    13731382
    13741383    We compute a Laplace transform::
    13751384
  • sage/calculus/desolvers.py

    diff --git a/sage/calculus/desolvers.py b/sage/calculus/desolvers.py
    a b  
    3535- ``eulers_method_2x2_plot`` - Plots the sequence of points obtained
    3636  from Euler's method.
    3737
    38 AUTHORS: 
     38AUTHORS:
    3939
    4040- David Joyner (3-2006) - Initial version of functions
    4141
     
    6868def desolve(de, dvar, ics=None, ivar=None, show_method=False, contrib_ode=False):
    6969    r"""
    7070    Solves a 1st or 2nd order linear ODE via maxima. Including IVP and BVP.
    71    
     71
    7272    *Use* ``desolve? <tab>`` *if the output in truncated in notebook.*
    7373
    7474    INPUT:
    75    
     75
    7676    - ``de`` - an expression or equation representing the ODE
    77    
     77
    7878    - ``dvar`` - the dependent variable (hereafter called ``y``)
    7979
    8080    - ``ics`` - (optional) the initial or boundary conditions
    8181
    8282      - for a first-order equation, specify the initial ``x`` and ``y``
    83      
     83
    8484      - for a second-order equation, specify the initial ``x``, ``y``,
    8585        and ``dy/dx``, i.e. write `[x_0, y(x_0), y'(x_0)]`
    86      
     86
    8787      - for a second-order boundary solution, specify initial and
    88         final ``x`` and ``y`` boundary conditions, i.e. write `[x_0, y(x_0), x_1, y(x_1)]`.
    89        
    90       - gives an error if the solution is not SymbolicEquation (as happens for
     88        final ``x`` and ``y`` boundary conditions, i.e. write
     89        `[x_0, y(x_0), x_1, y(x_1)]`.
     90
     91      - gives an error if the solution is not SymbolicEquation (as happens for
    9192        example for Clairaut equation)
    92        
     93
    9394    - ``ivar`` - (optional) the independent variable (hereafter called
    9495      x), which must be specified if there is more than one
    9596      independent variable in the equation.
    96                    
     97
    9798    - ``show_method`` - (optional) if true, then Sage returns pair
    9899      ``[solution, method]``, where method is the string describing
    99100      method which has been used to get solution (Maxima uses the
     
    111112      contain singular solution, for example)
    112113
    113114    OUTPUT:
    114    
     115
    115116    In most cases returns SymbolicEquation which defines the solution
    116117    implicitly.  If the result is in the form y(x)=... (happens for
    117118    linear eqs.), returns the right-hand side only.  The possible
     
    124125        sage: y = function('y', x)
    125126        sage: desolve(diff(y,x) + y - 1, y)
    126127        (c + e^x)*e^(-x)
    127    
     128
    128129    ::
    129130
    130131        sage: f = desolve(diff(y,x) + y - 1, y, ics=[10,2]); f
     
    133134    ::
    134135
    135136        sage: plot(f)
    136        
     137
    137138    We can also solve second-order differential equations.::
    138    
     139
    139140        sage: x = var('x')
    140141        sage: y = function('y', x)
    141142        sage: de = diff(y,x,2) - y == x
    142143        sage: desolve(de, y)
    143144        k1*e^x + k2*e^(-x) - x
    144        
     145
    145146
    146147    ::
    147148
     
    166167
    167168    ::
    168169
    169         sage: desolve(de, y, [0,1,pi/2,4]) 
     170        sage: desolve(de, y, [0,1,pi/2,4])
    170171        4*sin(x) + cos(x)
    171172
    172173    ::
     
    178179
    179180        sage: desolve(diff(y,x)^2+x*diff(y,x)-y==0,y,contrib_ode=True,show_method=True)
    180181        [[y(x) == c^2 + c*x, y(x) == -1/4*x^2], 'clairault']
    181    
     182
    182183    For equations involving more variables we specify independent variable::
    183184
    184185        sage: a,b,c,n=var('a b c n')
     
    248249    condition at x=0, since this point is singlar point of the
    249250    equation. Anyway, if the solution should be bounded at x=0, then
    250251    k2=0.::
    251            
     252
    252253        sage: desolve(x^2*diff(y,x,x)+x*diff(y,x)+(x^2-4)*y==0,y)
    253254        k1*bessel_j(2, x) + k2*bessel_y(2, x)
    254    
     255
    255256    Difficult ODE produces error::
    256257
    257258        sage: desolve(sqrt(y)*diff(y,x)+e^(y)+cos(x)-sin(x+y)==0,y) # not tested
    258259        Traceback (click to the left for traceback)
    259260        ...
    260261        NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True."
    261        
     262
    262263    Difficult ODE produces error - moreover, takes a long time ::
    263264
    264265        sage: desolve(sqrt(y)*diff(y,x)+e^(y)+cos(x)-sin(x+y)==0,y,contrib_ode=True) # not tested
     
    269270        [[y(x) == c + log(x), y(x) == c*e^x], 'factor']
    270271
    271272    ::
    272    
     273
    273274        sage: desolve(diff(y,x)==(x+y)^2,y,contrib_ode=True,show_method=True)
    274275        [[[x == c - arctan(sqrt(t)), y(x) == -x - sqrt(t)], [x == c + arctan(sqrt(t)), y(x) == -x + sqrt(t)]], 'lagrange']
    275276
     
    281282        Traceback (click to the left for traceback)
    282283        ...
    283284        NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True."
    284        
     285
    285286    ::
    286287
    287288        sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,2],show_method=True) # not tested
     
    293294
    294295        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y)
    295296        (k2*x + k1)*e^(-x) + 1/2*sin(x)
    296        
     297
    297298    ::
    298        
     299
    299300        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,show_method=True)
    300301        [(k2*x + k1)*e^(-x) + 1/2*sin(x), 'variationofparameters']
    301        
     302
    302303    ::
    303        
     304
    304305        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,[0,3,1])
    305306        1/2*(7*x + 6)*e^(-x) + 1/2*sin(x)
    306        
     307
    307308    ::
    308        
     309
    309310        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,[0,3,1],show_method=True)
    310311        [1/2*(7*x + 6)*e^(-x) + 1/2*sin(x), 'variationofparameters']
    311        
     312
    312313    ::
    313        
     314
    314315        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,[0,3,pi/2,2])
    315316        3*((e^(1/2*pi) - 2)*x/pi + 1)*e^(-x) + 1/2*sin(x)
    316        
     317
    317318    ::
    318        
     319
    319320        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,[0,3,pi/2,2],show_method=True)
    320321        [3*((e^(1/2*pi) - 2)*x/pi + 1)*e^(-x) + 1/2*sin(x), 'variationofparameters']
    321        
     322
    322323    ::
    323        
     324
    324325        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y)
    325326        (k2*x + k1)*e^(-x)
    326        
     327
    327328    ::
    328        
     329
    329330        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,show_method=True)
    330331        [(k2*x + k1)*e^(-x), 'constcoeff']
    331        
     332
    332333    ::
    333        
     334
    334335        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,[0,3,1])
    335336        (4*x + 3)*e^(-x)
    336        
     337
    337338    ::
    338        
     339
    339340        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,[0,3,1],show_method=True)
    340341        [(4*x + 3)*e^(-x), 'constcoeff']
    341        
     342
    342343    ::
    343        
     344
    344345        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,[0,3,pi/2,2])
    345346        (2*(2*e^(1/2*pi) - 3)*x/pi + 3)*e^(-x)
    346        
     347
    347348    ::
    348        
     349
    349350        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,[0,3,pi/2,2],show_method=True)
    350351        [(2*(2*e^(1/2*pi) - 3)*x/pi + 3)*e^(-x), 'constcoeff']
    351        
     352
    352353    TESTS:
    353    
     354
    354355    Trac #9961 fixed (allow assumptions on the dependent variable in desolve)::
    355    
     356
    356357        sage: y=function('y',x); assume(x>0); assume(y>0)
    357358        sage: sage.calculus.calculus.maxima('domain:real')  # needed since Maxima 5.26.0 to get the answer as below
    358359        real
    359360        sage: desolve(x*diff(y,x)-x*sqrt(y^2+x^2)-y == 0, y, contrib_ode=True)
    360361        [x - arcsinh(y(x)/x) == c]
    361362
    362     Trac #10682 updated Maxima to 5.26, and it started to show a different solution in the complex domain for the ODE above::
     363    Trac #10682 updated Maxima to 5.26, and it started to show a different
     364    solution in the complex domain for the ODE above::
    363365
    364         sage: sage.calculus.calculus.maxima('domain:complex')  # back to the dafault, complex, domain
     366        sage: sage.calculus.calculus.maxima('domain:complex')  # back to the default, complex, domain
    365367        complex
    366368        sage: desolve(x*diff(y,x)-x*sqrt(y^2+x^2)-y == 0, y, contrib_ode=True)
    367         [1/2*(2*x^2*sqrt(x^(-2)) - 2*x*sqrt(x^(-2))*arcsinh(y(x)/sqrt(x^2))
    368         - 2*x*sqrt(x^(-2))*arcsinh(y(x)^2/(sqrt(y(x)^2)*x))
    369         + log(4*(2*x^2*sqrt((x^2*y(x)^2 + y(x)^4)/x^2)*sqrt(x^(-2)) + x^2 + 2*y(x)^2)/x^2))/(x*sqrt(x^(-2))) == c]
     369        [1/2*(2*x^2*sqrt(x^(-2)) - 2*x*sqrt(x^(-2))*arcsinh(y(x)/sqrt(x^2))
     370         - 2*x*sqrt(x^(-2))*arcsinh(y(x)^2/(sqrt(y(x)^2)*x))
     371         + log(4*(2*x^2*sqrt((x^2*y(x)^2 + y(x)^4)/x^2)*sqrt(x^(-2))
     372         + x^2 + 2*y(x)^2)/x^2))/(x*sqrt(x^(-2))) == c]
    370373
    371374    Trac #6479 fixed::
    372375
     
    376379        x
    377380
    378381    ::
    379    
     382
    380383        sage: desolve( diff(y,x,x) == 0, y, [0,1,1])
    381384        x + 1
    382385
     
    394397        sage: x=var('x'); f=function('f',x); k=var('k'); assume(k>0)
    395398        sage: desolve(diff(f,x,2)/f==k,f,ivar=x)
    396399        k1*e^(sqrt(k)*x) + k2*e^(-sqrt(k)*x)
    397    
    398        
     400
     401
    399402    AUTHORS:
    400  
     403
    401404    - David Joyner (1-2006)
    402405
    403406    - Robert Bradshaw (10-2008)
     
    419422            raise ValueError, "Unable to determine independent variable, please specify."
    420423        ivar = ivars[0]
    421424    def sanitize_var(exprs):
    422         return exprs.replace("'"+dvar_str+"("+ivar_str+")",dvar_str)   
     425        return exprs.replace("'"+dvar_str+"("+ivar_str+")",dvar_str)
    423426    de00 = de._maxima_()
    424427    P = de00.parent()
    425428    dvar_str=P(dvar.operator()).str()
     
    427430    de00 = de00.str()
    428431    de0 = sanitize_var(de00)
    429432    ode_solver="ode2"
    430     cmd="(TEMP:%s(%s,%s,%s), if TEMP=false then TEMP else substitute(%s=%s(%s),TEMP))"%(ode_solver,de0,dvar_str,ivar_str,dvar_str,dvar_str,ivar_str)   
     433    cmd="(TEMP:%s(%s,%s,%s), if TEMP=false then TEMP else substitute(%s=%s(%s),TEMP))"%(ode_solver,de0,dvar_str,ivar_str,dvar_str,dvar_str,ivar_str)
    431434    # we produce string like this
    432435    # ode2('diff(y,x,2)+2*'diff(y,x,1)+y-cos(x),y(x),x)
    433436    soln = P(cmd)
     
    436439        if contrib_ode:
    437440            ode_solver="contrib_ode"
    438441            P("load('contrib_ode)")
    439             cmd="(TEMP:%s(%s,%s,%s), if TEMP=false then TEMP else substitute(%s=%s(%s),TEMP))"%(ode_solver,de0,dvar_str,ivar_str,dvar_str,dvar_str,ivar_str)   
     442            cmd="(TEMP:%s(%s,%s,%s), if TEMP=false then TEMP else substitute(%s=%s(%s),TEMP))"%(ode_solver,de0,dvar_str,ivar_str,dvar_str,dvar_str,ivar_str)
    440443            # we produce string like this
    441444            # (TEMP:contrib_ode(x*('diff(y,x,1))^2-(x*y+1)*'diff(y,x,1)+y,y,x), if TEMP=false then TEMP else substitute(y=y(x),TEMP))
    442445            soln = P(cmd)
    443446            if str(soln).strip() == 'false':
    444                 raise NotImplementedError, "Maxima was unable to solve this ODE."   
     447                raise NotImplementedError, "Maxima was unable to solve this ODE."
    445448        else:
    446449            raise NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True."
    447450
     
    462465            # (TEMP:ic2(ode2('diff(y,x,2)+2*'diff(y,x,1)+y-cos(x),y,x),x=0,y=3,'diff(y,x)=1),substitute(y=y(x),TEMP))
    463466            soln=P(cmd)
    464467        if len(ics) == 3:
    465             #fixed ic2 command from Maxima - we have to ensure that %k1, %k2 do not depend on variables, should be removed when fixed in Maxima 
     468            #fixed ic2 command from Maxima - we have to ensure that %k1, %k2 do not depend on variables, should be removed when fixed in Maxima
    466469            P("ic2_sage(soln,xa,ya,dya):=block([programmode:true,backsubst:true,singsolve:true,temp,%k2,%k1,TEMP_k], \
    467470                noteqn(xa), noteqn(ya), noteqn(dya), boundtest('%k1,%k1), boundtest('%k2,%k2), \
    468471                temp: lhs(soln) - rhs(soln), \
     
    479482            # (TEMP:ic2(ode2('diff(y,x,2)+2*'diff(y,x,1)+y-cos(x),y,x),x=0,y=3,'diff(y,x)=1),substitute(y=y(x),TEMP))
    480483            soln=P(cmd)
    481484            if str(soln).strip() == 'false':
    482                 raise NotImplementedError, "Maxima was unable to solve this IVP. Remove the initial condition to get the general solution."   
     485                raise NotImplementedError, "Maxima was unable to solve this IVP. Remove the initial condition to get the general solution."
    483486        if len(ics) == 4:
    484             #fixed bc2 command from Maxima - we have to ensure that %k1, %k2 do not depend on variables, should be removed when fixed in Maxima 
     487            #fixed bc2 command from Maxima - we have to ensure that %k1, %k2 do not depend on variables, should be removed when fixed in Maxima
    485488            P("bc2_sage(soln,xa,ya,xb,yb):=block([programmode:true,backsubst:true,singsolve:true,temp,%k1,%k2,TEMP_k], \
    486489                noteqn(xa), noteqn(ya), noteqn(xb), noteqn(yb), boundtest('%k1,%k1), boundtest('%k2,%k2), \
    487490                TEMP_k:solve([subst([xa,ya],soln), subst([xb,yb],soln)], [%k1,%k2]), \
     
    489492                temp: maplist(lambda([zz], subst(zz,soln)),TEMP_k), \
    490493                if length(temp)=1 then return(first(temp)) else return(temp))")
    491494            cmd="bc2_sage(%s(%s,%s,%s),%s,%s=%s,%s,%s=%s)"%(ode_solver,de00,dvar_str,ivar_str,P(ivar==ics[0]).str(),dvar_str,P(ics[1]).str(),P(ivar==ics[2]).str(),dvar_str,P(ics[3]).str())
    492             cmd="(TEMP:%s,substitute(%s=%s(%s),TEMP))"%(cmd,dvar_str,dvar_str,ivar_str)           
     495            cmd="(TEMP:%s,substitute(%s=%s(%s),TEMP))"%(cmd,dvar_str,dvar_str,ivar_str)
    493496            cmd=sanitize_var(cmd)
    494497            # we produce string like this
    495498            # (TEMP:bc2(ode2('diff(y,x,2)+2*'diff(y,x,1)+y-cos(x),y,x),x=0,y=3,x=%pi/2,y=2),substitute(y=y(x),TEMP))
    496499            soln=P(cmd)
    497500            if str(soln).strip() == 'false':
    498                 raise NotImplementedError, "Maxima was unable to solve this BVP. Remove the initial condition to get the general solution."   
     501                raise NotImplementedError, "Maxima was unable to solve this BVP. Remove the initial condition to get the general solution."
    499502
    500503    soln=soln.sage()
    501504    if is_SymbolicEquation(soln) and soln.lhs() == dvar:
     
    504507        soln = soln.rhs()
    505508    if show_method:
    506509        return [soln,maxima_method.str()]
    507     else:   
     510    else:
    508511        return soln
    509512
    510513
     
    549552#    #return maxima.eval(cmd)
    550553#    return de0.desolve(vars[0]).rhs()
    551554
    552    
     555
    553556def desolve_laplace(de, dvar, ics=None, ivar=None):
    554557    """
    555558    Solves an ODE using laplace transforms. Initials conditions are optional.
    556559
    557560    INPUT:
    558    
     561
    559562    - ``de`` - a lambda expression representing the ODE (eg, de =
    560563      diff(y,x,2) == diff(y,x)+sin(x))
    561564
     
    573576    Solution of the ODE as symbolic expression
    574577
    575578    EXAMPLES::
    576    
     579
    577580        sage: u=function('u',x)
    578581        sage: eq = diff(u,x) - exp(-x) - u == 0
    579582        sage: desolve_laplace(eq,u)
    580583        1/2*(2*u(0) + 1)*e^x - 1/2*e^(-x)
    581    
     584
    582585    We can use initial conditions::
    583586
    584587        sage: desolve_laplace(eq,u,ics=[0,3])
    585588        -1/2*e^(-x) + 7/2*e^x
    586589
    587     The initial conditions do not persist in the system (as they persisted 
     590    The initial conditions do not persist in the system (as they persisted
    588591    in previous versions)::
    589592
    590593        sage: desolve_laplace(eq,u)
     
    596599        sage: eq = diff(f,x) + f == 0
    597600        sage: desolve_laplace(eq,f,[0,1])
    598601        e^(-x)
    599    
     602
    600603    ::
    601    
     604
    602605        sage: x = var('x')
    603606        sage: f = function('f', x)
    604607        sage: de = diff(f,x,x) - 2*diff(f,x) + f
     
    624627
    625628        sage: soln(t=3)
    626629        e^(-3) + 1
    627    
    628     AUTHORS: 
     630
     631    AUTHORS:
    629632
    630633    - David Joyner (1-2006,8-2007)
    631    
     634
    632635    - Robert Marik (10-2009)
    633636    """
    634637    #This is the original code from David Joyner (inputs and outputs strings)
     
    659662    ## verbatim copy from desolve - end
    660663
    661664    def sanitize_var(exprs):  # 'y(x) -> y(x)
    662         return exprs.replace("'"+str(dvar),str(dvar))   
     665        return exprs.replace("'"+str(dvar),str(dvar))
    663666    de0=de._maxima_()
    664667    P = de0.parent()
    665668    cmd = sanitize_var("desolve("+de0.str()+","+str(dvar)+")")
    666669    soln=P(cmd).rhs()
    667670    if str(soln).strip() == 'false':
    668         raise NotImplementedError, "Maxima was unable to solve this ODE."   
     671        raise NotImplementedError, "Maxima was unable to solve this ODE."
    669672    soln=soln.sage()
    670673    if ics!=None:
    671674        d = len(ics)
    672675        for i in range(0,d-1):
    673             soln=eval('soln.substitute(diff(dvar,ivar,i)('+str(ivar)+'=ics[0])==ics[i+1])') 
     676            soln=eval('soln.substitute(diff(dvar,ivar,i)('+str(ivar)+'=ics[0])==ics[i+1])')
    674677    return soln
    675    
    676    
     678
     679
    677680def desolve_system(des, vars, ics=None, ivar=None):
    678681    """
    679682    Solves any size system of 1st order ODE's. Initials conditions are optional.
    680    
    681     Onedimensional systems are passed to :meth:`desolve_laplace`. 
     683
     684    Onedimensional systems are passed to :meth:`desolve_laplace`.
    682685
    683686    INPUT:
    684    
     687
    685688    - ``des`` - list of ODEs
    686    
     689
    687690    - ``vars`` - list of dependent variables
    688    
     691
    689692    - ``ics`` - (optional) list of initial values for ivar and vars
    690    
     693
    691694    - ``ivar`` - (optional) the independent variable, which must be
    692695      specified if there is more than one independent variable in the
    693696      equation.
     
    702705        sage: desolve_system([de1, de2], [x,y])
    703706        [x(t) == (x(0) - 1)*cos(t) - (y(0) - 1)*sin(t) + 1,
    704707         y(t) == (x(0) - 1)*sin(t) + (y(0) - 1)*cos(t) + 1]
    705          
     708
    706709    Now we give some initial conditions::
    707    
     710
    708711        sage: sol = desolve_system([de1, de2], [x,y], ics=[0,1,2]); sol
    709712        [x(t) == -sin(t) + 1, y(t) == cos(t) + 1]
    710    
     713
    711714    ::
    712    
     715
    713716        sage: solnx, solny = sol[0].rhs(), sol[1].rhs()
    714717        sage: plot([solnx,solny],(0,1))  # not tested
    715718        sage: parametric_plot((solnx,solny),(0,1))  # not tested
     
    717720    TESTS:
    718721
    719722    Trac #9823 fixed::
    720    
     723
    721724        sage: t = var('t')
    722725        sage: x = function('x', t)
    723726        sage: de1 = diff(x,t) + 1 == 0
    724727        sage: desolve_system([de1], [x])
    725728        -t + x(0)
    726729
    727     AUTHORS: 
     730    AUTHORS:
    728731
    729732    - Robert Bradshaw (10-2008)
    730733    """
     
    756759        for dvar, ic in zip(dvars, ics[:1]):
    757760            dvar.atvalue(ivar==ivar_ic, dvar)
    758761    return soln
    759        
     762
    760763
    761764def desolve_system_strings(des,vars,ics=None):
    762765    r"""
     
    788791        sage: s = var('s')
    789792        sage: function('x', s)
    790793        x(s)
    791    
     794
    792795    ::
    793    
     796
    794797        sage: function('y', s)
    795798        y(s)
    796    
     799
    797800    ::
    798    
     801
    799802        sage: de1 = lambda z: diff(z[0],s) + z[1] - 1
    800803        sage: de2 = lambda z: diff(z[1],s) - z[0] + 1
    801804        sage: des = [de1([x(s),y(s)]),de2([x(s),y(s)])]
    802805        sage: vars = ["s","x","y"]
    803806        sage: desolve_system_strings(des,vars)
    804         ["(1-'y(0))*sin(s)+('x(0)-1)*cos(s)+1", "('x(0)-1)*sin(s)+('y(0)-1)*cos(s)+1"]
    805    
     807        ["(1-'y(0))*sin(s)+('x(0)-1)*cos(s)+1",
     808         "('x(0)-1)*sin(s)+('y(0)-1)*cos(s)+1"]
     809
    806810    ::
    807    
     811
    808812        sage: ics = [0,1,-1]
    809813        sage: soln = desolve_system_strings(des,vars,ics); soln
    810814        ['2*sin(s)+1', '1-2*cos(s)']
    811    
     815
    812816    ::
    813    
     817
    814818        sage: solnx, solny = map(SR, soln)
    815819        sage: RR(solnx(s=3))
    816820        1.28224001611973
    817    
     821
    818822    ::
    819    
     823
    820824        sage: P1 = plot([solnx,solny],(0,1))
    821825        sage: P2 = parametric_plot((solnx,solny),(0,1))
    822826
    823827    Now type show(P1), show(P2) to view these.
    824                      
    825828
    826     AUTHORS:
     829
     830    AUTHORS:
    827831
    828832    - David Joyner (3-2006, 8-2007)
    829833    """
     
    856860    last column.
    857861
    858862    *For pedagogical purposes only.*
    859    
     863
    860864    EXAMPLES::
    861    
     865
    862866        sage: from sage.calculus.desolvers import eulers_method
    863867        sage: x,y = PolynomialRing(QQ,2,"xy").gens()
    864868        sage: eulers_method(5*x+y-5,0,1,1/2,1)
     
    912916        sage: P2 = line(pts)
    913917        sage: (P1+P2).show()
    914918
    915     AUTHORS: 
     919    AUTHORS:
    916920
    917921    - David Joyner
    918922    """
     
    948952    next (last) column.
    949953
    950954    *For pedagogical purposes only.*
    951    
     955
    952956    EXAMPLES::
    953957
    954958        sage: from sage.calculus.desolvers import eulers_method_2x2
    955959        sage: t, x, y = PolynomialRing(QQ,3,"txy").gens()
    956960        sage: f = x+y+t; g = x-y
    957961        sage: eulers_method_2x2(f,g, 0, 0, 0, 1/3, 1,algorithm="none")
    958         [[0, 0, 0], [1/3, 0, 0], [2/3, 1/9, 0], [1, 10/27, 1/27], [4/3, 68/81, 4/27]]
     962        [[0, 0, 0], [1/3, 0, 0], [2/3, 1/9, 0], [1, 10/27, 1/27],
     963         [4/3, 68/81, 4/27]]
    959964
    960965    ::
    961966
     
    994999           3/4                 0.63                   -0.0078               -0.031                 0.11
    9951000             1                 0.63                     0.020                0.079                0.071
    9961001
    997     To numerically approximate y(1), where `y''+ty'+y=0`, `y(0)=1`, `y'(0)=0`:: 
     1002    To numerically approximate y(1), where `y''+ty'+y=0`, `y(0)=1`, `y'(0)=0`::
    9981003
    9991004        sage: t,x,y=PolynomialRing(RR,3,"txy").gens()
    10001005        sage: f = y; g = -x-y*t
     
    10061011           3/4                 0.88                     -0.15                -0.62                -0.10
    10071012             1                 0.75                     -0.17                -0.68               -0.015
    10081013
    1009     AUTHORS: 
     1014    AUTHORS:
    10101015
    10111016    - David Joyner
    10121017    """
     
    10361041    `x(a)=x_0`, `y' = g(t,x,y)`, `y(a) = y_0`.
    10371042
    10381043    *For pedagogical purposes only.*
    1039    
     1044
    10401045    EXAMPLES::
    10411046
    10421047        sage: from sage.calculus.desolvers import eulers_method_2x2_plot
     
    10661071def desolve_rk4_determine_bounds(ics,end_points=None):
    10671072    """
    10681073    Used to determine bounds for numerical integration.
    1069    
     1074
    10701075    - If end_points is None, the interval for integration is from ics[0]
    10711076      to ics[0]+10
    10721077
    1073     - If end_points is a or [a], the interval for integration is from min(ics[0],a)
    1074       to max(ics[0],a)
     1078    - If end_points is a or [a], the interval for integration is from
     1079      min(ics[0],a) to max(ics[0],a)
    10751080
    10761081    - If end_points is [a,b], the interval for integration is from min(ics[0],a)
    10771082      to max(ics[0],b)
     
    10811086        sage: from sage.calculus.desolvers import desolve_rk4_determine_bounds
    10821087        sage: desolve_rk4_determine_bounds([0,2],1)
    10831088        (0, 1)
    1084    
     1089
    10851090    ::
    10861091
    1087         sage: desolve_rk4_determine_bounds([0,2]) 
     1092        sage: desolve_rk4_determine_bounds([0,2])
    10881093        (0, 10)
    1089    
     1094
    10901095    ::
    10911096
    10921097        sage: desolve_rk4_determine_bounds([0,2],[-2])
    10931098        (-2, 0)
    1094    
     1099
    10951100    ::
    10961101
    10971102        sage: desolve_rk4_determine_bounds([0,2],[-2,4])
    10981103        (-2, 4)
    1099    
     1104
    11001105    """
    1101     if end_points is None: 
     1106    if end_points is None:
    11021107        return((ics[0],ics[0]+10))
    11031108    if not isinstance(end_points,list):
    11041109        end_points=[end_points]
     
    11061111        return (min(ics[0],end_points[0]),max(ics[0],end_points[0]))
    11071112    else:
    11081113        return (min(ics[0],end_points[0]),max(ics[0],end_points[1]))
    1109    
     1114
    11101115
    11111116def desolve_rk4(de, dvar, ics=None, ivar=None, end_points=None, step=0.1, output='list', **kwds):
    11121117    """
     
    11141119    equation. See also ``ode_solver``.
    11151120
    11161121    INPUT:
    1117    
    1118     input is similar to ``desolve`` command. The differential equation can be 
     1122
     1123    input is similar to ``desolve`` command. The differential equation can be
    11191124    written in a form close to the plot_slope_field or desolve command
    11201125
    11211126    - Variant 1 (function in two variables)
    1122      
    1123       - ``de`` - right hand side, i.e. the function `f(x,y)` from ODE `y'=f(x,y)`
    1124      
     1127
     1128      - ``de`` - right hand side, i.e. the function `f(x,y)` from ODE
     1129        `y'=f(x,y)`
     1130
    11251131      - ``dvar`` - dependent variable (symbolic variable declared by var)
    11261132
    11271133    - Variant 2 (symbolic equation)
    11281134
    11291135      - ``de`` - equation, including term with ``diff(y,x)``
    11301136
    1131       - ``dvar``` - dependent variable (declared as funciton of independent variable)
     1137      - ``dvar``` - dependent variable (declared as a function of independent
     1138                    variables)
    11321139
    1133     - Other parameters 
     1140    - Other parameters
    11341141
    1135       - ``ivar`` - should be specified, if there are more variables or if the equation is autonomous
     1142      - ``ivar`` - should be specified, if there are more variables or if the
     1143                   equation is autonomous
    11361144
    11371145      - ``ics`` - initial conditions in the form [x0,y0]
    1138      
     1146
    11391147      - ``end_points`` - the end points of the interval
    11401148
    1141         - if end_points is a or [a], we integrate on between min(ics[0],a) and max(ics[0],a)
     1149        - if end_points is a or [a], we integrate on between min(ics[0],a) and
     1150          max(ics[0],a)
     1151
    11421152        - if end_points is None, we use end_points=ics[0]+10
    11431153
    1144         - if end_points is [a,b] we integrate on between min(ics[0],a) and max(ics[0],b)
     1154        - if end_points is [a,b] we integrate on between min(ics[0],a) and
     1155          max(ics[0],b)
    11451156
    1146       - ``step`` - (optional, default:0.1) the length of the step (positive number)
    1147  
     1157      - ``step`` - (optional, default:0.1) the length of the step (positive
     1158                   number)
     1159
    11481160      - ``output`` - (optional, default: 'list') one of 'list',
    11491161        'plot', 'slope_field' (graph of the solution with slope field)
    11501162
    1151     OUTPUT: 
     1163    OUTPUT:
    11521164
    11531165    Returns a list of points, or plot produced by list_plot,
    11541166    optionally with slope field.
     
    11661178
    11671179    Variant 1 for input - we can pass ODE in the form used by
    11681180    desolve function In this example we integrate bakwards, since
    1169     ``end_points < ics[0]``:: 
     1181    ``end_points < ics[0]``::
    11701182
    1171         sage: y=function('y',x) 
    1172         sage: desolve_rk4(diff(y,x)+y*(y-1) == x-2,y,ics=[1,1],step=0.5, end_points=0) 
     1183        sage: y=function('y',x)
     1184        sage: desolve_rk4(diff(y,x)+y*(y-1) == x-2,y,ics=[1,1],step=0.5, end_points=0)
    11731185        [[0.0, 8.90425710896], [0.5, 1.90932794536], [1, 1]]
    11741186
    11751187    Here we show how to plot simple pictures. For more advanced
     
    11781190
    11791191        sage: x,y=var('x y')
    11801192        sage: P=desolve_rk4(y*(2-y),y,ics=[0,.1],ivar=x,output='slope_field',end_points=[-4,6],thickness=3)
    1181                
     1193
    11821194    ALGORITHM:
    11831195
    11841196    4th order Runge-Kutta method. Wrapper for command ``rk`` in
    11851197    Maxima's dynamics package.  Perhaps could be faster by using
    11861198    fast_float instead.
    11871199
    1188     AUTHORS: 
     1200    AUTHORS:
    11891201
    11901202    - Robert Marik (10-2009)
    11911203    """
    1192     if ics is None: 
     1204    if ics is None:
    11931205        raise ValueError, "No initial conditions, specify with ics=[x0,y0]."
    11941206
    11951207    if ivar is None:
     
    12461258        YMIN=sol[0][1]
    12471259        XMAX=XMIN
    12481260        YMAX=YMIN
    1249         for s,t in sol: 
     1261        for s,t in sol:
    12501262            if s>XMAX:XMAX=s
    12511263            if s<XMIN:XMIN=s
    12521264            if t>YMAX:YMAX=t
     
    12601272    r"""
    12611273    Solves numerically system of first-order ordinary differential
    12621274    equations using the 4th order Runge-Kutta method. Wrapper for
    1263     Maxima command ``rk``. See also ``ode_solver``. 
     1275    Maxima command ``rk``. See also ``ode_solver``.
    12641276
    12651277    INPUT:
    1266    
     1278
    12671279    input is similar to desolve_system and desolve_rk4 commands
    1268    
     1280
    12691281    - ``des`` - right hand sides of the system
    12701282
    12711283    - ``vars`` - dependent variables
     
    12751287      missing
    12761288
    12771289    - ``ics`` - initial conditions in the form [x0,y01,y02,y03,....]
    1278    
     1290
    12791291    - ``end_points`` - the end points of the interval
    1280    
    1281       - if end_points is a or [a], we integrate on between min(ics[0],a) and max(ics[0],a)
     1292
     1293      - if end_points is a or [a], we integrate on between min(ics[0],a) and
     1294        max(ics[0],a)
     1295
    12821296      - if end_points is None, we use end_points=ics[0]+10
    12831297
    1284       - if end_points is [a,b] we integrate on between min(ics[0],a) and max(ics[0],b)
     1298      - if end_points is [a,b] we integrate on between min(ics[0],a) and
     1299        max(ics[0],b)
    12851300
    12861301    - ``step`` -- (optional, default: 0.1) the length of the step
    12871302
     
    13001315        sage: P=desolve_system_rk4([x*(1-y),-y*(1-x)],[x,y],ics=[0,0.5,2],ivar=t,end_points=20)
    13011316        sage: Q=[ [i,j] for i,j,k in P]
    13021317        sage: LP=list_plot(Q)
    1303        
     1318
    13041319        sage: Q=[ [j,k] for i,j,k in P]
    13051320        sage: LP=list_plot(Q)
    1306    
     1321
    13071322    ALGORITHM:
    13081323
    13091324    4th order Runge-Kutta method. Wrapper for command ``rk`` in Maxima's
    13101325    dynamics package.  Perhaps could be faster by using ``fast_float``
    13111326    instead.
    13121327
    1313     AUTHOR: 
     1328    AUTHOR:
    13141329
    13151330    - Robert Marik (10-2009)
    13161331    """
    13171332
    1318     if ics is None: 
     1333    if ics is None:
    13191334        raise ValueError, "No initial conditions, specify with ics=[x0,y01,y02,...]."
    13201335
    13211336    ivars = set([])
     
    13351350    x0=ics[0]
    13361351    icss = [ics[i]._maxima_().str() for i in range(1,len(ics))]
    13371352    icstr = "[" + ",".join(icss) + "]"
    1338     step=abs(step) 
     1353    step=abs(step)
    13391354
    13401355    maxima("load('dynamics)")
    13411356    lower_bound,upper_bound=desolve_rk4_determine_bounds(ics,end_points)
     
    13611376, rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0
    13621377, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0):
    13631378    r"""
    1364     Solves numerically a system of first-order ordinary differential equations 
     1379    Solves numerically a system of first-order ordinary differential equations
    13651380    using ``odeint`` from scipy.integrate module.
    13661381
    13671382    INPUT:
    1368    
     1383
    13691384    - ``des``  -- right hand sides of the system
    13701385
    1371     - ``ics``  -- initial conditions 
     1386    - ``ics``  -- initial conditions
    13721387
    13731388    - ``times`` -- a sequence of time points in which the solution must be found
    13741389
     
    13771392
    13781393    - ``ivar`` -- independent variable, optional.
    13791394
    1380     - ``compute_jac`` -- boolean. If True, the Jacobian of des is computed and 
     1395    - ``compute_jac`` -- boolean. If True, the Jacobian of des is computed and
    13811396      used during the integration of Stiff Systems. Default value is False.
    13821397
    1383     Other Parameters (taken from the documentation of odeint function from 
     1398    Other Parameters (taken from the documentation of odeint function from
    13841399      scipy.integrate module)
    13851400
    13861401    - ``rtol``, ``atol`` : float
     
    14101425    - ``hmin`` : float, (0: solver-determined)
    14111426      The minimum absolute step size allowed.
    14121427
    1413     - ``ixpr`` : boolean. 
     1428    - ``ixpr`` : boolean.
    14141429      Whether to generate extra printing at method switches.
    14151430
    14161431    - ``mxstep`` : integer, (0: solver-determined)
     
    14251440
    14261441    - ``mxords`` : integer, (0: solver-determined)
    14271442      Maximum order to be allowed for the stiff (BDF) method.
    1428    
     1443
    14291444    OUTPUT:
    1430    
     1445
    14311446    Returns a list with the solution of the system at each time in times.
    1432    
     1447
    14331448    EXAMPLES:
    1434    
     1449
    14351450    Lotka Volterra Equations::
    14361451
    14371452        sage: from sage.calculus.desolvers import desolve_odeint
     
    14401455        sage: sol=desolve_odeint(f,[0.5,2],srange(0,10,0.1),[x,y])
    14411456        sage: p=line(zip(sol[:,0],sol[:,1]))
    14421457        sage: p.show()
    1443    
     1458
    14441459    Lorenz Equations::
    14451460
    14461461        sage: x,y,z=var('x,y,z')
     
    14491464        sage: rho=28
    14501465        sage: beta=8/3
    14511466        sage: # The Lorenz equations
    1452         sage: lorenz=[sigma*(y-x),x*(rho-z)-y,x*y-beta*z] 
     1467        sage: lorenz=[sigma*(y-x),x*(rho-z)-y,x*y-beta*z]
    14531468        sage: # Time and initial conditions
    14541469        sage: times=srange(0,50.05,0.05)
    14551470        sage: ics=[0,1,1]
    14561471        sage: sol=desolve_odeint(lorenz,ics,times,[x,y,z],rtol=1e-13,atol=1e-14)
    1457    
     1472
    14581473    One-dimensional Stiff system::
    14591474
    14601475        sage: y= var('y')
     
    14651480        sage: sol=desolve_odeint(f,ic,t,y,rtol=1e-9,atol=1e-10,compute_jac=True)
    14661481        sage: p=points(zip(t,sol))
    14671482        sage: p.show()
    1468    
    1469     Another Stiff system with some optional parameters with no 
     1483
     1484    Another Stiff system with some optional parameters with no
    14701485    default value::
    14711486
    14721487        sage: y1,y2,y3=var('y1,y2,y3')
     
    14781493        sage: t=srange(0,10,0.01)
    14791494        sage: v=[y1,y2,y3]
    14801495        sage: sol=desolve_odeint(f,ci,t,v,rtol=1e-3,atol=1e-4,h0=0.1,hmax=1,hmin=1e-4,mxstep=1000,mxords=17)
    1481    
     1496
    14821497    AUTHOR:
    1483    
     1498
    14841499    - Oriol Castejon (05-2010)
    14851500    """
    14861501
     
    15471562            def Dfun(y,t):
    15481563                v = list(y[:])
    15491564                v.append(t)
    1550                 return [[element(*v) for element in row] for row in J]     
    1551        
     1565                return [[element(*v) for element in row] for row in J]
     1566
    15521567    sol=odeint(func, ics, times, args=args, Dfun=Dfun, rtol=rtol, atol=atol,
    15531568        tcrit=tcrit, h0=h0, hmax=hmax, hmin=hmin, ixpr=ixpr, mxstep=mxstep,
    15541569        mxhnil=mxhnil, mxordn=mxordn, mxords=mxords, printmessg=printmessg)
  • sage/symbolic/relation.py

    diff --git a/sage/symbolic/relation.py b/sage/symbolic/relation.py
    a b  
    6262    0 == -10*a + b - 144
    6363
    6464Subtract two symbolic equations::
    65                
     65
    6666    sage: var('a,b')
    67     (a, b)           
     67    (a, b)
    6868    sage: m = 144 == 20 * a + b
    6969    sage: n = 136 == 10 * a + b
    7070    sage: m - n
     
    188188    sage: assume(x>0, y < 2)
    189189    sage: assumptions()
    190190    [x > 0, y < 2]
    191     sage: (y < 2).forget() 
     191    sage: (y < 2).forget()
    192192    sage: assumptions()
    193193    [x > 0]
    194194    sage: forget()
     
    222222    sage: eq = (x == 2)
    223223    sage: eq._maple_init_()
    224224    'x = 2'
    225    
     225
    226226Comparison::
    227227
    228228    sage: x = var('x')
     
    262262
    263263    sage: f = x^5 + a
    264264    sage: solve(f==0,x)
    265     [x == (-a)^(1/5)*e^(2/5*I*pi), x == (-a)^(1/5)*e^(4/5*I*pi), x == (-a)^(1/5)*e^(-4/5*I*pi), x == (-a)^(1/5)*e^(-2/5*I*pi), x == (-a)^(1/5)]
     265    [x == (-a)^(1/5)*e^(2/5*I*pi), x == (-a)^(1/5)*e^(4/5*I*pi),
     266     x == (-a)^(1/5)*e^(-4/5*I*pi), x == (-a)^(1/5)*e^(-2/5*I*pi),
     267     x == (-a)^(1/5)]
    266268
    267269You can also do arithmetic with inequalities, as illustrated
    268270below::
     
    297299    sage: eqn = x^3 + 2/3 >= x
    298300    sage: loads(dumps(eqn))
    299301    x^3 + 2/3 >= x
    300     sage: loads(dumps(eqn)) == eqn 
     302    sage: loads(dumps(eqn)) == eqn
    301303    True
    302304
    303305AUTHORS:
     
    433435    Used internally by the symbolic solve command to convert the output
    434436    of Maxima's solve command to a list of solutions in Sage's symbolic
    435437    package.
    436    
     438
    437439    EXAMPLES:
    438440
    439441    We derive the (monic) quadratic formula::
    440    
     442
    441443        sage: var('x,a,b')
    442444        (x, a, b)
    443445        sage: solve(x^2 + a*x + b == 0, x)
    444446        [x == -1/2*a - 1/2*sqrt(a^2 - 4*b), x == -1/2*a + 1/2*sqrt(a^2 - 4*b)]
    445    
     447
    446448    Behind the scenes when the above is evaluated the function
    447449    :func:`string_to_list_of_solutions` is called with input the
    448450    string `s` below::
    449    
     451
    450452        sage: s = '[x=-(sqrt(a^2-4*b)+a)/2,x=(sqrt(a^2-4*b)-a)/2]'
    451453        sage: sage.symbolic.relation.string_to_list_of_solutions(s)
    452454         [x == -1/2*a - 1/2*sqrt(a^2 - 4*b), x == -1/2*a + 1/2*sqrt(a^2 - 4*b)]
     
    464466def solve(f, *args, **kwds):
    465467    r"""
    466468    Algebraically solve an equation or system of equations (over the
    467     complex numbers) for given variables. Inequalities and systems 
    468     of inequalities are also supported.   
    469    
     469    complex numbers) for given variables. Inequalities and systems
     470    of inequalities are also supported.
     471
    470472    INPUT:
    471    
     473
    472474    -  ``f`` - equation or system of equations (given by a
    473475       list or tuple)
    474    
     476
    475477    -  ``*args`` - variables to solve for.
    476    
    477     -  ``solution_dict`` - bool (default: False); if True or non-zero, 
     478
     479    -  ``solution_dict`` - bool (default: False); if True or non-zero,
    478480       return a list of dictionaries containing the solutions. If there
    479481       are no solutions, return an empty list (rather than a list containing
    480482       an empty dictionary). Likewise, if there's only a single solution,
    481483       return a list containing one dictionary with that solution.
    482    
     484
    483485    EXAMPLES::
    484    
     486
    485487        sage: x, y = var('x, y')
    486488        sage: solve([x+y==6, x-y==4], x, y)
    487489        [[x == 5, y == 1]]
    488490        sage: solve([x^2+y^2 == 1, y^2 == x^3 + x + 1], x, y)
    489         [[x == -1/2*I*sqrt(3) - 1/2, y == -sqrt(-1/2*I*sqrt(3) + 3/2)], 
    490          [x == -1/2*I*sqrt(3) - 1/2, y == sqrt(-1/2*I*sqrt(3) + 3/2)], 
    491          [x == 1/2*I*sqrt(3) - 1/2, y == -sqrt(1/2*I*sqrt(3) + 3/2)], 
    492          [x == 1/2*I*sqrt(3) - 1/2, y == sqrt(1/2*I*sqrt(3) + 3/2)], 
    493          [x == 0, y == -1], 
     491        [[x == -1/2*I*sqrt(3) - 1/2, y == -sqrt(-1/2*I*sqrt(3) + 3/2)],
     492         [x == -1/2*I*sqrt(3) - 1/2, y == sqrt(-1/2*I*sqrt(3) + 3/2)],
     493         [x == 1/2*I*sqrt(3) - 1/2, y == -sqrt(1/2*I*sqrt(3) + 3/2)],
     494         [x == 1/2*I*sqrt(3) - 1/2, y == sqrt(1/2*I*sqrt(3) + 3/2)],
     495         [x == 0, y == -1],
    494496         [x == 0, y == 1]]
    495497        sage: solve([sqrt(x) + sqrt(y) == 5, x + y == 10], x, y)
    496         [[x == -5/2*I*sqrt(5) + 5, y == 5/2*I*sqrt(5) + 5], [x == 5/2*I*sqrt(5) + 5, y == -5/2*I*sqrt(5) + 5]]
     498        [[x == -5/2*I*sqrt(5) + 5, y == 5/2*I*sqrt(5) + 5],
     499         [x == 5/2*I*sqrt(5) + 5, y == -5/2*I*sqrt(5) + 5]]
    497500        sage: solutions=solve([x^2+y^2 == 1, y^2 == x^3 + x + 1], x, y, solution_dict=True)
    498501        sage: for solution in solutions: print solution[x].n(digits=3), ",", solution[y].n(digits=3)
    499502        -0.500 - 0.866*I , -1.27 + 0.341*I
     
    503506        0.000 , -1.00
    504507        0.000 , 1.00
    505508
    506     Whenever possible, answers will be symbolic, but with systems of 
    507     equations, at times approximations will be given, due to the 
     509    Whenever possible, answers will be symbolic, but with systems of
     510    equations, at times approximations will be given, due to the
    508511    underlying algorithm in Maxima::
    509512
    510513        sage: sols = solve([x^3==y,y^2==x],[x,y]); sols[-1], sols[0]
    511         ([x == 0, y == 0], [x == (0.309016994375 + 0.951056516295*I),  y == (-0.809016994375 - 0.587785252292*I)])
     514        ([x == 0, y == 0],
     515         [x == (0.309016994375 + 0.951056516295*I),
     516          y == (-0.809016994375 - 0.587785252292*I)])
    512517        sage: sols[0][0].rhs().pyobject().parent()
    513518        Complex Double Field
    514519
     
    516521    for symbolic expressions, which defaults to exact answers only::
    517522
    518523        sage: solve([y^6==y],y)
    519         [y == e^(2/5*I*pi), y == e^(4/5*I*pi), y == e^(-4/5*I*pi), y == e^(-2/5*I*pi), y == 1, y == 0]
     524        [y == e^(2/5*I*pi), y == e^(4/5*I*pi), y == e^(-4/5*I*pi),
     525         y == e^(-2/5*I*pi), y == 1, y == 0]
    520526        sage: solve( [y^6 == y], y)==solve( y^6 == y, y)
    521527        True
    522528
     
    540546        TypeError: 5 is not a valid variable.
    541547
    542548    If we ask for dictionaries containing the solutions, we get them::
    543    
     549
    544550        sage: solve([x^2-1],x,solution_dict=True)
    545551        [{x: -1}, {x: 1}]
    546552        sage: solve([x^2-4*x+4],x,solution_dict=True)
     
    550556        x: 2, y: 4
    551557        x: -2, y: 4
    552558
    553     If there is a parameter in the answer, that will show up as
    554     a new variable.  In the following example, ``r1`` is a real free
    555     variable (because of the ``r``)::
     559    If there is a parameter in the answer, that will show up as a new variable.
     560    In the following example, ``r1`` is a real free variable
     561    (because of the ``r``)::
    556562
    557563        sage: solve([x+y == 3, 2*x+2*y == 6],x,y)
    558564        [[x == -r1 + 3, y == r1]]
     
    574580    solutions are returned. E.g., note that the first
    575581    ``3==3`` evaluates to ``True``, not to a
    576582    symbolic equation.
    577    
     583
    578584    ::
    579    
     585
    580586        sage: solve([3==3, 1.00000000000000*x^3 == 0], x)
    581587        [x == 0]
    582588        sage: solve([1.00000000000000*x^3 == 0], x)
    583589        [x == 0]
    584    
     590
    585591    Here, the first equation evaluates to ``False``, so
    586592    there are no solutions::
    587    
     593
    588594        sage: solve([1==3, 1.00000000000000*x^3 == 0], x)
    589595        []
    590    
     596
    591597    Completely symbolic solutions are supported::
    592    
     598
    593599        sage: var('s,j,b,m,g')
    594600        (s, j, b, m, g)
    595601        sage: sys = [ m*(1-s) - b*s*j, b*s*j-g*j ];
     
    608614    Use use_grobner if no solution is obtained from to_poly_solve::
    609615
    610616       sage: x,y=var('x y'); c1(x,y)=(x-5)^2+y^2-16; c2(x,y)=(y-3)^2+x^2-9
    611        sage: solve([c1(x,y),c2(x,y)],[x,y])                               
    612        [[x == -9/68*sqrt(55) + 135/68, y == -15/68*sqrt(5)*sqrt(11) + 123/68], [x == 9/68*sqrt(55) + 135/68, y == 15/68*sqrt(5)*sqrt(11) + 123/68]]
    613        
     617       sage: solve([c1(x,y),c2(x,y)],[x,y])
     618       [[x == -9/68*sqrt(55) + 135/68, y == -15/68*sqrt(5)*sqrt(11) + 123/68],
     619        [x == 9/68*sqrt(55) + 135/68, y == 15/68*sqrt(5)*sqrt(11) + 123/68]]
     620
    614621    TESTS::
    615622
    616623        sage: solve([sin(x)==x,y^2==x],x,y)
     
    672679        variables = args
    673680    else:
    674681        variables = tuple(args[0])
    675    
     682
    676683    for v in variables:
    677684        if not is_SymbolicVariable(v):
    678685            raise TypeError, "%s is not a valid variable."%v
     
    732739    Return all solutions to an equation or list of equations modulo the
    733740    given integer modulus. Each equation must involve only polynomials
    734741    in 1 or many variables.
    735    
     742
    736743    By default the solutions are returned as `n`-tuples, where `n`
    737744    is the number of variables appearing anywhere in the given
    738745    equations. The variables are in alphabetical order.
    739    
     746
    740747    INPUT:
    741    
    742    
     748
     749
    743750    -  ``eqns`` - equation or list of equations
    744    
     751
    745752    -  ``modulus`` - an integer
    746    
     753
    747754    -  ``solution_dict`` - bool (default: False); if True or non-zero,
    748755       return a list of dictionaries containing the solutions. If there
    749756       are no solutions, return an empty list (rather than a list containing
    750757       an empty dictionary). Likewise, if there's only a single solution,
    751758       return a list containing one dictionary with that solution.
    752759
    753    
     760
    754761    EXAMPLES::
    755    
     762
    756763        sage: var('x,y')
    757764        (x, y)
    758765        sage: solve_mod([x^2 + 2 == x, x^2 + y == y^2], 14)
    759766        [(4, 2), (4, 6), (4, 9), (4, 13)]
    760767        sage: solve_mod([x^2 == 1, 4*x  == 11], 15)
    761768        [(14,)]
    762    
     769
    763770    Fermat's equation modulo 3 with exponent 5::
    764    
     771
    765772        sage: var('x,y,z')
    766773        (x, y, z)
    767774        sage: solve_mod([x^5 + y^5 == z^5], 3)
    768         [(0, 0, 0), (0, 1, 1), (0, 2, 2), (1, 0, 1), (1, 1, 2), (1, 2, 0), (2, 0, 2), (2, 1, 0), (2, 2, 1)]
     775        [(0, 0, 0), (0, 1, 1), (0, 2, 2), (1, 0, 1), (1, 1, 2), (1, 2, 0),
     776         (2, 0, 2), (2, 1, 0), (2, 2, 1)]
    769777
    770     We can solve with respect to a bigger modulus if it consists only of small prime factors::
     778    We can solve with respect to a bigger modulus if it consists only of small
     779    prime factors::
    771780
    772781        sage: [d] = solve_mod([5*x + y == 3, 2*x - 3*y == 9], 3*5*7*11*19*23*29, solution_dict = True)
    773782        sage: d[x]
     
    780789    is large::
    781790
    782791        sage: sorted(solve_mod([x^2 == 41], 10^20))
    783         [(4538602480526452429,), (11445932736758703821,), (38554067263241296179,),
    784         (45461397519473547571,), (54538602480526452429,), (61445932736758703821,),
    785         (88554067263241296179,), (95461397519473547571,)]
     792        [(4538602480526452429,), (11445932736758703821,),
     793         (38554067263241296179,), (45461397519473547571,),
     794         (54538602480526452429,), (61445932736758703821,),
     795         (88554067263241296179,), (95461397519473547571,)]
    786796
    787797    We solve a simple equation modulo 2::
    788    
     798
    789799        sage: x,y = var('x,y')
    790800        sage: solve_mod([x == y], 2)
    791801        [(0, 0), (1, 1)]
     
    821831        [(8, 5, 6)]
    822832
    823833    Confirm that modulus 1 now behaves as it should::
    824    
     834
    825835        sage: x, y = var("x y")
    826836        sage: solve_mod([x==1], 1)
    827837        [(0,)]
    828838        sage: solve_mod([2*x^2+x*y, -x*y+2*y^2+x-2*y, -2*x^2+2*x*y-y^2-x-y], 1)
    829839        [(0, 0)]
    830    
     840
    831841
    832842    """
    833843    from sage.rings.all import Integer, Integers, crt_basis
     
    844854        raise ValueError, "the modulus must be a positive integer"
    845855    vars = list(set(sum([list(e.variables()) for e in eqns], [])))
    846856    vars.sort(key=repr)
    847    
     857
    848858    if modulus == 1: # degenerate case
    849859        ans = [tuple(Integers(1)(0) for v in vars)]
    850860        return ans
     
    852862    factors = modulus.factor()
    853863    crt_basis = vector(Integers(modulus), crt_basis([p**i for p,i in factors]))
    854864    solutions = []
    855    
     865
    856866    has_solution = True
    857867    for p,i in factors:
    858868        solution =_solve_mod_prime_power(eqns, p, i, vars)
     
    861871        else:
    862872            has_solution = False
    863873            break
    864            
     874
    865875
    866876    ans = []
    867877    if has_solution:
     
    882892    Internal help function for solve_mod, does little checking since it expects
    883893    solve_mod to do that
    884894
    885     Return all solutions to an equation or list of equations modulo p^m. 
     895    Return all solutions to an equation or list of equations modulo p^m.
    886896    Each equation must involve only polynomials
    887897    in 1 or many variables.
    888898
    889899    The solutions are returned as `n`-tuples, where `n`
    890900    is the number of variables in vars.
    891        
     901
    892902    INPUT:
    893    
    894    
     903
     904
    895905    -  ``eqns`` - equation or list of equations
    896    
     906
    897907    -  ``p`` - a prime
    898908
    899909    -  ``i`` - an integer > 0
    900910
    901911    -  ``vars`` - a list of variables to solve for
    902    
    903    
     912
     913
    904914    EXAMPLES::
    905    
     915
    906916        sage: var('x,y')
    907917        (x, y)
    908918        sage: solve_mod([x^2 + 2 == x, x^2 + y == y^2], 14)
    909919        [(4, 2), (4, 6), (4, 9), (4, 13)]
    910920        sage: solve_mod([x^2 == 1, 4*x  == 11], 15)
    911921        [(14,)]
    912    
     922
    913923    Fermat's equation modulo 3 with exponent 5::
    914    
     924
    915925        sage: var('x,y,z')
    916926        (x, y, z)
    917927        sage: solve_mod([x^5 + y^5 == z^5], 3)
    918         [(0, 0, 0), (0, 1, 1), (0, 2, 2), (1, 0, 1), (1, 1, 2), (1, 2, 0), (2, 0, 2), (2, 1, 0), (2, 2, 1)]
     928        [(0, 0, 0), (0, 1, 1), (0, 2, 2), (1, 0, 1), (1, 1, 2), (1, 2, 0),
     929         (2, 0, 2), (2, 1, 0), (2, 2, 1)]
    919930
    920931    We solve a simple equation modulo 2::
    921    
     932
    922933        sage: x,y = var('x,y')
    923934        sage: solve_mod([x == y], 2)
    924935        [(0, 0), (1, 1)]
    925936
    926        
     937
    927938    .. warning::
    928    
     939
    929940       Currently this constructs possible solutions by building up
    930941       from the smallest prime factor of the modulus.  The interface
    931942       is good, but the algorithm is horrible if the modulus isn't the
     
    937948       interface.
    938949
    939950    TESTS:
    940    
     951
    941952    Confirm we can reproduce the first few terms of OEIS A187719::
    942    
     953
    943954        sage: from sage.symbolic.relation import _solve_mod_prime_power
    944955        sage: [sorted(_solve_mod_prime_power([x^2==41], 10, i, [x]))[0][0] for i in [1..13]]
    945         [1, 21, 71, 1179, 2429, 47571, 1296179, 8703821, 26452429, 526452429, 
    946         13241296179, 19473547571, 2263241296179]
     956        [1, 21, 71, 1179, 2429, 47571, 1296179, 8703821, 26452429, 526452429,
     957         13241296179, 19473547571, 2263241296179]
    947958
    948959    """
    949960    from sage.rings.all import Integers, PolynomialRing
    950961    from sage.modules.all import vector
    951962    from sage.misc.all import cartesian_product_iterator
    952    
     963
    953964    mrunning = 1
    954965    ans = []
    955966    for mi in xrange(m):
     
    970981
    971982def solve_ineq_univar(ineq):
    972983    """
    973     Function solves rational inequality in one variable. 
    974    
     984    Function solves rational inequality in one variable.
     985
    975986    INPUT:
    976    
     987
    977988    - ``ineq`` - inequality in one variable
    978    
     989
    979990    OUTPUT:
    980    
     991
    981992    - ``list`` -- output is list of solutions as a list of simple inequalities
    982       output [A,B,C] means (A or B or C) each A, B, C is again a list and 
    983       if A=[a,b], then A means (a and b). The list is empty if there is no 
     993      output [A,B,C] means (A or B or C) each A, B, C is again a list and
     994      if A=[a,b], then A means (a and b). The list is empty if there is no
    984995      solution.
    985        
     996
    986997    EXAMPLES::
    987    
     998
    988999        sage: from sage.symbolic.relation import solve_ineq_univar
    9891000        sage: solve_ineq_univar(x-1/x>0)
    9901001        [[x > -1, x < 0], [x > 1]]
    991        
     1002
    9921003        sage: solve_ineq_univar(x^2-1/x>0)
    9931004        [[x < 0], [x > 1]]
    994        
     1005
    9951006        sage: solve_ineq_univar((x^3-1)*x<=0)
    9961007        [[x >= 0, x <= 1]]
    9971008
    9981009    ALGORITHM:
    999    
     1010
    10001011    Calls Maxima command solve_rat_ineq
    10011012
    10021013    AUTHORS:
    1003    
     1014
    10041015    - Robert Marik (01-2010)
    10051016    """
    10061017    ineqvar = ineq.variables()
     
    10091020    ineq0 = ineq._maxima_()
    10101021    ineq0.parent().eval("if solve_rat_ineq_loaded#true then (solve_rat_ineq_loaded:true,load(\"solve_rat_ineq.mac\")) ")
    10111022    sol = ineq0.solve_rat_ineq().sage()
    1012     if repr(sol)=="all": 
    1013         from sage.rings.infinity import Infinity 
    1014         sol = [ineqvar[0]<Infinity]       
     1023    if repr(sol)=="all":
     1024        from sage.rings.infinity import Infinity
     1025        sol = [ineqvar[0]<Infinity]
    10151026    return sol
    10161027
    10171028def solve_ineq_fourier(ineq,vars=None):
    10181029    """
    1019     Solves sytem of inequalities using Maxima and fourier elimination 
     1030    Solves sytem of inequalities using Maxima and fourier elimination
    10201031
    10211032    Can be used for system of linear inequalities and for some types
    10221033    of nonlinear inequalities. For examples see the section EXAMPLES
    10231034    below and http://maxima.cvs.sourceforge.net/viewvc/maxima/maxima/share/contrib/fourier_elim/rtest_fourier_elim.mac
    1024    
    1025    
     1035
     1036
    10261037    INPUT:
    1027    
    1028     - ``ineq`` - list with system of inequalities
    10291038
    1030     - ``vars`` - optionally list with variables for fourier elimination.
     1039    - ``ineq`` - list with system of inequalities
     1040
     1041    - ``vars`` - optionally list with variables for fourier elimination.
    10311042
    10321043    OUTPUT:
    1033    
    1034     - ``list`` - output is list of solutions as a list of simple inequalities 
    1035       output [A,B,C] means (A or B or C) each A, B, C is again a list and 
    1036       if A=[a,b], then A means (a and b). The list is empty if there is no 
     1044
     1045    - ``list`` - output is list of solutions as a list of simple inequalities
     1046      output [A,B,C] means (A or B or C) each A, B, C is again a list and
     1047      if A=[a,b], then A means (a and b). The list is empty if there is no
    10371048      solution.
    1038        
     1049
    10391050    EXAMPLES::
    1040    
     1051
    10411052        sage: from sage.symbolic.relation import solve_ineq_fourier
    1042         sage: y=var('y') 
     1053        sage: y=var('y')
    10431054        sage: solve_ineq_fourier([x+y<9,x-y>4],[x,y])
    1044         [[y + 4 < x, x < -y + 9, y < (5/2)]]       
     1055        [[y + 4 < x, x < -y + 9, y < (5/2)]]
    10451056        sage: solve_ineq_fourier([x+y<9,x-y>4],[y,x])
    10461057        [[y < min(x - 4, -x + 9)]]
    1047        
     1058
    10481059        sage: solve_ineq_fourier([x^2>=0])
    10491060        [[x < +Infinity]]
    10501061
     
    10611072        [[y < x, 0 < y]]
    10621073
    10631074    ALGORITHM:
    1064    
     1075
    10651076    Calls Maxima command fourier_elim
    10661077
    10671078    AUTHORS:
    1068    
     1079
    10691080    - Robert Marik (01-2010)
    10701081    """
    10711082    if vars is None:
     
    10801091        if not atom(x) and op(x)=\"or\" then args(x) \
    10811092        else [x]")
    10821093    sol = sol.or_to_list().sage()
    1083     if repr(sol) == "[emptyset]": 
     1094    if repr(sol) == "[emptyset]":
    10841095        sol = []
    10851096    if repr(sol) == "[universalset]":
    1086         from sage.rings.infinity import Infinity 
     1097        from sage.rings.infinity import Infinity
    10871098        sol = [[i<Infinity for i in vars]]
    10881099    return sol
    10891100
    10901101def solve_ineq(ineq, vars=None):
    10911102    """
    1092     Solves inequalities and systems of inequalities using Maxima. 
    1093     Switches between rational inequalities 
    1094     (sage.symbolic.relation.solve_ineq_rational) 
    1095     and fourier elimination (sage.symbolic.relation.solve_ineq_fouried). 
    1096     See the documentation of these functions for more details. 
    1097    
     1103    Solves inequalities and systems of inequalities using Maxima.
     1104    Switches between rational inequalities
     1105    (sage.symbolic.relation.solve_ineq_rational)
     1106    and fourier elimination (sage.symbolic.relation.solve_ineq_fouried).
     1107    See the documentation of these functions for more details.
     1108
    10981109    INPUT:
    10991110
    11001111    - ``ineq`` - one inequality or a list of inequalities
     
    11101121      for some types of nonlinear inequalities. See
    11111122      http://maxima.cvs.sourceforge.net/viewvc/maxima/maxima/share/contrib/fourier_elim/rtest_fourier_elim.mac
    11121123      for a big gallery of problems covered by this algorithm.
    1113        
     1124
    11141125    - ``vars`` - optional parameter with list of variables. This list
    11151126      is used only if fourier elimination is used. If omitted or if
    11161127      rational inequality is solved, then variables are determined
    11171128      automatically.
    11181129
    1119     OUTPUT: 
     1130    OUTPUT:
    11201131
    11211132    - ``list`` -- output is list of solutions as a list of simple inequalities
    1122       output [A,B,C] means (A or B or C) each A, B, C is again a list and 
     1133      output [A,B,C] means (A or B or C) each A, B, C is again a list and
    11231134      if A=[a,b], then A means (a and b).
    11241135
    11251136    EXAMPLES::
    1126    
     1137
    11271138        sage: from sage.symbolic.relation import solve_ineq
    11281139
    11291140    Inequalities in one variable. The variable is detected automatically::
    11301141
    1131         sage: solve_ineq(x^2-1>3)   
     1142        sage: solve_ineq(x^2-1>3)
    11321143        [[x < -2], [x > 2]]
    1133        
     1144
    11341145        sage: solve_ineq(1/(x-1)<=8)
    11351146        [[x < 1], [x >= (9/8)]]
    1136        
     1147
    11371148    System of inequalities with automatically detected inequalities::
    11381149
    11391150        sage: y=var('y')
     
    11501161        [[x < y, y < -x + 3, x < (3/2)]]
    11511162
    11521163    ALGORITHM:
    1153      
     1164
    11541165    Calls solve_ineq_fourier if inequalities are list and
    11551166    solve_ineq_univar of the inequality is symbolic expression. See
    11561167    the description of these commands for more details related to the
     
    11581169    there is no solution.
    11591170
    11601171    AUTHORS:
    1161    
     1172
    11621173    - Robert Marik (01-2010)
    11631174    """
    11641175    if isinstance(ineq,list):