Ticket #6479: trac_6479_marik_revision3.patch

File trac_6479_marik_revision3.patch, 33.9 KB (added by robert.marik, 12 years ago)

this replaces previous patches and installs on the top of patch for trac #385

  • sage/calculus/all.py

    # HG changeset patch
    # User Robert Marik <marik@mendelu.cz>
    # Date 1256678359 -3600
    # Node ID 3fe9af889529137ee0fbb660a8e3c46ce7541ecd
    # Parent  164faaa675c4895b80270a29326faf398a6b0c79
    trac 6479
    
    diff -r 164faaa675c4 -r 3fe9af889529 sage/calculus/all.py
    a b  
    1010
    1111from desolvers import (desolve, desolve_laplace, desolve_system,
    1212                       eulers_method, eulers_method_2x2,
    13                        eulers_method_2x2_plot)
     13                       eulers_method_2x2_plot, desolve_rk4, desolve_system_rk4)
    1414
    1515from var import (var, function, clear_vars)
    1616
  • sage/calculus/desolvers.py

    diff -r 164faaa675c4 -r 3fe9af889529 sage/calculus/desolvers.py
    a b  
    1111* desolve_system -- Solves any size system of 1st order odes using Maxima.
    1212                    Initials conditions are optional.
    1313
     14* desolve_rk4 -- Solves numerically IVP for one first order equation,
     15                 returns list of points or plot
     16
     17* desolve_system_rk4 -- Solves numerically IVP for system of first order
     18                        equations, returns list of points
     19
    1420* eulers_method -- Approximate solution to a 1st order DE, presented as a table.
    1521
    1622* eulers_method_2x2 -- Approximate solution to a 1st order system of DEs,
     
    2228AUTHORS: David Joyner (3-2006)     -- Initial version of functions
    2329         Marshall Hampton (7-2007) -- Creation of Python module and testing
    2430         Robert Bradshaw (10-2008) -- Some interface cleanup.
    25 
    26 Other functions for solving DEs are given in functions/elementary.py.
     31         Robert Marik (10-2009) -- Some bugfixes and enhancements
    2732
    2833"""
    2934
    3035##########################################################################
    31 #  Copyright (C) 2006 David Joyner <wdjoyner@gmail.com>, Marshall Hampton
     36#  Copyright (C) 2006 David Joyner <wdjoyner@gmail.com>, Marshall Hampton,
     37#  Robert Marik <marik@mendelu.cz>
    3238#
    3339#  Distributed under the terms of the GNU General Public License (GPL):
    3440#
     
    3844from sage.interfaces.maxima import MaximaElement, Maxima
    3945from sage.plot.all import line
    4046from sage.symbolic.expression import is_SymbolicEquation
     47from sage.symbolic.ring import is_SymbolicVariable
     48from sage.calculus.functional import diff
    4149
    4250maxima = Maxima()
    4351
    44 def desolve(de, dvar, ics=None, ivar=None):
     52def desolve(de, dvar, ics=None, ivar=None, show_method=False, contrib_ode=False):
    4553    """
    46     Solves a 1st or 2nd order linear ODE via maxima. Initials conditions
    47     are not given.
     54    Solves a 1st or 2nd order linear ODE via maxima. Including IVP and BVP.
    4855
    4956    INPUT:
    5057        de    -- an expression or equation representing the ODE
     
    5360                    for a first-order equation, specify the initial x and y
    5461                    for a second-order equation, specify the initial x, y, and dy/dx
    5562                    for a second-order boundary solution, specify initial and final x and y
     63                    initial conditions gives error if the solution is not SymbolisEquation
     64                    (as happens for example for clairot equation)
    5665        ivar  -- (optional) the independent variable  (hereafter called x),
    5766                    which must be specified if there is more than one
    5867                    independent variable in the equation.
     68        show_method  -- (optional) if true, then Sage returns pair [solution, method], where method
     69                           is the string describing method which has been used to get solution
     70                           (Maxima uses the following order for first order equations: linear, separable,
     71                           exact (including exact with integrating factor), homogeneous, bernoulli,
     72                           generalized homogeneous) - use carefully in class, see below for the example
     73                           of the equation which is separable but this property is not recognized
     74                           by Maxima and equation is solved as exact.
     75        contrib_ode  -- (optional) if true, desolve allows to solve clairot, lagrange, riccati and
     76                           some other equations. May take a long time and thus turned off by default.
     77                           Initial conditions can be used only if the result is one SymbolicEquation
     78                           (does not contain singular solution, for example)
     79
     80
     81    OUTPUT:
     82        In most cases returns SymbolicEquation which defines the solution implicitly.
     83        If the result is in the form y(x)=... (happens for linear eqs.),
     84        returns the right-hand side only.
     85        The possible constant solutions of separable ODE's are omitted.
     86
    5987
    6088    EXAMPLES:
    6189        sage: x = var('x')
     
    6795        sage: plot(f)
    6896       
    6997    We can also solve second-order differential equations.
     98   
    7099        sage: x = var('x')
    71100        sage: y = function('y', x)
    72101        sage: de = diff(y,x,2) - y == x
    73102        sage: desolve(de, y)
    74103        k1*e^x + k2*e^(-x) - x
    75104        sage: f = desolve(de, y, [10,2,1]); f
    76         1/2*(y(10) + 12)*e^(x - 10) + 1/2*(e^10*y(10) + 8*e^10)*e^(-x) - x
    77         sage: f(x=10).expand().simplify()
    78         y(10)
    79         sage: diff(f,x)(x=10).expand().simplify()
     105        -x + 5*e^(-x + 10) + 7*e^(x - 10)
     106        sage: f(x=10)
     107        2
     108        sage: diff(f,x)(x=10)
    80109        1
    81110
     111        sage: de = diff(y,x,2) + y == 0
     112        sage: desolve(de, y)
     113        k1*sin(x) + k2*cos(x)
     114        sage: desolve(de, y, [0,1,pi/2,4])
     115        4*sin(x) + cos(x)
     116
     117        sage: desolve(y*diff(y,x)+sin(x)==0,y)
     118        -1/2*y(x)^2 == c - cos(x)
     119
     120    Clairot equation: general and sungular solutions
     121
     122        sage: desolve(diff(y,x)^2+x*diff(y,x)-y==0,y,contrib_ode=True,show_method=True)
     123        [[y(x) == c^2 + c*x, y(x) == -1/4*x^2], 'clairault']
     124 
     125    For equations involving more variables we specify independent variable
     126
     127        sage: a,b,c,n=var('a b c n')
     128        sage: desolve(x^2*diff(y,x)==a+b*x^n+c*x^2*y^2,y,ivar=x,contrib_ode=True)
     129        [[y(x) == 0, (b*x^(n - 2) + a/x^2)*c^2*u == 0]]
     130        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)
     131        [[[y(x) == 0, (b*x^(n - 2) + a/x^2)*c^2*u == 0]], 'riccati']
     132
     133
     134    Higher orded, not involving independent variable
     135
     136        sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y).expand()
     137        1/6*y(x)^3 + k1*y(x) == k2 + x
     138        sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,1,3]).expand()
     139        1/6*y(x)^3 - 5/3*y(x) == x - 3/2
     140        sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,1,3],show_method=True)
     141        [1/6*y(x)^3 - 5/3*y(x) == x - 3/2, 'freeofx']
     142
     143    Separable equations - Sage returns solution in implicit form
     144
     145        sage: desolve(diff(y,x)*sin(y) == cos(x),y)
     146        -cos(y(x)) == c + sin(x)
     147        sage: desolve(diff(y,x)*sin(y) == cos(x),y,show_method=True)
     148        [-cos(y(x)) == c + sin(x), 'separable']
     149        sage: desolve(diff(y,x)*sin(y) == cos(x),y,[pi/2,1])
     150        -cos(y(x)) == sin(x) - cos(1) - 1
     151
     152    Linear equation - Sage returns the expression on the right hand side only
     153
     154        sage: desolve(diff(y,x)+(y) == cos(x),y)
     155        1/2*((sin(x) + cos(x))*e^x + 2*c)*e^(-x)
     156        sage: desolve(diff(y,x)+(y) == cos(x),y,show_method=True)
     157        [1/2*((sin(x) + cos(x))*e^x + 2*c)*e^(-x), 'linear']
     158        sage: desolve(diff(y,x)+(y) == cos(x),y,[0,1])
     159        1/2*(e^x*sin(x) + e^x*cos(x) + 1)*e^(-x)
     160
     161    This ODE with separated variables is solved as exact
     162    Explanation - factor does not split exp(x-y) in maxima into exp(x)*exp(-y)
     163
     164        sage: desolve(diff(y,x)==exp(x-y),y,show_method=True)
     165        [-e^x + e^y(x) == c, 'exact']
     166
     167    You can solve Bessel equations. You can also use initial
     168    conditions, but you cannot put (sometimes desired) initial
     169    condition at x=0, since this point is singlar point of the
     170    eruation. Anyway, if the solution should be bounded at x=0, then
     171    k2=0.
     172           
     173        sage: desolve(x^2*diff(y,x,x)+x*diff(y,x)+(x^2-4)*y==0,y)
     174        k1*bessel_j(2, x) + k2*bessel_y(2, x)
     175   
     176    This example asks (silly) question, whether sqrt(3) is an integer.
     177    has been reported to Maxima developers.
     178        # sage: desolve(x^2*diff(y,x,x)+x*diff(y,x)+(x^2-3)*y==0,y)
     179
     180
     181    Produces error (difficult ODE)
     182        #sage: desolve(sqrt(y)*diff(y,x)+e^(y)+cos(x)-sin(x+y)==0,y)
     183        #Traceback (click to the left for traceback)
     184        #...
     185        #NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True."
     186
     187        #Produces error (difficult ODE) - moreover, takes a long time
     188        #sage: desolve(sqrt(y)*diff(y,x)+e^(y)+cos(x)-sin(x+y)==0,y,contrib_ode=True)
     189
     190    Some more types od ODE's
     191
     192        sage: desolve(x*diff(y,x)^2-(1+x*y)*diff(y,x)+y==0,y,contrib_ode=True,show_method=True)
     193        [[y(x) == c + log(x), y(x) == c*e^x], 'factor']
     194
     195        sage: desolve(diff(y,x)==(x+y)^2,y,contrib_ode=True,show_method=True)
     196        [[[x == c - arctan(sqrt(t)), y(x) == -x - sqrt(t)], [x == c + arctan(sqrt(t)), y(x) == -x + sqrt(t)]], 'lagrange']
     197
     198    These two examples produce error (as expected, Maxima 5.18 cannot solve equations from initial conditions)
     199    Current Maxima 5.18 returns false answer in this case!
     200
     201        #sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,2]).expand()
     202        #Traceback (click to the left for traceback)
     203        #...
     204        #NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True."
     205        #sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,2],show_method=True)
     206        #Traceback (click to the left for traceback)
     207        #...
     208        #NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True."
     209
     210
     211    Second order linear ODE
     212
     213        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y)
     214        (k2*x + k1)*e^(-x) + 1/2*sin(x)
     215        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,show_method=True)
     216        [(k2*x + k1)*e^(-x) + 1/2*sin(x), 'variationofparameters']
     217        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,[0,3,1])
     218        1/2*(7*x + 6)*e^(-x) + 1/2*sin(x)
     219        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,[0,3,1],show_method=True)
     220        [1/2*(7*x + 6)*e^(-x) + 1/2*sin(x), 'variationofparameters']
     221        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,[0,3,pi/2,2])
     222        3*((e^(1/2*pi) - 2)*x/pi + 1)*e^(-x) + 1/2*sin(x)
     223        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,[0,3,pi/2,2],show_method=True)
     224        [3*((e^(1/2*pi) - 2)*x/pi + 1)*e^(-x) + 1/2*sin(x), 'variationofparameters']
     225
     226        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y)
     227        (k2*x + k1)*e^(-x)
     228        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,show_method=True)
     229        [(k2*x + k1)*e^(-x), 'constcoeff']
     230        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,[0,3,1])
     231        (4*x + 3)*e^(-x)
     232        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,[0,3,1],show_method=True)
     233        [(4*x + 3)*e^(-x), 'constcoeff']
     234        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,[0,3,pi/2,2])
     235        (2*(2*e^(1/2*pi) - 3)*x/pi + 3)*e^(-x)
     236        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,[0,3,pi/2,2],show_method=True)
     237        [(2*(2*e^(1/2*pi) - 3)*x/pi + 3)*e^(-x), 'constcoeff']
     238
     239       
     240    This equation can be solved within Maxima but not within Sage
     241    needs assumptions assume(x>0,y>0) and works in maxima, but not in Sage
     242        # sage: assume(x>0)
     243        # sage: assume(y>0)
     244        # sage: desolve(x*diff(y,x)-x*sqrt(y^2+x^2)-y,y,show_method=True)   
     245
     246    This equation can be solved within Maxima but not within Sage
     247    Perhaps upgrading Maxima help
     248        # sage: desolve(x*(x+1)*diff(y,x,2)+(x+5)*diff(y,x,1)+(-4)*y,y,contrib_ode=True)
     249
     250       
     251    TESTS:
     252   
     253    Trac #6479 fixed
     254
     255        sage: x = var('x')
     256        sage: y = function('y', x)
     257        sage: desolve( diff(y,x,x) == 0, y, [0,0,1])
     258        x
     259        sage: desolve( diff(y,x,x) == 0, y, [0,1,1])
     260        x + 1
     261
     262
    82263    AUTHOR: David Joyner (1-2006)
    83264            Robert Bradshaw (10-2008)
     265            Robert Marik (10-2009)
     266
    84267    """
    85268    if is_SymbolicEquation(de):
    86269        de = de.lhs() - de.rhs()
     270    if is_SymbolicVariable(dvar):
     271        raise ValueError, "You have to declare dependent variable as a function, eg. y=function('y',x)"
    87272    # for backwards compatibility
    88273    if isinstance(dvar, list):
    89274        dvar, ivar = dvar
     
    93278        if len(ivars) != 1:
    94279            raise ValueError, "Unable to determine independent variable, please specify."
    95280        ivar = ivars[0]
    96     de0 = de._maxima_()
    97     soln = de0.ode2(dvar, ivar)
     281    def sanitize_var(exprs):
     282        return exprs.replace("'"+dvar_str+"("+ivar_str+")",dvar_str)   
     283    dvar_str=dvar.operator()._maxima_().str()
     284    ivar_str=ivar._maxima_().str()
     285    de00 = de._maxima_().str()
     286    de0 = sanitize_var(de00)
     287    ode_solver="ode2"
     288    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)   
     289    # we produce string like this
     290    # ode2('diff(y,x,2)+2*'diff(y,x,1)+y-cos(x),y(x),x)
     291    soln = maxima(cmd)
     292
    98293    if str(soln).strip() == 'false':
    99         raise NotImplementedError, "Maxima was unable to solve this system."
    100     def to_eqns(lhs, exprs):
    101         eqns = []
    102         for lhs, expr in zip(lhs, exprs):
    103             if is_SymbolicEquation(expr):
    104                 eqns.append(expr)
    105             else:
    106                 if lhs == dvar and len(exprs) == 2:
    107                     ivar_ic = exprs[0] # figure this out...
    108                     lhs = lhs.subs(lhs.default_variable()==ivar_ic)
    109                 eqns.append(lhs == expr)
    110         return eqns
    111     if ics is not None:
     294        if contrib_ode:
     295            ode_solver="contrib_ode"
     296            maxima("load('contrib_ode)")
     297            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)   
     298            # we produce string like this
     299            # (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))
     300            soln = maxima(cmd)
     301            if str(soln).strip() == 'false':
     302                raise NotImplementedError, "Maxima was unable to solve this ODE."   
     303        else:
     304            raise NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True."
     305
     306    if show_method:
     307        maxima_method=maxima("method")
     308
     309    if (ics is not None):
     310        if not is_SymbolicEquation(soln.sage()):
     311             raise NotImplementedError, "Maxima was unable to use initial condition for this equation (%s)"%(maxima_method.str())   
    112312        if len(ics) == 2:
    113             ics = to_eqns([ivar, dvar], ics)
    114             soln = soln.ic1(ics[0], ics[1])
     313            tempic=(ivar==ics[0])._maxima_().str()
     314            tempic=tempic+","+(dvar==ics[1])._maxima_().str()
     315            cmd="(TEMP:ic1(%s(%s,%s,%s),%s),substitute(%s=%s(%s),TEMP))"%(ode_solver,de00,dvar_str,ivar_str,tempic,dvar_str,dvar_str,ivar_str)
     316            cmd=sanitize_var(cmd)
     317            # we produce string like this
     318            # (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))
     319            soln=maxima(cmd)
    115320        if len(ics) == 3:
    116             ics = to_eqns([ivar, dvar, dvar.derivative(ivar)], ics)
    117             soln = soln.ic2(*ics)
     321            #fixed ic2 command from Maxima - we have to ensure that %k1, %k2 do not depend on variables, should be removed when fixed in Maxima
     322            maxima("ic2_sage(soln,xa,ya,dya):=block([programmode:true,backsubst:true,singsolve:true,temp,%k2,%k1,TEMP_k], \
     323                noteqn(xa), noteqn(ya), noteqn(dya), boundtest('%k1,%k1), boundtest('%k2,%k2), \
     324                temp: lhs(soln) - rhs(soln), \
     325                TEMP_k:solve([subst([xa,ya],soln), subst([dya,xa], lhs(dya)=-subst(0,lhs(dya),diff(temp,lhs(xa)))/diff(temp,lhs(ya)))],[%k1,%k2]), \
     326                if not freeof(lhs(ya),TEMP_k) or not freeof(lhs(xa),TEMP_k) then return (false), \
     327                temp: maplist(lambda([zz], subst(zz,soln)), TEMP_k), \
     328                if length(temp)=1 then return(first(temp)) else return(temp))")
     329            tempic=(ivar==ics[0])._maxima_().str()
     330            tempic=tempic+","+(dvar==ics[1])._maxima_().str()
     331            tempic=tempic+",'diff("+dvar_str+","+ivar_str+")="+(ics[2])._maxima_().str()
     332            cmd="(TEMP:ic2_sage(%s(%s,%s,%s),%s),substitute(%s=%s(%s),TEMP))"%(ode_solver,de00,dvar_str,ivar_str,tempic,dvar_str,dvar_str,ivar_str)
     333            cmd=sanitize_var(cmd)
     334            # we produce string like this
     335            # (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))
     336            soln=maxima(cmd)
     337            if str(soln).strip() == 'false':
     338                raise NotImplementedError, "Maxima was unable to solve this IVP. Remove the initial contition to get the general solution."   
    118339        if len(ics) == 4:
    119             ics = to_eqns([ivar, dvar, ivar, dvar], ics)
    120             soln = soln.bc(*ics)
    121     if soln.lhs() == dvar:
     340            #fixed bc2 command from Maxima - we have to ensure that %k1, %k2 do not depend on variables, should be removed when fixed in Maxima
     341            maxima("bc2_sage(soln,xa,ya,xb,yb):=block([programmode:true,backsubst:true,singsolve:true,temp,%k1,%k2,TEMP_k], \
     342                noteqn(xa), noteqn(ya), noteqn(xb), noteqn(yb), boundtest('%k1,%k1), boundtest('%k2,%k2), \
     343                TEMP_k:solve([subst([xa,ya],soln), subst([xb,yb],soln)], [%k1,%k2]), \
     344                if not freeof(lhs(ya),TEMP_k) or not freeof(lhs(xa),TEMP_k) then return (false), \
     345                temp: maplist(lambda([zz], subst(zz,soln)),TEMP_k), \
     346                if length(temp)=1 then return(first(temp)) else return(temp))")
     347            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())
     348            cmd="(TEMP:%s,substitute(%s=%s(%s),TEMP))"%(cmd,dvar_str,dvar_str,ivar_str)           
     349            cmd=sanitize_var(cmd)
     350            # we produce string like this
     351            # (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))
     352            soln=maxima(cmd)
     353            if str(soln).strip() == 'false':
     354                raise NotImplementedError, "Maxima was unable to solve this BVP. Remove the initial contition to get the general solution."   
     355
     356    soln=soln.sage()
     357    if is_SymbolicEquation(soln) and soln.lhs() == dvar:
     358        # Remark: Here we do not check that the right hand side does not depend on dvar.
     359        # This probably will not hapen for soutions obtained via ode2, anyway.
    122360        soln = soln.rhs()
    123     return soln.sage()
     361    if show_method:
     362        return [soln,maxima_method.str()]
     363    else:   
     364        return soln
     365
    124366
    125367#def desolve_laplace2(de,vars,ics=None):
    126368##     """
     
    162404#    #cmd = "desolve("+de+","+vars[1]+"("+vars[0]+"));"
    163405#    #return maxima.eval(cmd)
    164406#    return de0.desolve(vars[0]).rhs()
     407
    165408   
    166 def desolve_laplace(de,vars,ics=None):
     409def desolve_laplace(de, dvar, ics=None, ivar=None):
    167410    """
    168     Solves an ODE using laplace transforms via maxima. Initials conditions
    169     are optional.
     411    Solves an ODE using laplace transforms. Initials conditions are optional.
    170412
    171413    INPUT:
    172414        de    -- a lambda expression representing the ODE
    173                  (eg, de = "diff(f(x),x,2)=diff(f(x),x)+sin(x)")
    174         vars  -- a list of strings representing the variables
    175                  (eg, vars = ["x","f"], if x is the independent
    176                   variable and f is the dependent variable)
     415                 (eg, de = diff(y,x,2) == diff(y,x)+sin(x))
     416        dvar  -- the dependent variable (eg y)
     417        ivar  -- (optional) the independent variable  (hereafter called x),
     418                 which must be specified if there is more than one
     419                 independent variable in the equation.
    177420        ics   -- a list of numbers representing initial conditions,
    178                  with symbols allowed which are represented by strings
    179421                 (eg, f(0)=1, f'(0)=2 is ics = [0,1,2])
    180422
     423    OUTPUT:
     424        Solution of the ODE as symbolic expression
     425
    181426    EXAMPLES:
    182         sage: from sage.calculus.desolvers import desolve_laplace
     427        sage: u=function('u',x)
     428        sage: eq = diff(u,x) - exp(-x) - u == 0
     429        sage: desolve_laplace(eq,u)
     430        1/2*(2*u(0) + 1)*e^x - 1/2*e^(-x)
     431   
     432    We can use initial conditions
     433
     434        sage: desolve_laplace(eq,u,ics=[0,3])
     435        -1/2*e^(-x) + 7/2*e^x
     436
     437    The initial conditions do not persist in the system (as they persisted
     438    in previous versions)
     439
     440        sage: desolve_laplace(eq,u)
     441        1/2*(2*u(0) + 1)*e^x - 1/2*e^(-x)
     442
     443        sage: f=function('f', x)
     444        sage: eq = diff(f,x) + f == 0
     445        sage: desolve_laplace(eq,f,[0,1])
     446        e^(-x)
     447
    183448        sage: x = var('x')
    184         sage: function('f', x)
    185         f(x)
    186         sage: de = lambda y: diff(y,x,x) - 2*diff(y,x) + y
    187         sage: desolve_laplace(de(f(x)),["x","f"])
    188         "x*%e^x*('at('diff('f(x),x,1),x=0))-'f(0)*x*%e^x+'f(0)*%e^x"
    189         sage: desolve_laplace(de(f(x)),["x","f"],[0,1,2])
    190          'x*%e^x+%e^x'
     449        sage: f = function('f', x)
     450        sage: de = diff(f,x,x) - 2*diff(f,x) + f
     451        sage: desolve_laplace(de,f)
     452        -x*e^x*f(0) + x*e^x*D[0](f)(0) + e^x*f(0)
     453        sage: desolve_laplace(de,f,ics=[0,1,2])
     454        x*e^x + e^x
    191455
    192     WARNING:
    193         The second Sage command in the above example sets the values of f(0) and f'(0)
    194         in Maxima, so subsequent ODEs involving these variables will have these initial conditions
    195         automatically imposed.
     456    TESTS:
     457        #4839 fixed
     458
     459        sage: t=var('t')
     460        sage: x=function('x', t)
     461        sage: soln=desolve_laplace(diff(x,t)+x==1, x, ics=[0,2])
     462        sage: soln
     463        e^(-t) + 1
     464        sage: soln(t=3)
     465        e^(-3) + 1
    196466
    197467    AUTHOR: David Joyner (1-2006,8-2007)
     468            Robert Marik (10-2009)
    198469    """
     470    #This is the original code from David Joyner (inputs and outputs strings)
     471    #maxima("de:"+de._repr_()+"=0;")
     472    #if ics!=None:
     473    #    d = len(ics)
     474    #    for i in range(0,d-1):
     475    #        ic = "atvalue(diff("+vars[1]+"("+vars[0]+"),"+str(vars[0])+","+str(i)+"),"+str(vars[0])+"="+str(ics[0])+","+str(ics[1+i])+")"
     476    #        maxima(ic)
     477    #
     478    #cmd = "desolve("+de._repr_()+","+vars[1]+"("+vars[0]+"));"
     479    #return maxima(cmd).rhs()._maxima_init_()
     480
     481    ## verbatim copy from desolve - begin
     482    if is_SymbolicEquation(de):
     483        de = de.lhs() - de.rhs()
     484    if is_SymbolicVariable(dvar):
     485        raise ValueError, "You have to declare dependent variable as a function, eg. y=function('y',x)"
     486    # for backwards compatibility
     487    if isinstance(dvar, list):
     488        dvar, ivar = dvar
     489    elif ivar is None:
     490        ivars = de.variables()
     491        ivars = [t for t in ivars if t != dvar]
     492        if len(ivars) != 1:
     493            raise ValueError, "Unable to determine independent variable, please specify."
     494        ivar = ivars[0]
     495    ## verbatim copy from desolve - end
     496
     497    def sanitize_var(exprs):  # 'y(x) -> y(x)
     498        return exprs.replace("'"+str(dvar),str(dvar))   
     499    de0=de._maxima_()
     500    cmd = sanitize_var("desolve("+de0.str()+","+str(dvar)+")")
     501    soln=maxima(cmd).rhs()
     502    if str(soln).strip() == 'false':
     503        raise NotImplementedError, "Maxima was unable to solve this ODE."   
     504    soln=soln.sage()
    199505    if ics!=None:
    200506        d = len(ics)
    201507        for i in range(0,d-1):
    202             ic = "atvalue(diff('"+vars[1]+"("+vars[0]+"),"+str(vars[0])+","+str(i)+"),"+str(vars[0])+"="+str(ics[0])+","+str(ics[1+i])+")"
    203             maxima(ic)
    204             #print i,ic
    205     cmd = "desolve("+de._maxima_init_()+",'"+vars[1]+"("+vars[0]+"));"
    206     return maxima(cmd).rhs()._maxima_init_()
     508            soln=eval('soln.substitute(diff(dvar,ivar,i)('+str(ivar)+'=ics[0])==ics[i+1])')
     509    return soln
     510   
    207511   
    208512def desolve_system(des, vars, ics=None, ivar=None):
    209513    """
    210     Solves any size system of 1st order odes using maxima. Initials conditions
    211     are optional.
     514    Solves any size system of 1st order ODE's. Initials conditions are optional.
    212515
    213516    INPUT:
    214517        des   -- list of ODEs
     
    266569
    267570def desolve_system_strings(des,vars,ics=None):
    268571    """
    269     Solves any size system of 1st order odes using maxima. Initials conditions
    270     are optional.
     572    Solves any size system of 1st order ODE's. Initials conditions are optional.
    271573
    272574
    273575    INPUT:
     
    503805    Q1 = line([[x[0],x[1]] for x in soln], rgbcolor=(1/4,1/8,3/4))
    504806    Q2 = line([[x[0],x[2]] for x in soln], rgbcolor=(1/2,1/8,1/4))
    505807    return [Q1,Q2]
     808
     809def desolve_rk4_determine_bounds(ics,end_points=None):
     810    """
     811    Used to determine bounds for numerical integration.
     812   
     813    If end_points is None, the interval for integration is from ics[0]
     814    to ics[0]+10
     815
     816    If end_points is a or [a], the interval for integration is from min(ics[0],a)
     817    to max(ics[0],a)
     818
     819    If end_points is [a,b], the interval for integration is from min(ics[0],a)
     820    to max(ics[0],b)
     821
     822   EXAMPLES:
     823        sage: from sage.calculus.desolvers import desolve_rk4_determine_bounds
     824        sage: desolve_rk4_determine_bounds([0,2],1)
     825        (0, 1)
     826        sage: desolve_rk4_determine_bounds([0,2]) 
     827        (0, 10)
     828        sage: desolve_rk4_determine_bounds([0,2],[-2])
     829        (-2, 0)
     830        sage: desolve_rk4_determine_bounds([0,2],[-2,4])
     831        (-2, 4)
     832   
     833    """
     834    if end_points is None:
     835        return((ics[0],ics[0]+10))
     836    if not isinstance(end_points,list):
     837        end_points=[end_points]
     838    if len(end_points)==1:
     839        return (min(ics[0],end_points[0]),max(ics[0],end_points[0]))
     840    else:
     841        return (min(ics[0],end_points[0]),max(ics[0],end_points[1]))
     842   
     843
     844def desolve_rk4(de, dvar, ics=None, ivar=None, end_points=None, step=0.1, output='list', **kwds):
     845    """
     846    Solves numerically one first-order ordinary differential equation.
     847
     848    INPUT:
     849        input is similar to desolve command. The differential equation can be
     850        written in a form close to the plot_slope_field or desolve command
     851
     852        Variant 1 (function in two variables)
     853        de   -- right hand side, i.e. the function f(x,y) from ODE y'=f(x,y)
     854        dvar -- dependent variable (symbolic variable declared by var)
     855
     856        Variant 2 (symbolic equation)
     857        de   -- equation, including term with diff(y,x)
     858        dvar -- dependent variable (declared as funciton of independent variable)
     859
     860        Other parameters
     861        ivar  -- should be specified, if there are more variables or if the equation is autonomous
     862        ics   -- initial conditions in the form [x0,y0]
     863        end_points  -- the end points of the interval
     864           if end_points is a or [a], we integrate on between min(ics[0],a) and max(ics[0],a)
     865           if end_points is None, we use end_points=ics[0]+10
     866           if end_points is [a,b] we integrate on between min(ics[0],a) and max(ics[0],b)
     867        step   -- the length of the step (positive number)
     868        output -- 'list', 'plot', 'slope_field' (graph of the solution with slope field)
     869
     870
     871    OUTPUT:
     872        Returns a list of points, or plot produced by list_plot, optionally with slope field.
     873
     874
     875    EXAMPLES:
     876        sage: from sage.calculus.desolvers import desolve_rk4
     877
     878    Variant 2 for input - more common in numerics
     879
     880        sage: x,y=var('x y')
     881        sage: desolve_rk4(x*y*(2-y),y,ics=[0,1],end_points=1,step=0.5)
     882        [[0, 1], [0.5, 1.12419127425], [1.0, 1.46159016229]]
     883
     884    Variant 1 for input - we can pass ODE in the form used by
     885    desolve function In this example we integrate bakwards, since
     886    end_points < ics[0]
     887
     888        sage: y=function('y',x)
     889        sage: desolve_rk4(diff(y,x)+y*(y-1) == x-2,y,ics=[1,1],step=0.5, end_points=0)
     890        [[0.0, 8.90425710896], [0.5, 1.90932794536], [1, 1]]
     891
     892    Here we show how to plot simple pictures. For more advanced aplications use
     893    list_plot instead. To see the resulting picture use show(P) in Sage notebook.
     894
     895        sage: x,y=var('x y')
     896        sage: P=desolve_rk4(y*(2-y),y,ics=[0,.1],ivar=x,output='slope_field',end_points=[-4,6],thickness=3)
     897               
     898    ALGORITHM:
     899        4th order Runge-Kutta method. Wrapper for command rk in Maxima's dynamics package.
     900        Perhaps could be faster by using fast_float instead.
     901
     902    AUTHOR:
     903        Robert Marik (10-2009)
     904    """
     905    if ics is None:
     906        raise ValueError, "No initial conditions, specify with ics=[x0,y0]."
     907
     908    if ivar is None:
     909        ivars = de.variables()
     910        ivars = [t for t in ivars if t != dvar]
     911        if len(ivars) != 1:
     912            raise ValueError, "Unable to determine independent variable, please specify."
     913        ivar = ivars[0]
     914
     915    if not is_SymbolicVariable(dvar):
     916        from sage.calculus.var import var
     917        from sage.calculus.all import diff
     918        from sage.symbolic.relation import solve
     919        if is_SymbolicEquation(de):
     920            de = de.lhs() - de.rhs()
     921        dummy_dvar=var('dummy_dvar')
     922        # consider to add warning if the solution is not unique
     923        de=solve(de,diff(dvar,ivar),solution_dict=True)
     924        if len(de) != 1:
     925            raise NotImplementedError, "Sorry, cannot find explicit formula for right-hand side of the ODE."
     926        de=de[0][diff(dvar,ivar)].subs(dvar==dummy_dvar)
     927    else:
     928        dummy_dvar=dvar
     929
     930    step=abs(step)
     931    de0=de._maxima_()
     932    maxima("load('dynamics)")
     933    lower_bound,upper_bound=desolve_rk4_determine_bounds(ics,end_points)
     934    sol_1, sol_2 = [],[]
     935    if lower_bound<ics[0]:
     936        cmd="rk(%s,%s,%s,[%s,%s,%s,%s])\
     937        "%(de0.str(),str(dummy_dvar),str(ics[1]),str(ivar),str(ics[0]),lower_bound,-step)
     938        sol_1=maxima(cmd).sage()
     939        sol_1.pop(0)
     940        sol_1.reverse()
     941    if upper_bound>ics[0]:
     942        cmd="rk(%s,%s,%s,[%s,%s,%s,%s])\
     943        "%(de0.str(),str(dummy_dvar),str(ics[1]),str(ivar),str(ics[0]),upper_bound,step)
     944        sol_2=maxima(cmd).sage()
     945        sol_2.pop(0)
     946    sol=sol_1
     947    sol.extend([[ics[0],ics[1]]])
     948    sol.extend(sol_2)
     949
     950    if output=='list':
     951        return sol
     952    from sage.plot.plot import list_plot
     953    from sage.plot.plot_field import plot_slope_field
     954    R = list_plot(sol,plotjoined=True,**kwds)
     955    if output=='plot':
     956        return R
     957    if output=='slope_field':
     958        XMIN=sol[0][0]
     959        YMIN=sol[0][1]
     960        XMAX=XMIN
     961        YMAX=YMIN
     962        for s,t in sol:
     963            if s>XMAX:XMAX=s
     964            if s<XMIN:XMIN=s
     965            if t>YMAX:YMAX=t
     966            if t<YMIN:YMIN=t
     967        return plot_slope_field(de,(ivar,XMIN,XMAX),(dummy_dvar,YMIN,YMAX))+R
     968
     969    raise ValueError, "Option output should be 'list', 'plot' or 'slope_field'."
     970
     971
     972def desolve_system_rk4(des, vars, ics=None, ivar=None, end_points=None, step=0.1):
     973    """
     974    Solves numerically system of first-order ordinary differential
     975    equations using the 4th order Runge-Kutta method. Wrapper for
     976    Maxima command rk. 
     977
     978    INPUT:
     979        input is similar to desolve_system and desolve_rk4 commands
     980        des  -- right hand sides of the system
     981        vars -- dependent variables
     982        ivar -- should be specified, if there are more variables or if the equation is autonomous
     983               and the independent variable is missing
     984        ics  -- initial conditions in the form [x0,y01,y02,y03,....]
     985        end_points  -- the end points of the interval
     986           if end_points is a or [a], we integrate on between min(ics[0],a) and max(ics[0],a)
     987           if end_points is None, we use end_points=ics[0]+10
     988           if end_points is [a,b] we integrate on between min(ics[0],a) and max(ics[0],b)
     989        step -- the length of the step
     990
     991
     992    OUTPUT:
     993        Returns a list of points.
     994
     995
     996    EXAMPLES:
     997    Lotka Volterra system
     998
     999        sage: from sage.calculus.desolvers import desolve_system_rk4
     1000        sage: x,y,t=var('x y t')
     1001        sage: P=desolve_system_rk4([x*(1-y),-y*(1-x)],[x,y],ics=[0,0.5,2],ivar=t,end_points=20)
     1002        sage: Q=[ [i,j] for i,j,k in P]
     1003        sage: LP=list_plot(Q)
     1004       
     1005        sage: Q=[ [j,k] for i,j,k in P]
     1006        sage: LP=list_plot(Q)
     1007   
     1008    ALGORITHM:
     1009        4th order Runge-Kutta method. Wrapper for command rk in Maxima's dynamics package.
     1010        Perhaps could be faster by using fast_float instead.
     1011
     1012    AUTHOR: Robert Marik (10-2009)
     1013    """
     1014
     1015    if ics is None:
     1016        raise ValueError, "No initial conditions, specify with ics=[x0,y01,y02,...]."
     1017
     1018    ivars = set([])
     1019
     1020    for de in des:
     1021        ivars = ivars.union(set(de.variables()))
     1022    if ivar is None:
     1023        ivars = ivars - set(vars)
     1024        if len(ivars) != 1:
     1025            raise ValueError, "Unable to determine independent variable, please specify."
     1026        ivar = list(ivars)[0]
     1027
     1028    dess = [de._maxima_().str() for de in des]
     1029    desstr = "[" + ",".join(dess) + "]"
     1030    varss = [varsi._maxima_().str() for varsi in vars]
     1031    varstr = "[" + ",".join(varss) + "]"
     1032    x0=ics[0]
     1033    icss = [ics[i]._maxima_().str() for i in range(1,len(ics))]
     1034    icstr = "[" + ",".join(icss) + "]"
     1035    step=abs(step)
     1036
     1037    maxima("load('dynamics)")
     1038    lower_bound,upper_bound=desolve_rk4_determine_bounds(ics,end_points)
     1039    sol_1, sol_2 = [],[]
     1040    if lower_bound<ics[0]:
     1041        cmd="rk(%s,%s,%s,[%s,%s,%s,%s])\
     1042        "%(desstr,varstr,icstr,str(ivar),str(x0),lower_bound,-step)
     1043        sol_1=maxima(cmd).sage()
     1044        sol_1.pop(0)
     1045        sol_1.reverse()
     1046    if upper_bound>ics[0]:
     1047        cmd="rk(%s,%s,%s,[%s,%s,%s,%s])\
     1048        "%(desstr,varstr,icstr,str(ivar),str(x0),upper_bound,step)
     1049        sol_2=maxima(cmd).sage()
     1050        sol_2.pop(0)
     1051    sol=sol_1
     1052    sol.append(ics)
     1053    sol.extend(sol_2)
     1054
     1055    return sol