Ticket #8616: trac_8616_symbolic_sage_correct.patch

File trac_8616_symbolic_sage_correct.patch, 52.2 KB (added by yuri.k, 11 years ago)

  • sage/calculus/all.py 

    # HG changeset patch
# User Yuri Karadzhov <yuri.karadzhov@gmail.com>
# Date 1269796851 0
# Node ID 8f637725000793f5850716f65617f0394f20f0f6
# Parent  46b1ee5997ac7e7faafdd86751fe7d121b6c1385
trac 8616: Symbolic type checking and expression parcing module

Provides unified interface for standard python types and sage specific types. Treats everything as symbolic expression which allows to check its type, take operator and operands and extract subexpressions by given types.

Provides autodetect dependent and independent variables in desolve()

diff -r 46b1ee5997ac -r 8f6377250007 sage/calculus/all.py
    a b  
    88
    99from functions import (wronskian,jacobian)
    1010
    11 from desolvers import (desolve, desolve_laplace, desolve_system,
    12                        eulers_method, eulers_method_2x2, 
     11from desolvers import (desolve, dsolve, desolve_laplace, dsolve_laplace, desolve_system,
     12                       eulers_method, eulers_method_2x2,
    1313                       eulers_method_2x2_plot, desolve_rk4, desolve_system_rk4)
    1414
    1515from var import (var, function, clear_vars)
    … …  
    2727    - a symbolic expression.
    2828
    2929    EXAMPLES::
    30    
     30
    3131        sage: a = symbolic_expression(3/2); a
    3232        3/2
    3333        sage: type(a)
    … …  
    5353        sage: symbolic_expression(E)
    5454        x*y + y^2 + y == x^3 + x^2 - 10*x - 10
    5555        sage: symbolic_expression(E) in SR
    56         True 
     56        True
    5757    """
    5858    from sage.symbolic.expression import Expression
    5959    from sage.symbolic.ring import SR

  • sage/calculus/desolvers.py 

    diff -r 46b1ee5997ac -r 8f6377250007 sage/calculus/desolvers.py
    a b  
    3232- ``eulers_method_2x2_plot`` - Plots the sequence of points obtained
    3333  from Euler's method.
    3434
    35 AUTHORS: 
     35AUTHORS:
    3636
    3737- David Joyner (3-2006) - Initial version of functions
    3838
    … …  
    4242
    4343- Robert Marik (10-2009) - Some bugfixes and enhancements
    4444
     45- Yuri Karadzhov (2010-03-27) - Autodetection of dependent and independent vars
     46
    4547"""
    4648
    4749##########################################################################
    … …  
    5860from sage.symbolic.expression import is_SymbolicEquation
    5961from sage.symbolic.ring import is_SymbolicVariable
    6062from sage.calculus.functional import diff
     63from sage.symbolic.mtype import *
    6164
    6265maxima = Maxima()
    6366
    64 def desolve(de, dvar, ics=None, ivar=None, show_method=False, contrib_ode=False):
     67def desolve(de, dvar=None, ics=None, ivar=None, show_method=False, contrib_ode=False):
    6568    """
    6669    Solves a 1st or 2nd order linear ODE via maxima. Including IVP and BVP.
    67    
     70
    6871    *Use* ``desolve? <tab>`` *if the output in truncated in notebook.*
    6972
    7073    INPUT:
    71    
     74
    7275    - ``de`` - an expression or equation representing the ODE
    73    
    74     - ``dvar`` - the dependent variable (hereafter called ``y``)
     76
     77    - ``dvar`` - (optional) the dependent variable (hereafter called ``y``)
    7578
    7679    - ``ics`` - (optional) the initial or boundary conditions
    7780
    … …  
    7982
    8083      - for a second-order equation, specify the initial ``x``, ``y``,
    8184        and ``dy/dx``
    82      
     85
    8386      - for a second-order boundary solution, specify initial and
    8487        final ``x`` and ``y`` initial conditions gives error if the
    8588        solution is not SymbolicEquation (as happens for example for
    8689        clairot equation)
    87        
     90
    8891    - ``ivar`` - (optional) the independent variable (hereafter called
    8992      x), which must be specified if there is more than one
    9093      independent variable in the equation.
    91                    
     94
    9295    - ``show_method`` - (optional) if true, then Sage returns pair
    9396      ``[solution, method]``, where method is the string describing
    9497      method which has been used to get solution (Maxima uses the
    … …  
    107110
    108111
    109112    OUTPUT:
    110    
     113
    111114    In most cases returns SymbolicEquation which defines the solution
    112115    implicitly.  If the result is in the form y(x)=... (happens for
    113116    linear eqs.), returns the right-hand side only.  The possible
    … …  
    118121
    119122        sage: x = var('x')
    120123        sage: y = function('y', x)
    121         sage: desolve(diff(y,x) + y - 1, y)
     124        sage: desolve(diff(y,x) + y - 1)
    122125        (c + e^x)*e^(-x)
    123    
     126
    124127    ::
    125128
    126         sage: f = desolve(diff(y,x) + y - 1, y, ics=[10,2]); f
     129        sage: f = desolve(diff(y,x) + y - 1, ics=[10,2]); f
    127130        (e^10 + e^x)*e^(-x)
    128131
    129132    ::
    130133
    131134        sage: plot(f)
    132        
     135
    133136    We can also solve second-order differential equations.::
    134    
     137
    135138        sage: x = var('x')
    136139        sage: y = function('y', x)
    137140        sage: de = diff(y,x,2) - y == x
    138         sage: desolve(de, y)
     141        sage: desolve(de)
    139142        k1*e^x + k2*e^(-x) - x
    140        
     143
    141144
    142145    ::
    143146
    … …  
    150153
    151154    ::
    152155
     156        sage: f = desolve(de, ics=[10,2,1]); f
     157        -x + 5*e^(-x + 10) + 7*e^(x - 10)
     158        sage: f(x=10)
     159        2
     160        sage: diff(f,x)(x=10)
     161        1
     162
     163    ::
     164
    153165        sage: de = diff(y,x,2) + y == 0
    154         sage: desolve(de, y)
     166        sage: desolve(de)
    155167        k1*sin(x) + k2*cos(x)
    156         sage: desolve(de, y, [0,1,pi/2,4])
     168        sage: desolve(de, ics=[0,1,pi/2,4])
    157169        4*sin(x) + cos(x)
    158        
    159         sage: desolve(y*diff(y,x)+sin(x)==0,y)
     170
     171        sage: desolve(y*diff(y,x)+sin(x)==0)
    160172        -1/2*y(x)^2 == c - cos(x)
    161173
    162174    Clairot equation: general and singular solutions::
    163175
    164         sage: desolve(diff(y,x)^2+x*diff(y,x)-y==0,y,contrib_ode=True,show_method=True)
     176        sage: desolve(diff(y,x)^2+x*diff(y,x)-y==0,contrib_ode=True,show_method=True)
    165177        [[y(x) == c^2 + c*x, y(x) == -1/4*x^2], 'clairault']
    166    
    167     For equations involving more variables we specify independent variable::
     178
     179    For equations involving more variables we shouldn't specify independent variable if they are obvious::
    168180
    169181        sage: a,b,c,n=var('a b c n')
    170         sage: desolve(x^2*diff(y,x)==a+b*x^n+c*x^2*y^2,y,ivar=x,contrib_ode=True)
     182        sage: desolve(x^2*diff(y,x)==a+b*x^n+c*x^2*y^2,contrib_ode=True)
    171183        [[y(x) == 0, (b*x^(n - 2) + a/x^2)*c^2*u == 0]]
    172         sage: desolve(x^2*diff(y,x)==a+b*x^n+c*x^2*y^2,y,ivar=x,contrib_ode=True,show_method=True)
     184        sage: desolve(x^2*diff(y,x)==a+b*x^n+c*x^2*y^2,contrib_ode=True,show_method=True)
    173185        [[[y(x) == 0, (b*x^(n - 2) + a/x^2)*c^2*u == 0]], 'riccati']
    174186
    175187
    176188    Higher orded, not involving independent variable::
    177189
    178         sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y).expand()
     190        sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0).expand()
    179191        1/6*y(x)^3 + k1*y(x) == k2 + x
    180         sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,1,3]).expand()
     192        sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,ics=[0,1,1,3]).expand()
    181193        1/6*y(x)^3 - 5/3*y(x) == x - 3/2
    182         sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,1,3],show_method=True)
     194        sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,ics=[0,1,1,3],show_method=True)
    183195        [1/6*y(x)^3 - 5/3*y(x) == x - 3/2, 'freeofx']
    184196
    185197    Separable equations - Sage returns solution in implicit form::
    186198
    187         sage: desolve(diff(y,x)*sin(y) == cos(x),y)
     199        sage: desolve(diff(y,x)*sin(y) == cos(x))
    188200        -cos(y(x)) == c + sin(x)
    189         sage: desolve(diff(y,x)*sin(y) == cos(x),y,show_method=True)
     201        sage: desolve(diff(y,x)*sin(y) == cos(x),show_method=True)
    190202        [-cos(y(x)) == c + sin(x), 'separable']
    191         sage: desolve(diff(y,x)*sin(y) == cos(x),y,[pi/2,1])
     203        sage: desolve(diff(y,x)*sin(y) == cos(x),ics=[pi/2,1])
    192204        -cos(y(x)) == sin(x) - cos(1) - 1
    193205
    194206    Linear equation - Sage returns the expression on the right hand side only::
    195207
    196         sage: desolve(diff(y,x)+(y) == cos(x),y)
     208        sage: desolve(diff(y,x)+(y) == cos(x))
    197209        1/2*((sin(x) + cos(x))*e^x + 2*c)*e^(-x)
    198         sage: desolve(diff(y,x)+(y) == cos(x),y,show_method=True)
     210        sage: desolve(diff(y,x)+(y) == cos(x),show_method=True)
    199211        [1/2*((sin(x) + cos(x))*e^x + 2*c)*e^(-x), 'linear']
    200         sage: desolve(diff(y,x)+(y) == cos(x),y,[0,1])
     212        sage: desolve(diff(y,x)+(y) == cos(x),ics=[0,1])
    201213        1/2*(e^x*sin(x) + e^x*cos(x) + 1)*e^(-x)
    202214
    203215    This ODE with separated variables is solved as
    204216    exact. Explanation - factor does not split `e^{x-y}` in Maxima
    205217    into `e^{x}e^{y}`::
    206218
    207         sage: desolve(diff(y,x)==exp(x-y),y,show_method=True)
     219        sage: desolve(diff(y,x)==exp(x-y),show_method=True)
    208220        [-e^x + e^y(x) == c, 'exact']
    209221
    210222    You can solve Bessel equations. You can also use initial
    … …  
    212224    condition at x=0, since this point is singlar point of the
    213225    equation. Anyway, if the solution should be bounded at x=0, then
    214226    k2=0.::
    215            
    216         sage: desolve(x^2*diff(y,x,x)+x*diff(y,x)+(x^2-4)*y==0,y)
     227
     228        sage: desolve(x^2*diff(y,x,x)+x*diff(y,x)+(x^2-4)*y==0)
    217229        k1*bessel_j(2, x) + k2*bessel_y(2, x)
    218    
     230
    219231    Difficult ODE produces error::
    220232
    221233        sage: desolve(sqrt(y)*diff(y,x)+e^(y)+cos(x)-sin(x+y)==0,y) # not tested
    222234        Traceback (click to the left for traceback)
    223235        ...
    224236        NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True."
    225        
     237
    226238    Difficult ODE produces error - moreover, takes a long time ::
    227239
    228240        sage: desolve(sqrt(y)*diff(y,x)+e^(y)+cos(x)-sin(x+y)==0,y,contrib_ode=True) # not tested
    … …  
    233245        [[y(x) == c + log(x), y(x) == c*e^x], 'factor']
    234246
    235247    ::
    236    
     248
    237249        sage: desolve(diff(y,x)==(x+y)^2,y,contrib_ode=True,show_method=True)
    238250        [[[x == c - arctan(sqrt(t)), y(x) == -x - sqrt(t)], [x == c + arctan(sqrt(t)), y(x) == -x + sqrt(t)]], 'lagrange']
    239251
    … …  
    245257        Traceback (click to the left for traceback)
    246258        ...
    247259        NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True."
    248        
     260
    249261    ::
    250262
    251263        sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,2],show_method=True) # not tested
    … …  
    267279        3*((e^(1/2*pi) - 2)*x/pi + 1)*e^(-x) + 1/2*sin(x)
    268280        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,[0,3,pi/2,2],show_method=True)
    269281        [3*((e^(1/2*pi) - 2)*x/pi + 1)*e^(-x) + 1/2*sin(x), 'variationofparameters']
    270        
     282
    271283        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y)
    272284        (k2*x + k1)*e^(-x)
    273285        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,show_method=True)
    … …  
    281293        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,[0,3,pi/2,2],show_method=True)
    282294        [(2*(2*e^(1/2*pi) - 3)*x/pi + 3)*e^(-x), 'constcoeff']
    283295
    284        
     296
    285297    This equation can be solved within Maxima but not within Sage. It
    286298    needs assumptions assume(x>0,y>0) and works in Maxima, but not in Sage::
    287299
    … …  
    289301        sage: assume(y>0) # not tested
    290302        sage: desolve(x*diff(y,x)-x*sqrt(y^2+x^2)-y,y,show_method=True) # not tested
    291303
    292        
     304
    293305    TESTS:
    294    
     306
    295307    Trac #6479 fixed::
    296308
    297309        sage: x = var('x')
    … …  
    300312        x
    301313
    302314    ::
    303    
     315
    304316        sage: desolve( diff(y,x,x) == 0, y, [0,1,1])
    305317        x + 1
    306318
    307319
    308320    AUTHORS:
    309  
     321
    310322    - David Joyner (1-2006)
    311323
    312324    - Robert Bradshaw (10-2008)
    313325
    314326    - Robert Marik (10-2009)
    315327
     328    - Yuri Karadzhov (2010-03-27)
     329
    316330    """
    317331    if is_SymbolicEquation(de):
    318332        de = de.lhs() - de.rhs()
    319     if is_SymbolicVariable(dvar):
    320         raise ValueError, "You have to declare dependent variable as a function, eg. y=function('y',x)"
    321     # for backwards compatibility
    322     if isinstance(dvar, list):
    323         dvar, ivar = dvar
    324     elif ivar is None:
    325         ivars = de.variables()
    326         ivars = [t for t in ivars if t is not dvar]
     333    diff_list = map(ops,indets(de,'_diff'))
     334    if len(diff_list) == 0:
     335        raise ValueError, "First argument is not differential equation."
     336    if dvar is None:
     337        dvars = set([x[0] for x in diff_list])
     338        if len(dvars) != 1:
     339            raise ValueError, "Unable to determine dependent variable, please specify."
     340        dvar=dvars.pop()
     341    else:
     342        if is_SymbolicVariable(dvar):
     343            raise ValueError, "You have to declare dependent variable as a function, eg. y=function('y',x)"
     344        # for backwards compatibility
     345        if isinstance(dvar, list):
     346            dvar, ivar = dvar
     347    if ivar is None:
     348        ivars = set()
     349        for el in diff_list:
     350            ivars = ivars.union(el[1:])
    327351        if len(ivars) != 1:
    328352            raise ValueError, "Unable to determine independent variable, please specify."
    329         ivar = ivars[0]
     353        ivar = ivars.pop()
    330354    def sanitize_var(exprs):
    331         return exprs.replace("'"+dvar_str+"("+ivar_str+")",dvar_str)   
     355        return exprs.replace("'"+dvar_str+"("+ivar_str+")",dvar_str)
    332356    dvar_str=dvar.operator()._maxima_().str()
    333357    ivar_str=ivar._maxima_().str()
    334358    de00 = de._maxima_().str()
    335359    de0 = sanitize_var(de00)
    336360    ode_solver="ode2"
    337     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)   
     361    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)
    338362    # we produce string like this
    339363    # ode2('diff(y,x,2)+2*'diff(y,x,1)+y-cos(x),y(x),x)
    340364    soln = maxima(cmd)
    … …  
    343367        if contrib_ode:
    344368            ode_solver="contrib_ode"
    345369            maxima("load('contrib_ode)")
    346             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)   
     370            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)
    347371            # we produce string like this
    348372            # (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))
    349373            soln = maxima(cmd)
    350374            if str(soln).strip() == 'false':
    351                 raise NotImplementedError, "Maxima was unable to solve this ODE."   
     375                raise NotImplementedError, "Maxima was unable to solve this ODE."
    352376        else:
    353377            raise NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True."
    354378
    … …  
    357381
    358382    if (ics is not None):
    359383        if not is_SymbolicEquation(soln.sage()):
    360              raise NotImplementedError, "Maxima was unable to use initial condition for this equation (%s)"%(maxima_method.str())   
     384             raise NotImplementedError, "Maxima was unable to use initial condition for this equation (%s)"%(maxima_method.str())
    361385        if len(ics) == 2:
    362386            tempic=(ivar==ics[0])._maxima_().str()
    363387            tempic=tempic+","+(dvar==ics[1])._maxima_().str()
    … …  
    367391            # (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))
    368392            soln=maxima(cmd)
    369393        if len(ics) == 3:
    370             #fixed ic2 command from Maxima - we have to ensure that %k1, %k2 do not depend on variables, should be removed when fixed in Maxima 
     394            #fixed ic2 command from Maxima - we have to ensure that %k1, %k2 do not depend on variables, should be removed when fixed in Maxima
    371395            maxima("ic2_sage(soln,xa,ya,dya):=block([programmode:true,backsubst:true,singsolve:true,temp,%k2,%k1,TEMP_k], \
    372396                noteqn(xa), noteqn(ya), noteqn(dya), boundtest('%k1,%k1), boundtest('%k2,%k2), \
    373397                temp: lhs(soln) - rhs(soln), \
    … …  
    384408            # (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))
    385409            soln=maxima(cmd)
    386410            if str(soln).strip() == 'false':
    387                 raise NotImplementedError, "Maxima was unable to solve this IVP. Remove the initial contition to get the general solution."   
     411                raise NotImplementedError, "Maxima was unable to solve this IVP. Remove the initial contition to get the general solution."
    388412        if len(ics) == 4:
    389             #fixed bc2 command from Maxima - we have to ensure that %k1, %k2 do not depend on variables, should be removed when fixed in Maxima 
     413            #fixed bc2 command from Maxima - we have to ensure that %k1, %k2 do not depend on variables, should be removed when fixed in Maxima
    390414            maxima("bc2_sage(soln,xa,ya,xb,yb):=block([programmode:true,backsubst:true,singsolve:true,temp,%k1,%k2,TEMP_k], \
    391415                noteqn(xa), noteqn(ya), noteqn(xb), noteqn(yb), boundtest('%k1,%k1), boundtest('%k2,%k2), \
    392416                TEMP_k:solve([subst([xa,ya],soln), subst([xb,yb],soln)], [%k1,%k2]), \
    … …  
    394418                temp: maplist(lambda([zz], subst(zz,soln)),TEMP_k), \
    395419                if length(temp)=1 then return(first(temp)) else return(temp))")
    396420            cmd="bc2_sage(%s(%s,%s,%s),%s,%s=%s,%s,%s=%s)"%(ode_solver,de00,dvar_str,ivar_str,(ivar==ics[0])._maxima_().str(),dvar_str,(ics[1])._maxima_().str(),(ivar==ics[2])._maxima_().str(),dvar_str,(ics[3])._maxima_().str())
    397             cmd="(TEMP:%s,substitute(%s=%s(%s),TEMP))"%(cmd,dvar_str,dvar_str,ivar_str)           
     421            cmd="(TEMP:%s,substitute(%s=%s(%s),TEMP))"%(cmd,dvar_str,dvar_str,ivar_str)
    398422            cmd=sanitize_var(cmd)
    399423            # we produce string like this
    400424            # (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))
    401425            soln=maxima(cmd)
    402426            if str(soln).strip() == 'false':
    403                 raise NotImplementedError, "Maxima was unable to solve this BVP. Remove the initial contition to get the general solution."   
     427                raise NotImplementedError, "Maxima was unable to solve this BVP. Remove the initial contition to get the general solution."
    404428
    405429    soln=soln.sage()
    406430    if is_SymbolicEquation(soln) and soln.lhs() == dvar:
    … …  
    409433        soln = soln.rhs()
    410434    if show_method:
    411435        return [soln,maxima_method.str()]
    412     else:   
     436    else:
    413437        return soln
    414438
     439def dsolve(de, ics=None, dvar=None, ivar=None, show_method=False, contrib_ode=False):
     440    """
     441    ALIAS for desolve
     442    sage: y=function('y',x)
     443    sage: desolve(de,[0,1,pi/2,4])
     444    4*sin(x) + cos(x)
     445    """
     446    return desolve(de, dvar, ics, ivar, show_method, contrib_ode)
    415447
    416448#def desolve_laplace2(de,vars,ics=None):
    417449##     """
    … …  
    454486#    #return maxima.eval(cmd)
    455487#    return de0.desolve(vars[0]).rhs()
    456488
    457    
    458 def desolve_laplace(de, dvar, ics=None, ivar=None):
     489
     490def desolve_laplace(de, dvar=None, ics=None, ivar=None):
    459491    """
    460492    Solves an ODE using laplace transforms. Initials conditions are optional.
    461493
    462494    INPUT:
    463    
     495
    464496    - ``de`` - a lambda expression representing the ODE (eg, de =
    465497      diff(y,x,2) == diff(y,x)+sin(x))
    466498
    467     - ``dvar`` - the dependent variable (eg y)
     499    - ``dvar`` - (optional) the dependent variable (eg y)
    468500
    469501    - ``ivar`` - (optional) the independent variable (hereafter called
    470502      x), which must be specified if there is more than one
    … …  
    478510    Solution of the ODE as symbolic expression
    479511
    480512    EXAMPLES::
    481    
     513
    482514        sage: u=function('u',x)
    483515        sage: eq = diff(u,x) - exp(-x) - u == 0
    484516        sage: desolve_laplace(eq,u)
    485517        1/2*(2*u(0) + 1)*e^x - 1/2*e^(-x)
    486    
     518
    487519    We can use initial conditions::
    488520
    489521        sage: desolve_laplace(eq,u,ics=[0,3])
    490522        -1/2*e^(-x) + 7/2*e^x
    491523
    492     The initial conditions do not persist in the system (as they persisted 
     524    The initial conditions do not persist in the system (as they persisted
    493525    in previous versions)::
    494526
    495527        sage: desolve_laplace(eq,u)
    … …  
    501533        sage: eq = diff(f,x) + f == 0
    502534        sage: desolve_laplace(eq,f,[0,1])
    503535        e^(-x)
    504    
     536
    505537    ::
    506    
     538
    507539        sage: x = var('x')
    508540        sage: f = function('f', x)
    509541        sage: de = diff(f,x,x) - 2*diff(f,x) + f
    … …  
    523555        e^(-t) + 1
    524556        sage: soln(t=3)
    525557        e^(-3) + 1
    526    
    527     AUTHORS: 
     558
     559    AUTHORS:
    528560
    529561    - David Joyner (1-2006,8-2007)
    530    
     562
    531563    - Robert Marik (10-2009)
     564
     565    - Yuri Karadzhov (2010-03-27)
    532566    """
    533567    #This is the original code from David Joyner (inputs and outputs strings)
    534568    #maxima("de:"+de._repr_()+"=0;")
    … …  
    544578    ## verbatim copy from desolve - begin
    545579    if is_SymbolicEquation(de):
    546580        de = de.lhs() - de.rhs()
    547     if is_SymbolicVariable(dvar):
    548         raise ValueError, "You have to declare dependent variable as a function, eg. y=function('y',x)"
    549     # for backwards compatibility
    550     if isinstance(dvar, list):
    551         dvar, ivar = dvar
    552     elif ivar is None:
    553         ivars = de.variables()
    554         ivars = [t for t in ivars if t != dvar]
     581    diff_list = map(ops,indets(de,'_diff'))
     582    if len(diff_list) == 0:
     583        raise ValueError, "First argument is not differential equation."
     584    if dvar is None:
     585        dvars = set([x[0] for x in diff_list])
     586        if len(dvars) != 1:
     587            raise ValueError, "Unable to determine dependent variable, please specify."
     588        dvar=dvars.pop()
     589    else:
     590        if is_SymbolicVariable(dvar):
     591            raise ValueError, "You have to declare dependent variable as a function, eg. y=function('y',x)"
     592        # for backwards compatibility
     593        if isinstance(dvar, list):
     594            dvar, ivar = dvar
     595    if ivar is None:
     596        ivars = set()
     597        for el in diff_list:
     598            ivars = ivars.union(el[1:])
    555599        if len(ivars) != 1:
    556600            raise ValueError, "Unable to determine independent variable, please specify."
    557         ivar = ivars[0]
     601        ivar = ivars.pop()
    558602    ## verbatim copy from desolve - end
    559603
    560604    def sanitize_var(exprs):  # 'y(x) -> y(x)
    561         return exprs.replace("'"+str(dvar),str(dvar))   
     605        return exprs.replace("'"+str(dvar),str(dvar))
    562606    de0=de._maxima_()
    563607    cmd = sanitize_var("desolve("+de0.str()+","+str(dvar)+")")
    564608    soln=maxima(cmd).rhs()
    565609    if str(soln).strip() == 'false':
    566         raise NotImplementedError, "Maxima was unable to solve this ODE."   
     610        raise NotImplementedError, "Maxima was unable to solve this ODE."
    567611    soln=soln.sage()
    568612    if ics!=None:
    569613        d = len(ics)
    570614        for i in range(0,d-1):
    571             soln=eval('soln.substitute(diff(dvar,ivar,i)('+str(ivar)+'=ics[0])==ics[i+1])') 
     615            soln=eval('soln.substitute(diff(dvar,ivar,i)('+str(ivar)+'=ics[0])==ics[i+1])')
    572616    return soln
    573    
    574    
     617
     618def dsolve_laplace(de, ics=None, dvar=None, ivar=None):
     619    """
     620    ALIAS for desolve_laplace
     621    sage: x = var('x')
     622        sage: f = function('f', x)
     623        sage: de = diff(f,x,x) - 2*diff(f,x) + f
     624        sage: desolve_laplace(de)
     625        -x*e^x*f(0) + x*e^x*D[0](f)(0) + e^x*f(0)
     626        sage: desolve_laplace(de,[0,1,2])
     627        x*e^x + e^x
     628    """
     629    return desolve_laplace(de, dvar, ics, ivar)
     630
    575631def desolve_system(des, vars, ics=None, ivar=None):
    576632    """
    577633    Solves any size system of 1st order ODE's. Initials conditions are optional.
    578634
    579635    INPUT:
    580    
     636
    581637    - ``des`` - list of ODEs
    582    
     638
    583639    - ``vars`` - list of dependent variables
    584    
     640
    585641    - ``ics`` - (optional) list of initial values for ivar and vars
    586    
     642
    587643    - ``ivar`` - (optional) the independent variable, which must be
    588644      specified if there is more than one independent variable in the
    589645      equation.
    … …  
    598654        sage: desolve_system([de1, de2], [x,y])
    599655        [x(t) == (x(0) - 1)*cos(t) - (y(0) - 1)*sin(t) + 1,
    600656         y(t) == (x(0) - 1)*sin(t) + (y(0) - 1)*cos(t) + 1]
    601          
     657
    602658    Now we give some initial conditions::
    603    
     659
    604660        sage: sol = desolve_system([de1, de2], [x,y], ics=[0,1,2]); sol
    605661        [x(t) == -sin(t) + 1, y(t) == cos(t) + 1]
    606662        sage: solnx, solny = sol[0].rhs(), sol[1].rhs()
    607663        sage: plot([solnx,solny],(0,1))
    608664        sage: parametric_plot((solnx,solny),(0,1))
    609665
    610     AUTHORS: 
     666    AUTHORS:
    611667
    612668    - Robert Bradshaw (10-2008)
    613669    """
    … …  
    637693        for dvar, ic in zip(dvars, ics[:1]):
    638694            dvar.atvalue(ivar==ivar_ic, dvar)
    639695    return soln
    640        
     696
    641697
    642698def desolve_system_strings(des,vars,ics=None):
    643699    r"""
    … …  
    687743        sage: P2 = parametric_plot((solnx,solny),(0,1))
    688744
    689745    Now type show(P1), show(P2) to view these.
    690                      
    691746
    692     AUTHORS:
     747
     748    AUTHORS:
    693749
    694750    - David Joyner (3-2006, 8-2007)
    695751    """
    … …  
    721777    last column.
    722778
    723779    *For pedagogical purposes only.*
    724    
     780
    725781    EXAMPLES::
    726    
     782
    727783        sage: from sage.calculus.desolvers import eulers_method
    728784        sage: x,y = PolynomialRing(QQ,2,"xy").gens()
    729785        sage: eulers_method(5*x+y-5,0,1,1/2,1)
    … …  
    777833        sage: P2 = line(pts)
    778834        sage: (P1+P2).show()
    779835
    780     AUTHORS: 
     836    AUTHORS:
    781837
    782838    - David Joyner
    783839    """
    … …  
    812868    next (last) column.
    813869
    814870    *For pedagogical purposes only.*
    815    
     871
    816872    EXAMPLES::
    817873
    818874        sage: from sage.calculus.desolvers import eulers_method_2x2
    … …  
    858914           3/4                 0.63                   -0.0078               -0.031                 0.11
    859915             1                 0.63                     0.020                0.079                0.071
    860916
    861     To numerically approximate y(1), where `y''+ty'+y=0`, `y(0)=1`, `y'(0)=0`:: 
     917    To numerically approximate y(1), where `y''+ty'+y=0`, `y(0)=1`, `y'(0)=0`::
    862918
    863919        sage: t,x,y=PolynomialRing(RR,3,"txy").gens()
    864920        sage: f = y; g = -x-y*t
    … …  
    870926           3/4                 0.88                     -0.15                -0.62                -0.10
    871927             1                 0.75                     -0.17                -0.68               -0.015
    872928
    873     AUTHORS: 
     929    AUTHORS:
    874930
    875931    - David Joyner
    876932    """
    … …  
    900956    `x(a)=x_0`, `y' = g(t,x,y)`, `y(a) = y_0`.
    901957
    902958    *For pedagogical purposes only.*
    903    
     959
    904960    EXAMPLES::
    905961
    906962        sage: from sage.calculus.desolvers import eulers_method_2x2_plot
    … …  
    930986def desolve_rk4_determine_bounds(ics,end_points=None):
    931987    """
    932988    Used to determine bounds for numerical integration.
    933    
     989
    934990    - If end_points is None, the interval for integration is from ics[0]
    935991      to ics[0]+10
    936992
    … …  
    9451001        sage: from sage.calculus.desolvers import desolve_rk4_determine_bounds
    9461002        sage: desolve_rk4_determine_bounds([0,2],1)
    9471003        (0, 1)
    948    
     1004
    9491005    ::
    9501006
    951         sage: desolve_rk4_determine_bounds([0,2]) 
     1007        sage: desolve_rk4_determine_bounds([0,2])
    9521008        (0, 10)
    953    
     1009
    9541010    ::
    9551011
    9561012        sage: desolve_rk4_determine_bounds([0,2],[-2])
    9571013        (-2, 0)
    958    
     1014
    9591015    ::
    9601016
    9611017        sage: desolve_rk4_determine_bounds([0,2],[-2,4])
    9621018        (-2, 4)
    963    
     1019
    9641020    """
    965     if end_points is None: 
     1021    if end_points is None:
    9661022        return((ics[0],ics[0]+10))
    9671023    if not isinstance(end_points,list):
    9681024        end_points=[end_points]
    … …  
    9701026        return (min(ics[0],end_points[0]),max(ics[0],end_points[0]))
    9711027    else:
    9721028        return (min(ics[0],end_points[0]),max(ics[0],end_points[1]))
    973    
     1029
    9741030
    9751031def desolve_rk4(de, dvar, ics=None, ivar=None, end_points=None, step=0.1, output='list', **kwds):
    9761032    """
    … …  
    9781034    equation. See also ``ode_solver``.
    9791035
    9801036    INPUT:
    981    
    982     input is similar to ``desolve`` command. The differential equation can be 
     1037
     1038    input is similar to ``desolve`` command. The differential equation can be
    9831039    written in a form close to the plot_slope_field or desolve command
    9841040
    9851041    - Variant 1 (function in two variables)
    986      
     1042
    9871043      - ``de`` - right hand side, i.e. the function `f(x,y)` from ODE `y'=f(x,y)`
    988      
     1044
    9891045      - ``dvar`` - dependent variable (symbolic variable declared by var)
    9901046
    9911047    - Variant 2 (symbolic equation)
    … …  
    9941050
    9951051      - ``dvar``` - dependent variable (declared as funciton of independent variable)
    9961052
    997     - Other parameters 
     1053    - Other parameters
    9981054
    9991055      - ``ivar`` - should be specified, if there are more variables or if the equation is autonomous
    10001056
    10011057      - ``ics`` - initial conditions in the form [x0,y0]
    1002      
     1058
    10031059      - ``end_points`` - the end points of the interval
    10041060
    10051061        - if end_points is a or [a], we integrate on between min(ics[0],a) and max(ics[0],a)
    … …  
    10081064        - if end_points is [a,b] we integrate on between min(ics[0],a) and max(ics[0],b)
    10091065
    10101066      - ``step`` - (optional, default:0.1) the length of the step (positive number)
    1011  
     1067
    10121068      - ``output`` - (optional, default: 'list') one of 'list',
    10131069        'plot', 'slope_field' (graph of the solution with slope field)
    10141070
    1015     OUTPUT: 
     1071    OUTPUT:
    10161072
    10171073    Returns a list of points, or plot produced by list_plot,
    10181074    optionally with slope field.
    … …  
    10301086
    10311087    Variant 1 for input - we can pass ODE in the form used by
    10321088    desolve function In this example we integrate bakwards, since
    1033     ``end_points < ics[0]``:: 
     1089    ``end_points < ics[0]``::
    10341090
    1035         sage: y=function('y',x) 
    1036         sage: desolve_rk4(diff(y,x)+y*(y-1) == x-2,y,ics=[1,1],step=0.5, end_points=0) 
     1091        sage: y=function('y',x)
     1092        sage: desolve_rk4(diff(y,x)+y*(y-1) == x-2,y,ics=[1,1],step=0.5, end_points=0)
    10371093        [[0.0, 8.90425710896], [0.5, 1.90932794536], [1, 1]]
    10381094
    10391095    Here we show how to plot simple pictures. For more advanced
    … …  
    10421098
    10431099        sage: x,y=var('x y')
    10441100        sage: P=desolve_rk4(y*(2-y),y,ics=[0,.1],ivar=x,output='slope_field',end_points=[-4,6],thickness=3)
    1045                
     1101
    10461102    ALGORITHM:
    10471103
    10481104    4th order Runge-Kutta method. Wrapper for command ``rk`` in
    10491105    Maxima's dynamics package.  Perhaps could be faster by using
    10501106    fast_float instead.
    10511107
    1052     AUTHORS: 
     1108    AUTHORS:
    10531109
    10541110    - Robert Marik (10-2009)
    10551111    """
    1056     if ics is None: 
     1112    if ics is None:
    10571113        raise ValueError, "No initial conditions, specify with ics=[x0,y0]."
    10581114
    10591115    if ivar is None:
    … …  
    11101166        YMIN=sol[0][1]
    11111167        XMAX=XMIN
    11121168        YMAX=YMIN
    1113         for s,t in sol: 
     1169        for s,t in sol:
    11141170            if s>XMAX:XMAX=s
    11151171            if s<XMIN:XMIN=s
    11161172            if t>YMAX:YMAX=t
    … …  
    11241180    r"""
    11251181    Solves numerically system of first-order ordinary differential
    11261182    equations using the 4th order Runge-Kutta method. Wrapper for
    1127     Maxima command ``rk``. See also ``ode_solver``. 
     1183    Maxima command ``rk``. See also ``ode_solver``.
    11281184
    11291185    INPUT:
    1130    
     1186
    11311187    input is similar to desolve_system and desolve_rk4 commands
    1132    
     1188
    11331189    - ``des`` - right hand sides of the system
    11341190
    11351191    - ``vars`` - dependent variables
    … …  
    11391195      missing
    11401196
    11411197    - ``ics`` - initial conditions in the form [x0,y01,y02,y03,....]
    1142    
     1198
    11431199    - ``end_points`` - the end points of the interval
    1144    
     1200
    11451201      - if end_points is a or [a], we integrate on between min(ics[0],a) and max(ics[0],a)
    11461202      - if end_points is None, we use end_points=ics[0]+10
    11471203
    … …  
    11641220        sage: P=desolve_system_rk4([x*(1-y),-y*(1-x)],[x,y],ics=[0,0.5,2],ivar=t,end_points=20)
    11651221        sage: Q=[ [i,j] for i,j,k in P]
    11661222        sage: LP=list_plot(Q)
    1167        
     1223
    11681224        sage: Q=[ [j,k] for i,j,k in P]
    11691225        sage: LP=list_plot(Q)
    1170    
     1226
    11711227    ALGORITHM:
    11721228
    11731229    4th order Runge-Kutta method. Wrapper for command ``rk`` in Maxima's
    11741230    dynamics package.  Perhaps could be faster by using ``fast_float``
    11751231    instead.
    11761232
    1177     AUTHOR: 
     1233    AUTHOR:
    11781234
    11791235    - Robert Marik (10-2009)
    11801236    """
    11811237
    1182     if ics is None: 
     1238    if ics is None:
    11831239        raise ValueError, "No initial conditions, specify with ics=[x0,y01,y02,...]."
    11841240
    11851241    ivars = set([])
    … …  
    11991255    x0=ics[0]
    12001256    icss = [ics[i]._maxima_().str() for i in range(1,len(ics))]
    12011257    icstr = "[" + ",".join(icss) + "]"
    1202     step=abs(step) 
     1258    step=abs(step)
    12031259
    12041260    maxima("load('dynamics)")
    12051261    lower_bound,upper_bound=desolve_rk4_determine_bounds(ics,end_points)

  • sage/symbolic/all.py 

    diff -r 46b1ee5997ac -r 8f6377250007 sage/symbolic/all.py
    a b  
    88
    99from sage.symbolic.relation import solve, solve_mod, solve_ineq
    1010from sage.symbolic.assumptions import assume, forget, assumptions
     11
     12from sage.symbolic.mtype import *

  • new file sage/symbolic/mtype.py 

    diff -r 46b1ee5997ac -r 8f6377250007 sage/symbolic/mtype.py
    - +  
     1r"""
     2Symbolic type checking and expression parcing module
     3
     4Provides unified interface for standard python types and sage specific types.
     5Treats everything as symbolic expression which allows to check its type, take
     6operator and operands and extract subexpressions by given types.
     7
     8AUTHORS:
     9
     10- Yuri Karadzhov (2010-03-27) initial version
     11
     12EXAMPLES:
     13
     14sage: x,a,b,c,t=var('x,a,b,c,t')
     15sage: g=function('g',x)
     16sage: f=function('f',t,x)
     17sage: data=[x+sin(x+3)+tan(x+f/2)+ln(g+5)+exp(2*x), diff(f,x,t)+diff(g,x), 3,
     18vector((1,2,x)), vector((x, a, b, f)), 1]
     19sage:data
     20[x + e^(2*x) + log(g(x) + 5) + sin(x + 3) + tan(x + 1/2*f(t, x)), D[0](g)(x)
     21+ D[0, 1](f)(t, x), 3, (1, 2, x), (x, a, b, f(t, x)), 1]
     22sage: mtype(data)
     23'list.collection'
     24sage: map(mtype, data)
     25['+.arithmetic', '+.arithmetic', 'integer.number', 'symbolic.3.vector',
     26'symbolic.4.vector', 'integer.number']
     27sage: indets(data, '_symbolic')
     28Warning: vector (1, 2, x) was converted to tuple
     29Warning: vector (x, a, b, f(t, x)) was converted to tuple
     30set([x, g(x), f(t, x), t, (1, 2, x), (x, a, b, f(t, x)), a, b])
     31sage: indets(data,'_integer _symbolic')
     32Warning: vector (1, 2, x) was converted to tuple
     33Warning: vector (x, a, b, f(t, x)) was converted to tuple
     34set([1, 2, 3, 5, t, g(x), (1, 2, x), b, a, f(t, x), (x, a, b, f(t, x)), x])
     35sage: indets(data,'_.3.vector')
     36Warning: vector (1, 2, x) was converted to tuple
     37set([(1, 2, x)])
     38sage: indets(data,'_.4.vector')
     39Warning: vector (x, a, b, f(t, x)) was converted to tuple
     40set([(x, a, b, f(t, x))])
     41sage: vec=indets(data,'_.vector')
     42Warning: vector (1, 2, x) was converted to tuple
     43Warning: vector (x, a, b, f(t, x)) was converted to tuple
     44sage: vec
     45set([(x, a, b, f(t, x)), (1, 2, x)])
     46sage: ops(vec)
     47[(x, a, b, f(t, x)), (1, 2, x)]
     48sage: nops(vec)
     492
     50sage: v=ops(vec,0)
     51sage: v
     52(x, a, b, f(t, x))
     53sage: mtype(v)
     54'tuple.collection'
     55sage: mtype(v,'_.collection')
     56True
     57sage: mtype(v,'_list')
     58False
     59sage: mtype(v,'_tuple.')
     60True
     61sage: op(v)
     62<function tuple_operator at 0x3ab5050>
     63sage: ops(v)
     64[x, a, b, f(t, x)]
     65sage: op(v)(ops(v))
     66(x, a, b, f(t, x))
     67sage: mtype(op(v)(ops(v)),'_tuple.')
     68True
     69sage: ops(v,-1)
     70f(t, x)
     71sage: fu=ops(v,-1)
     72sage: fu
     73f(t, x)
     74sage: ops(fu)
     75[t, x]
     76sage: op(fu)(ops(fu))
     77f(t, x)
     78"""
     79#*****************************************************************************
     80#       Copyright (C) 2010 Yuri Karadzhov <yuri.karadzhov@gmail.com>
     81#
     82#  Distributed under the terms of the GNU General Public License (GPL)
     83#                  http://www.gnu.org/licenses/
     84#*****************************************************************************
     85
     86
     87import operator
     88import sage.symbolic.operators
     89import sage.matrix.matrix_symbolic_dense
     90from sage.symbolic.ring import is_SymbolicVariable
     91
     92def mtype(expr, tp=None):
     93    """
     94    RETURN symbolic type of expression
     95
     96    INPUT:
     97
     98    - ``expr`` - an expression
     99    - ``tp`` - (optional) symbolic type or types
     100
     101    OUTPUT:
     102
     103    symbolic type of expression
     104
     105    EXAMPLES:
     106
     107    _ is wildcard for type
     108    _type means "anything"type"anything"
     109
     110    sage: x=var('x')
     111    sage: eq=x+1
     112    sage: mtype(eq)
     113    '+.relation'
     114
     115    sage: x=var('x')
     116    sage: eq=2*x
     117    sage: mtype(eq)
     118    '*.relation'
     119
     120    sage: x=var('x')
     121    sage: eq=x+1
     122    sage: mtype(eq,'_+')
     123    True
     124
     125    sage: x=var('x')
     126    sage: eq=x+1
     127    sage: mtype(eq,'+')
     128    False
     129
     130    sage: x=var('x')
     131    sage: eq=x+1
     132    sage: mtype(eq,'+.relation')
     133    True
     134
     135    sage: x=var('x')
     136    sage: eq=2*x
     137    sage: mtype(eq,'_+')
     138    False
     139
     140    sage: x=var('x')
     141    sage: eq=2*x
     142    sage: mtype(eq,'*')
     143    True
     144
     145    sage: a,b=var('a b')
     146    sage: mtype(a>=b,'>=.relation')
     147    True
     148    sage: mtype(a>=b,'>.relation')
     149    False
     150    sage: mtype(a>=b,'==.relation')
     151    False
     152    sage: mtype(a>=b,'_=')
     153    True
     154    sage: mtype(a>=b,'_>')
     155    True
     156    sage: mtype(a>=b,'_>.')
     157    False
     158    sage: mtype(a>=b,'_>=.')
     159    True
     160
     161    AUTHORS:
     162
     163    - Yuri Karadzhov (2010-03-27)
     164
     165    """
     166    if not tp is None:#improove wildcards
     167        tp=tp.split()
     168        ist=False
     169        for t in tp:
     170            if t.startswith('_'):
     171                ist=ist or (mtype(expr).find(t[1:])!=-1)
     172            else: ist=ist or t==mtype(expr)
     173        return ist
     174    if isinstance(expr,sage.symbolic.expression.Expression):#Improve type names
     175        if is_SymbolicVariable(expr):
     176            return 'symbolic.variable'
     177        t=expr.operator()
     178        if t is operator.eq:
     179            return '==.relation'
     180        if t is operator.ge:
     181            return '>=.relation'
     182        if t is operator.le:
     183            return '<=.relation'
     184        if t is operator.ne:
     185            return '!=.relation'
     186        if t is operator.gt:
     187            return '>.relation'
     188        if t is operator.lt:
     189            return '<.relation'
     190        if t is operator.add:
     191            return '+.arithmetic'
     192        if t is operator.mul:
     193            return '*.arithmetic'
     194        if t is operator.pow:
     195            return '^.arithmetic'
     196        if t is None:
     197            return mtype(expr.pyobject())
     198        if isinstance(t,sage.symbolic.operators.FDerivativeOperator):
     199            return 'diff'
     200        if isinstance(t,sage.symbolic.function_factory.SymbolicFunction):
     201            return 'symbolic.function'
     202        if isinstance(t,sage.functions.log.Function_log):
     203            return 'log.logarithmic.standard.function'
     204        if isinstance(t,sage.functions.log.Function_polylog):
     205            return 'polylog.logarithmic.standard.function'
     206        if isinstance(t,sage.functions.log.Function_dilog):
     207            return 'dilog.logarithmic.standard.function'
     208        if isinstance(t,sage.functions.log.Function_exp):
     209            return 'exp.logarithmic.standard.function'
     210        if isinstance(t,sage.functions.trig.Function_sin):
     211            return 'sin.trig.standard.function'
     212        if isinstance(t,sage.functions.trig.Function_cos):
     213            return 'cos.trig.standard.function'
     214        if isinstance(t,sage.functions.trig.Function_sec):
     215            return 'sec.trig.standard.function'
     216        if isinstance(t,sage.functions.trig.Function_csc):
     217            return 'csc.trig.standard.function'
     218        if isinstance(t,sage.functions.trig.Function_cot):
     219            return 'cot.trig.standard.function'
     220        if isinstance(t,sage.functions.trig.Function_tan):
     221            return 'tan.trig.standard.function'
     222        if isinstance(t,sage.functions.trig.Function_arcsin):
     223            return 'arcsin.trig.standard.function'
     224        if isinstance(t,sage.functions.trig.Function_arccos):
     225            return 'arccos.trig.standard.function'
     226        if isinstance(t,sage.functions.trig.Function_arcsec):
     227            return 'arcsec.trig.standard.function'
     228        if isinstance(t,sage.functions.trig.Function_arccsc):
     229            return 'arccsc.trig.standard.function'
     230        if isinstance(t,sage.functions.trig.Function_arccot):
     231            return 'arccot.trig.standard.function'
     232        if isinstance(t,sage.functions.trig.Function_arctan):
     233            return 'arctan.trig.standard.function'
     234        if isinstance(t,sage.functions.trig.Function_arctan2):
     235            return 'arctan2.trig.standard.function'
     236        if isinstance(t,sage.functions.hyperbolic.Function_sinh):
     237            return 'sinh.hyperbolic.standard.function'
     238        if isinstance(t,sage.functions.hyperbolic.Function_cosh):
     239            return 'cosh.hyperbolic.standard.function'
     240        if isinstance(t,sage.functions.hyperbolic.Function_sech):
     241            return 'sech.hyperbolic.standard.function'
     242        if isinstance(t,sage.functions.hyperbolic.Function_csch):
     243            return 'csch.hyperbolic.standard.function'
     244        if isinstance(t,sage.functions.hyperbolic.Function_coth):
     245            return 'coth.hyperbolic.standard.function'
     246        if isinstance(t,sage.functions.hyperbolic.Function_tanh):
     247            return 'tanh.hyperbolic.standard.function'
     248        if isinstance(t,sage.functions.hyperbolic.Function_arcsinh):
     249            return 'arcsinh.hyperbolic.standard.function'
     250        if isinstance(t,sage.functions.hyperbolic.Function_arccosh):
     251            return 'arccosh.hyperbolic.standard.function'
     252        if isinstance(t,sage.functions.hyperbolic.Function_arcsech):
     253            return 'arcsech.hyperbolic.standard.function'
     254        if isinstance(t,sage.functions.hyperbolic.Function_arccsch):
     255            return 'arccsch.hyperbolic.standard.function'
     256        if isinstance(t,sage.functions.hyperbolic.Function_arccoth):
     257            return 'arccoth.hyperbolic.standard.function'
     258        if isinstance(t,sage.functions.hyperbolic.Function_arctanh):
     259            return 'arctanh.hyperbolic.standard.function'
     260        #TODO complete list of types
     261        else: return 'unknown'
     262    else:
     263        if isinstance(expr,set):
     264            return 'set.collection'
     265        if isinstance(expr,frozenset):
     266            return 'frozenset.collection'
     267        if isinstance(expr,list):
     268            return 'list.collection'
     269        if isinstance(expr,tuple):
     270            return 'tuple.collection'
     271        if isinstance(expr,dict):
     272            return 'dict.collection'
     273        if isinstance(expr,str):
     274            return 'string'
     275        if isinstance(expr,sage.rings.integer.Integer):
     276            return 'integer.number'
     277        if isinstance(expr,sage.rings.rational.Rational):
     278            return 'rational.number'
     279        if isinstance(expr,sage.rings.real_mpfr.RealLiteral):
     280            return 'real.number'
     281        if isinstance(expr,sage.rings.real_mpfr.RealNumber):
     282            return 'real.number'
     283        if isinstance(expr,sage.rings.complex_number.ComplexNumber):
     284            return 'complex.number'
     285        if isinstance(expr,sage.matrix.matrix_integer_dense.Matrix_integer_dense):
     286            return 'integer.{0}.matrix'.format(expr.parent().dims()).replace(' ','')
     287        if isinstance(expr,sage.matrix.matrix_rational_dense.Matrix_rational_dense):
     288            return 'rational.{0}.matrix'.format(expr.parent().dims()).replace(' ','')
     289        if isinstance(expr,sage.matrix.matrix_symbolic_dense.Matrix_symbolic_dense):
     290            return 'symbolic.{0}.matrix'.format(expr.parent().dims()).replace(' ','')
     291        if isinstance(expr,sage.matrix.matrix_generic_dense.Matrix_generic_dense):
     292            mp=expr.parent()
     293            b=mp.base()
     294            d=mp.dims()
     295            if isinstance(b,sage.rings.real_mpfr.RealField_class):
     296                return 'real.{0}.matrix'.format(d).replace(' ','')
     297            if isinstance(b,sage.rings.complex_field.ComplexField_class):
     298                return 'real.{0}.matrix'.format(d).replace(' ','')
     299        if isinstance(expr,sage.modules.vector_integer_dense.Vector_integer_dense):
     300            return 'integer.{0}.vector'.format(expr.parent().dimension())
     301        if isinstance(expr,sage.modules.vector_rational_dense.Vector_rational_dense):
     302            return 'rational,{0}.vector'.format(expr.parent().dimension())
     303        if isinstance(expr,sage.modules.free_module_element.FreeModuleElement_generic_dense):
     304            vp=expr.parent()
     305            b=vp.base()
     306            d=vp.dimension()
     307            if isinstance(b,sage.rings.real_mpfr.RealField_class):
     308                return 'real.{0}.vector'.format(d)
     309            if isinstance(b,sage.rings.complex_field.ComplexField_class):
     310                return 'real.{0}.vector'.format(d)
     311            if isinstance(b,sage.symbolic.ring.SymbolicRing):
     312                return 'symbolic.{0}.vector'.format(d)
     313        #TODO complete list of types
     314        else: return 'unknown'
     315
     316
     317def stype(expr):
     318    """
     319    RETURN short type of expression
     320
     321    INPUT:
     322
     323    - ``expr`` - an expression
     324
     325    OUTPUT:
     326
     327    short symbolic type of expression
     328
     329    EXAMPLES:
     330
     331    sage: x=var('x')
     332    sage: eq=x+1
     333    sage: stype(eq)
     334    '+'
     335
     336    sage: x=var('x')
     337    sage: eq=2*x
     338    sage: stype(eq)
     339    '*'
     340
     341    sage: stype(1)
     342    'real'
     343
     344    AUTHORS:
     345
     346    - Yuri Karadzhov (2010-03-27)
     347
     348    """
     349    t = mtype(expr)
     350    pos = t.find('.')
     351    if pos!=-1: return t[0:pos]
     352    else: return t
     353
     354
     355def wtype(expr):
     356    """
     357    RETURN general type of expression
     358
     359    INPUT:
     360
     361    - ``expr`` - an expression
     362
     363    OUTPUT:
     364
     365    general type of expression
     366
     367    EXAMPLES:
     368
     369    sage: x=var('x')
     370    sage: eq=x+1
     371    sage: stype(eq)
     372    'relation'
     373
     374    sage: x=var('x')
     375    sage: eq=2*x
     376    sage: stype(eq)
     377    'relation'
     378
     379    sage: stype(1)
     380    'number'
     381    """
     382    t = mtype(expr)
     383    pos = t.rfind('.')
     384    if pos!=-1: return t[pos+1:]
     385    else: return t
     386
     387
     388def ops(expr, n=None):#make possible to call like this ops(expr,1:3):=ops(expr)[1:3]
     389    """
     390    RETURN expression operands
     391
     392    INPUT:
     393
     394    - ``expr`` - an expression
     395    - ``n`` - (optional) operand number
     396
     397    OUTPUT:
     398
     399    a list of the expression operands
     400
     401    EXAMPLES:
     402
     403    sage: a=var('a')
     404    sage: t=1,2,a
     405    sage: ops(t)
     406    [1, 2, a]
     407
     408    sage: x=var('x')
     409    sage: y=function('y',x)
     410    sage: eq=diff(y,x)+y+x+1
     411    sage: ops(eq)
     412    [D[0](y)(x), y(x), x, 1]
     413
     414    sage: x,t=var('x t')
     415    sage: y=function('y',x,t)
     416    sage: eq=diff(y,x,t,t)
     417    sage: ops(eq)
     418    [y(x,t),x,t,t]
     419
     420    sage: x=var('x')
     421    sage: ops(x)
     422    []
     423
     424    sage: ops(1)
     425    []
     426
     427    AUTHORS:
     428
     429    - Yuri Karadzhov (2010-03-27)
     430
     431    """
     432    if not n is None:
     433        return ops(expr)[n]
     434    if isinstance(expr,sage.symbolic.expression.Expression):
     435        t=expr.operator()
     436        if isinstance(t,sage.symbolic.operators.FDerivativeOperator):
     437            oprs=[apply(t.function(),expr.operands())]
     438            oprs.extend([expr.operands()[x] for x in t.parameter_set()])
     439            return oprs
     440        else: return expr.operands()
     441    else:
     442        if isinstance(expr,list):
     443            return expr
     444        if isinstance(expr,set):
     445            return list(expr)
     446        if isinstance(expr,tuple):
     447            return list(expr)
     448        if isinstance(expr,dict):
     449            return expr.items()
     450        if isinstance(expr,sage.matrix.matrix_integer_dense.Matrix_integer_dense):
     451            return list(expr)
     452        if isinstance(expr,sage.matrix.matrix_rational_dense.Matrix_rational_dense):
     453            return list(expr)
     454        if isinstance(expr,sage.matrix.matrix_generic_dense.Matrix_generic_dense):
     455            return list(expr)
     456        if isinstance(expr,sage.modules.vector_integer_dense.Vector_integer_dense):
     457            return list(expr)
     458        if isinstance(expr,sage.modules.vector_rational_dense.Vector_rational_dense):
     459            return list(expr)
     460        if isinstance(expr,sage.modules.free_module_element.FreeModuleElement_generic_dense):
     461            return list(expr)
     462        else: return []
     463
     464
     465#operands=ops
     466
     467
     468def nops(expr):
     469    """
     470    RETURN the number of operands
     471
     472    INPUT:
     473
     474    - ``expr`` - an expression
     475
     476    OUTPUT:
     477
     478    number of operands of the expression
     479
     480    EXAMPLES:
     481
     482    sage: a=var('a')
     483    sage: t=1,2,a
     484    sage: nops(t)
     485    3
     486
     487    sage: x=var('x')
     488    sage: y=function('y',x)
     489    sage: eq=diff(y,x)+y+x+1
     490    sage: nops(eq)
     491    4
     492
     493    AUTHORS:
     494
     495    - Yuri Karadzhov (2010-03-27)
     496
     497    """
     498    return len(ops(expr))
     499
     500
     501#Operators set for op function
     502
     503def identical_operator(x): return x
     504
     505def set_operator(x): return set(x)
     506
     507def tuple_operator(x): return tuple(x)
     508
     509def dict_operator(x): return dict(x)
     510
     511def matrix_operator(x): return matrix(x)
     512
     513def vector_operator(x): return vector(x)
     514
     515def diff_operator(x): return diff(x[0],x[1:])
     516
     517
     518def op(expr):
     519    """
     520    RETURN an expression operator
     521
     522    INPUT:
     523
     524    - ``expr`` - an expression
     525
     526    OUTPUT:
     527
     528    an expression operator
     529
     530    EXAMPLES:
     531
     532    op(expr) is the function that
     533    op(expr)(ops(expr)) equal to expr
     534
     535    unless
     536    sage: x,t=var(x,t)
     537    sage: f=function('f',x,t)
     538    sage: f
     539    f(x,t)
     540    sage: ops(f)
     541    [x,t]
     542    sage: op(f)(op(f))
     543    f(x,t)
     544    sage: f is f
     545    True
     546    sage: op(f)(ops(f)) is f
     547    False
     548
     549    but anyway
     550    sage: op(f)(ops(f)) - f
     551    0
     552
     553    AUTHORS:
     554
     555    - Yuri Karadzhov (2010-03-27)
     556
     557    """
     558    if isinstance(expr,sage.symbolic.expression.Expression):#get rid of lambda functions if it is possible
     559        t=expr.operator()
     560        if isinstance(t,sage.symbolic.operators.FDerivativeOperator):
     561            return diff_operator
     562        if not t is None: return lambda x: apply(t,x)
     563        else: return lambda x: expr
     564    else:
     565        if isinstance(expr,list):
     566            return identical_operator
     567        if isinstance(expr,set):
     568            return set_operator
     569        if isinstance(expr,tuple):
     570            return tuple_operator
     571        if isinstance(expr,dict):
     572            return dict_operator
     573        if isinstance(expr,sage.matrix.matrix_integer_dense.Matrix_integer_dense):
     574            return matrix_operator
     575        if isinstance(expr,sage.matrix.matrix_rational_dense.Matrix_rational_dense):
     576            return matrix_operator
     577        if isinstance(expr,sage.matrix.matrix_generic_dense.Matrix_generic_dense):
     578            return matrix_operator
     579        if isinstance(expr,sage.modules.vector_integer_dense.Vector_integer_dense):
     580            return vector_operator
     581        if isinstance(expr,sage.modules.vector_rational_dense.Vector_rational_dense):
     582            return vector_operator
     583        if isinstance(expr,sage.modules.free_module_element.FreeModuleElement_generic_dense):
     584            return vector_operator
     585        else: return lambda x: expr
     586
     587
     588#operator=op
     589
     590
     591def indets(expr,tp):
     592    """
     593    find indeterminates of given types
     594
     595    INPUT:
     596
     597    - ``expr`` - an expression
     598    - ``tp`` - symbolic type or types
     599
     600    OUTPUT:
     601
     602    set of indeterminants of given types
     603
     604    .. WARNING::
     605
     606    indets convert mutable types to apropriate immutable ones
     607    to avoid performance issues
     608
     609    EXAMPLES:
     610
     611    sage: x,a,b,c,t=var('x,a,b,c,t')
     612    sage: y=function('y',t,x)
     613    sage: eq=a*(diff(y,t,x))^2+b*diff(y,x)+3+x+ln(x+1)+c*x
     614    sage: indets(eq,'diff')
     615    set([D[0, 1](y)(t, x), D[1](y)(t, x)])
     616    sage: indets(eq,'_function')
     617    set([log(x + 1), y(t, x)])
     618    sage: indets(eq,'_function diff')
     619    set([D[0, 1](y)(t, x), D[1](y)(t, x), log(x + 1), y(t, x)])
     620
     621    sage: x,a,b,c,t=var('x,a,b,c,t')
     622    sage: g=function('g',x)
     623    sage: f=function('f',t,x)
     624    sage: eq=a*x+3+b*sin(x)+ln(f+x)+tan(1/x)+g+exp(c+t)
     625    sage: eq
     626    a*x + b*sin(x) + e^(c + t) + log(x + f(t, x)) + tan(1/x) + g(x) + 3
     627    sage: indets(eq,'_function')
     628    set([sin(x), log(x + f(t, x)), f(t, x), g(x), tan(1/x), e^(c + t)])
     629    sage: indets(eq,'_standard.function')
     630    set([sin(x), tan(1/x), log(x + f(t, x)), e^(c + t)])
     631    sage: indets(eq,'_trig.standard.function')
     632    set([sin(x), tan(1/x)])
     633    sage: indets(eq,'_trig')
     634    set([sin(x), tan(1/x)])
     635    sage: indets(eq,'_logarithmic.standard.function')
     636    set([log(x + f(t, x)), e^(c + t)])
     637    sage: indets(eq,'_logarithmic')
     638    set([log(x + f(t, x)), e^(c + t)])
     639    sage: indets(eq,'log.logarithmic.standard.function')
     640    set([log(x + f(t, x))])
     641    sage: indets(eq,'_log')
     642    set([log(x + f(t, x)), e^(c + t)])
     643    sage: indets(eq,'symbolic.function')
     644    set([f(t, x), g(x)])
     645    sage: indets(eq,'symbolic.variable')
     646    set([c, t, a, x, b])
     647    sage: indets(eq,'_symbolic')
     648    set([f(t, x), g(x), b, c, t, a, x])
     649
     650    sage: eq=log(x+1)+dilog(x+3)
     651    sage: indets(eq,'_log')
     652    set([dilog(x + 3), log(x + 1)])
     653    sage: indets(eq,'log.logarithmic.standard.function')
     654    set([log(x + 1)])
     655
     656    sage: x=var('x')
     657    sage: eq=(x+1)^2+ln(x+1)+3+x
     658    sage:indets(eq,'_+')
     659    set([x+1,(x+1)^2+ln(x+1)+3+x])
     660
     661    sage: eq=x + log((x + 5)/7) + 5
     662    sage: eq
     663    x + log(1/7*x + 5/7) + 5
     664    sage: indets(eq,'rational.number')
     665    set([5/7, 1/7])
     666    sage: indets(eq,'integer.number')
     667    set([5])
     668    sage: indets(eq,'_number')
     669    set([5/7, 5, 1/7])
     670
     671    AUTHORS:
     672
     673    - Yuri Karadzhov (2010-03-27)
     674
     675    """
     676    iset=set()
     677    if mtype(expr,tp):
     678        #Convert mutable to immutable ones
     679        if mtype(expr,'_.vector'):
     680            print 'Warning: vector {0} was converted to tuple'.format(expr)
     681            expr=tuple(expr)
     682        if mtype(expr,'set.collection'):
     683            print 'Warning: set {0} was converted to frozenset'.format(expr)
     684            expr=frozenset(expr)
     685        iset.add(expr)
     686    for el in ops(expr):
     687        iset=iset.union(indets(el,tp))
     688    return iset