Ticket #8616: trac_8616_symbolic_sage.2.patch

File trac_8616_symbolic_sage.2.patch, 31.8 KB (added by yuri.k, 11 years ago)

symbolic sage module

  • sage/calculus/all.py

    # HG changeset patch
    # User Yuri Karadzhov <yuri.karadzhov@gmail.com>
    # Date 1269721789 0
    # Node ID acadcfa52c47da58a72fe835959b7f429e9298af
    # Parent  46b1ee5997ac7e7faafdd86751fe7d121b6c1385
    Trac 8616: Add functionality. 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.
    
    diff -r 46b1ee5997ac -r acadcfa52c47 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 acadcfa52c47 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 acadcfa52c47 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 *