Ticket #6479: trac_6479_marik_revised_2.patch

File trac_6479_marik_revised_2.patch, 26.1 KB (added by robert.marik, 12 years ago)

Apply on the top of the patch trac_6479_marik_revised.patch and on the top pf the patch for Ticket #385 http://trac.sagemath.org/sage_trac/ticket/385

  • sage/calculus/all.py

    # HG changeset patch
    # User Robert Marik <marik@mendelu.cz>
    # Date 1255441029 -7200
    # Node ID 4b7bf272a1f094e98358c1fc1bcf9b0383b80ad7
    # Parent  8da994a4ee55a3736f814f7b3bedf5283342c7bf
    Some enhancements to differential equations
    
    diff -r 8da994a4ee55 -r 4b7bf272a1f0 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 8da994a4ee55 -r 4b7bf272a1f0 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
     
    3843from sage.interfaces.maxima import MaximaElement, Maxima
    3944from sage.plot.all import line
    4045from sage.symbolic.expression import is_SymbolicEquation
     46from sage.symbolic.ring import is_SymbolicVariable
     47from sage.calculus.functional import diff
    4148
    4249maxima = Maxima()
    4350
    4451def desolve(de, dvar, ics=None, ivar=None, show_method=False, contrib_ode=False):
    4552    """
    4653    Solves a 1st or 2nd order linear ODE via maxima. Including IVP and BVP.
    47     In most cases returns SymbolicEquation which defines the solution implicitly.
    48     If the result is in the form y(x)=... (happens for linear eqs.),
    49     returns the right-hand side only.
    50     The possible constant solutions of separable ODE's are omitted.
    5154
    5255    INPUT:
    5356        de    -- an expression or equation representing the ODE
     
    6972                           of the equation which is separable but this property is not recognized
    7073                           by Maxima and equation is solved as exact.
    7174        contrib_ode  -- (optional) if true, desolve allows to solve clairot, lagrange, riccati and
    72                            some other equations. May take a long time and thus turned off by default. Initial conditions cqan be used only if the result is one SymbolicEquation
     75                           some other equations. May take a long time and thus turned off by default.
     76                           Initial conditions can be used only if the result is one SymbolicEquation
    7377                           (does not contain singular solution, for example)
    74                      
     78
     79
     80    OUTPUT:
     81        In most cases returns SymbolicEquation which defines the solution implicitly.
     82        If the result is in the form y(x)=... (happens for linear eqs.),
     83        returns the right-hand side only.
     84        The possible constant solutions of separable ODE's are omitted.
     85
    7586
    7687    EXAMPLES:
    7788        sage: x = var('x')
     
    104115        sage: desolve(y*diff(y,x)+sin(x)==0,y)
    105116        -1/2*y(x)^2 == c - cos(x)
    106117
    107     ### Produces error (difficult ODE)
    108     #sage: desolve(sqrt(y)*diff(y,x)+e^(y)+cos(x)-sin(x+y)==0,y)
    109 
    110     ### Produces error (difficult ODE) - moreover, takes a long time
    111     #sage: desolve(sqrt(y)*diff(y,x)+e^(y)+cos(x)-sin(x+y)==0,y,contrib_ode=True)
    112 
    113         sage: desolve(x*diff(y,x)^2-(1+x*y)*diff(y,x)+y==0,y,contrib_ode=True)
    114         [y(x) == c + log(x), y(x) == c*e^x]
    115         sage: desolve(x*diff(y,x)^2-(1+x*y)*diff(y,x)+y==0,y,contrib_ode=True,show_method=True)
    116         [[y(x) == c + log(x), y(x) == c*e^x], 'factor']
    117         sage: desolve(diff(y,x)^2+x*diff(y,x)-y==0,y,contrib_ode=True)
    118         [y(x) == c^2 + c*x, y(x) == -1/4*x^2]
     118    Clairot equation: general and sungular solutions
    119119        sage: desolve(diff(y,x)^2+x*diff(y,x)-y==0,y,contrib_ode=True,show_method=True)
    120120        [[y(x) == c^2 + c*x, y(x) == -1/4*x^2], 'clairault']
    121         sage: desolve(diff(y,x)==(x+y)^2,y,contrib_ode=True)
    122         [[x == c - arctan(sqrt(t)), y(x) == -x - sqrt(t)], [x == c + arctan(sqrt(t)), y(x) == -x + sqrt(t)]]
    123         sage: desolve(diff(y,x)==(x+y)^2,y,contrib_ode=True,show_method=True)
    124         [[[x == c - arctan(sqrt(t)), y(x) == -x - sqrt(t)], [x == c + arctan(sqrt(t)), y(x) == -x + sqrt(t)]], 'lagrange']
    125 
     121 
    126122    For equations involving more variables we specify independent variable
    127123        sage: a,b,c,n=var('a b c n')
    128124        sage: desolve(x^2*diff(y,x)==a+b*x^n+c*x^2*y^2,y,ivar=x,contrib_ode=True)
     
    139135        sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,1,3],show_method=True)
    140136        [1/6*y(x)^3 - 5/3*y(x) == x - 3/2, 'freeofx']
    141137
    142     ### These two examples produce error (as expected, Sage cannot solve equations from initial conditions)
    143     # sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,2]).expand()
    144     # sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,2],show_method=True)
    145 
    146 
    147138    Separable equations - Sage returns solution in implicit form
    148139
    149140        sage: desolve(diff(y,x)*sin(y) == cos(x),y)
     
    152143        [-cos(y(x)) == c + sin(x), 'separable']
    153144        sage: desolve(diff(y,x)*sin(y) == cos(x),y,[pi/2,1])
    154145        -cos(y(x)) == sin(x) - cos(1) - 1
    155         sage: desolve(diff(y,x)*sin(y) == cos(x),y,[pi/2,1],show_method=True)
    156         [-cos(y(x)) == sin(x) - cos(1) - 1, 'separable']
    157         sage: desolve(diff(y,x)*(y) == cos(x),y)
    158         1/2*y(x)^2 == c + sin(x)
    159146
    160147    Linear equation - Sage returns the expression on the right hand side only
    161148        sage: desolve(diff(y,x)+(y) == cos(x),y)
     
    164151        [1/2*((sin(x) + cos(x))*e^x + 2*c)*e^(-x), 'linear']
    165152        sage: desolve(diff(y,x)+(y) == cos(x),y,[0,1])
    166153        1/2*(e^x*sin(x) + e^x*cos(x) + 1)*e^(-x)
    167         sage: desolve(diff(y,x)+(y) == cos(x),y,[0,1],show_method=True)
    168         [1/2*(e^x*sin(x) + e^x*cos(x) + 1)*e^(-x), 'linear']
     154
     155    This ODE with separated variables is solved as exact
     156    Explanation: factor does not split exp(x-y) in maxima into exp(x)*exp(-y)
     157        sage: desolve(diff(y,x)==exp(x-y),y,show_method=True)
     158        [-e^x + e^y(x) == c, 'exact']
     159
     160    You can solve Bessel equations. You can also use initial
     161    conditions, but you cannot put (sometimes desired) initial
     162    condition at x=0, since this point is singlar point of the
     163    eruation. Anyway, if the solution should be bounded at x=0, then
     164    k2=0.
     165           
     166        sage: desolve(x^2*diff(y,x,x)+x*diff(y,x)+(x^2-4)*y==0,y)
     167        k1*bessel_j(2, x) + k2*bessel_y(2, x)
     168   
     169    This example asks (silly) question, whether sqrt(3) is an integer.
     170    has been reported to Maxima developers.
     171        # sage: desolve(x^2*diff(y,x,x)+x*diff(y,x)+(x^2-3)*y==0,y)
     172
     173
     174    Produces error (difficult ODE)
     175        #sage: desolve(sqrt(y)*diff(y,x)+e^(y)+cos(x)-sin(x+y)==0,y)
     176        #Traceback (click to the left for traceback)
     177        #...
     178        #NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True."
     179
     180        #Produces error (difficult ODE) - moreover, takes a long time
     181        #sage: desolve(sqrt(y)*diff(y,x)+e^(y)+cos(x)-sin(x+y)==0,y,contrib_ode=True)
     182
     183        sage: desolve(x*diff(y,x)^2-(1+x*y)*diff(y,x)+y==0,y,contrib_ode=True,show_method=True)
     184        [[y(x) == c + log(x), y(x) == c*e^x], 'factor']
     185
     186        sage: desolve(diff(y,x)==(x+y)^2,y,contrib_ode=True,show_method=True)
     187        [[[x == c - arctan(sqrt(t)), y(x) == -x - sqrt(t)], [x == c + arctan(sqrt(t)), y(x) == -x + sqrt(t)]], 'lagrange']
     188
     189    These two examples produce error (as expected, Maxima 5.18 cannot solve equations from initial conditions)
     190    Current Maxima 5.18 returns false answer in this case!
     191        #sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,2]).expand()
     192        #Traceback (click to the left for traceback)
     193        #...
     194        #NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True."
     195        #sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,2],show_method=True)
     196        #Traceback (click to the left for traceback)
     197        #...
     198        #NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True."
     199
    169200
    170201    Second order linear
    171202        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y)
     
    194225        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,[0,3,pi/2,2],show_method=True)
    195226        [(2*(2*e^(1/2*pi) - 3)*x/pi + 3)*e^(-x), 'constcoeff']
    196227
    197     This ODE with separated variables is solved as exact
    198     Explanation: factor does not split exp(x-y) in maxima into exp(x)*exp(-y)
    199         sage: desolve(diff(y,x)==exp(x-y),y,show_method=True)
    200         [-e^x + e^y(x) == c, 'exact']
    201 
    202228       
    203229    This equation can be solved within Maxima but not within Sage
    204230    needs assumptions assume(x>0,y>0) and works in maxima, but not in Sage
    205     # sage: assume(x>0)
    206     # sage: assume(y>0)
    207     # sage: maxima("assume(x>0,y>0)")
    208     # sage: desolve(x*diff(y,x)-x*sqrt(y^2+x^2)-y,y,show_method=True)   
     231        # sage: assume(x>0)
     232        # sage: assume(y>0)
     233        # sage: desolve(x*diff(y,x)-x*sqrt(y^2+x^2)-y,y,show_method=True)   
    209234
    210235    This equation can be solved within Maxima but not within Sage
    211     # sage: desolve(x*(x+1)*diff(y,x,2)+(x+5)*diff(y,x,1)+(-4)*y,y,contrib_ode=True)
     236    Perhaps upgrading Maxima help
     237        # sage: desolve(x*(x+1)*diff(y,x,2)+(x+5)*diff(y,x,1)+(-4)*y,y,contrib_ode=True)
     238
     239       
     240    TESTS:
     241    Trac #6479 fixed
     242        sage: x = var('x')
     243        sage: y = function('y', x)
     244        sage: desolve( diff(y,x,x) == 0, y, [0,0,1])
     245        x
     246        sage: desolve( diff(y,x,x) == 0, y, [0,1,1])
     247        x + 1
     248
    212249
    213250    AUTHOR: David Joyner (1-2006)
    214251            Robert Bradshaw (10-2008)
    215     CHANGES: Robert Marik (10-2009)
     252            Robert Marik (10-2009)
    216253
    217     Remark: The method contrib_ode has not been tested with initial conditions,
    218             since I was not able to find equation which cannot be solved using
    219             ode2, can be solved using contrib_ode and the result is sufficiently
    220             simple to continue safely (one SymbolicEquation) - robert.marik
    221            
    222254    """
    223255    if is_SymbolicEquation(de):
    224256        de = de.lhs() - de.rhs()
     257    if is_SymbolicVariable(dvar):
     258        raise ValueError, "You have to declare dependent variable as a function, eg. y=function('y',x)"
    225259    # for backwards compatibility
    226260    if isinstance(dvar, list):
    227261        dvar, ivar = dvar
     
    272306            soln=maxima(cmd)
    273307        if len(ics) == 3:
    274308            #fixed ic2 command from Maxima - we have to ensure that %k1, %k2 do not depend on variables, should be removed when fixed in Maxima
    275             maxima("ic2_sage(soln,xa,ya,dya):=block([programmode:true,backsubst:true,singsolve:true,temp,%k2,%k1], noteqn(xa), noteqn(ya), noteqn(dya), boundtest('%k1,%k1), boundtest('%k2,%k2), temp: lhs(soln) - rhs(soln), 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]), if not freeof(lhs(ya),TEMP_k) or not freeof(lhs(xa),TEMP_k) then return (false), temp: maplist(lambda([zz], subst(zz,soln)), TEMP_k), if length(temp)=1 then return(first(temp)) else return(temp))")
     309            maxima("ic2_sage(soln,xa,ya,dya):=block([programmode:true,backsubst:true,singsolve:true,temp,%k2,%k1,TEMP_k], \
     310                noteqn(xa), noteqn(ya), noteqn(dya), boundtest('%k1,%k1), boundtest('%k2,%k2), \
     311                temp: lhs(soln) - rhs(soln), \
     312                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]), \
     313                if not freeof(lhs(ya),TEMP_k) or not freeof(lhs(xa),TEMP_k) then return (false), \
     314                temp: maplist(lambda([zz], subst(zz,soln)), TEMP_k), \
     315                if length(temp)=1 then return(first(temp)) else return(temp))")
    276316            tempic=(ivar==ics[0])._maxima_().str()
    277317            tempic=tempic+","+(dvar==ics[1])._maxima_().str()
    278318            tempic=tempic+",'diff("+dvar_str+","+ivar_str+")="+(ics[2])._maxima_().str()
     
    285325                raise NotImplementedError, "Maxima was unable to solve this IVP. Remove the initial contition to get the general solution."   
    286326        if len(ics) == 4:
    287327            #fixed bc2 command from Maxima - we have to ensure that %k1, %k2 do not depend on variables, should be removed when fixed in Maxima
    288             maxima("bc2_sage(soln,xa,ya,xb,yb):=block([programmode:true,backsubst:true,singsolve:true,temp,%k1,%k2], noteqn(xa), noteqn(ya), noteqn(xb), noteqn(yb), boundtest('%k1,%k1), boundtest('%k2,%k2),TEMP_k:solve([subst([xa,ya],soln), subst([xb,yb],soln)], [%k1,%k2]),if not freeof(lhs(ya),TEMP_k) or not freeof(lhs(xa),TEMP_k) then return (false),temp: maplist(lambda([zz], subst(zz,soln)),TEMP_k),if length(temp)=1 then return(first(temp)) else return(temp))")
     328            maxima("bc2_sage(soln,xa,ya,xb,yb):=block([programmode:true,backsubst:true,singsolve:true,temp,%k1,%k2,TEMP_k], \
     329                noteqn(xa), noteqn(ya), noteqn(xb), noteqn(yb), boundtest('%k1,%k1), boundtest('%k2,%k2), \
     330                TEMP_k:solve([subst([xa,ya],soln), subst([xb,yb],soln)], [%k1,%k2]), \
     331                if not freeof(lhs(ya),TEMP_k) or not freeof(lhs(xa),TEMP_k) then return (false), \
     332                temp: maplist(lambda([zz], subst(zz,soln)),TEMP_k), \
     333                if length(temp)=1 then return(first(temp)) else return(temp))")
    289334            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())
    290335            cmd="(TEMP:%s,substitute(%s=%s(%s),TEMP))"%(cmd,dvar_str,dvar_str,ivar_str)           
    291336            cmd=sanitize_var(cmd)
     
    305350    else:   
    306351        return soln
    307352
     353
    308354#def desolve_laplace2(de,vars,ics=None):
    309355##     """
    310356##     Solves an ODE using laplace transforms via maxima. Initial conditions
     
    346392#    #return maxima.eval(cmd)
    347393#    return de0.desolve(vars[0]).rhs()
    348394   
    349 def desolve_laplace(de,vars,ics=None):
     395def desolve_laplace(de, dvar, ics=None, ivar=None):
     396
    350397    """
    351     Solves an ODE using laplace transforms via maxima. Initials conditions
    352     are optional.
     398    Solves an ODE using laplace transforms. Initials conditions are optional.
    353399
    354400    INPUT:
    355401        de    -- a lambda expression representing the ODE
    356                  (eg, de = "diff(f(x),x,2)=diff(f(x),x)+sin(x)")
    357         vars  -- a list of strings representing the variables
    358                  (eg, vars = ["x","f"], if x is the independent
    359                   variable and f is the dependent variable)
     402                 (eg, de = diff(y,x,2) == diff(y,x)+sin(x))
     403        dvar  -- the dependent variable (eg y)
     404        ivar  -- (optional) the independent variable  (hereafter called x),
     405                 which must be specified if there is more than one
     406                 independent variable in the equation.
    360407        ics   -- a list of numbers representing initial conditions,
    361                  with symbols allowed which are represented by strings
    362408                 (eg, f(0)=1, f'(0)=2 is ics = [0,1,2])
    363409
    364410    EXAMPLES:
    365         sage: from sage.calculus.desolvers import desolve_laplace
     411        sage: u=function('u',x)
     412        sage: eq = diff(u,x) - exp(-x) - u == 0
     413        sage: desolve_laplace(eq,u)
     414        1/2*(2*u(0) + 1)*e^x - 1/2*e^(-x)
     415   
     416    We can use initial conditions
     417        sage: desolve_laplace(eq,u,ics=[0,3])
     418        -1/2*e^(-x) + 7/2*e^x
     419
     420    The initial conditions do not persist in the system (as they persisted
     421    in previous versions)
     422        sage: desolve_laplace(eq,u)
     423        1/2*(2*u(0) + 1)*e^x - 1/2*e^(-x)
     424
     425        sage: f=function('f', x)
     426        sage: eq = diff(f,x) + f == 0
     427        sage: desolve_laplace(eq,f,[0,1])
     428        e^(-x)
     429
    366430        sage: x = var('x')
    367         sage: function('f', x)
    368         f(x)
    369         sage: de = lambda y: diff(y,x,x) - 2*diff(y,x) + y
    370         sage: desolve_laplace(de(f(x)),["x","f"])
    371         "x*%e^x*('at('diff('f(x),x,1),x=0))-'f(0)*x*%e^x+'f(0)*%e^x"
    372         sage: desolve_laplace(de(f(x)),["x","f"],[0,1,2])
    373          'x*%e^x+%e^x'
    374 
    375     WARNING:
    376         The second Sage command in the above example sets the values of f(0) and f'(0)
    377         in Maxima, so subsequent ODEs involving these variables will have these initial conditions
    378         automatically imposed.
     431        sage: f = function('f', x)
     432        sage: de = diff(f,x,x) - 2*diff(f,x) + f
     433        sage: desolve_laplace(de,f)
     434        -x*e^x*f(0) + x*e^x*D[0](f)(0) + e^x*f(0)
     435        sage: desolve_laplace(de,f,ics=[0,1,2])
     436        x*e^x + e^x
    379437
    380438    AUTHOR: David Joyner (1-2006,8-2007)
     439            Robert Marik (10-2009)
    381440    """
     441    #This is the original code from David Joyner (inputs and outputs strings)
     442    #maxima("de:"+de._repr_()+"=0;")
     443    #if ics!=None:
     444    #    d = len(ics)
     445    #    for i in range(0,d-1):
     446    #        ic = "atvalue(diff("+vars[1]+"("+vars[0]+"),"+str(vars[0])+","+str(i)+"),"+str(vars[0])+"="+str(ics[0])+","+str(ics[1+i])+")"
     447    #        maxima(ic)
     448    #
     449    #cmd = "desolve("+de._repr_()+","+vars[1]+"("+vars[0]+"));"
     450    #return maxima(cmd).rhs()._maxima_init_()
     451
     452    ## verbatim copy from desolve - begin
     453    if is_SymbolicEquation(de):
     454        de = de.lhs() - de.rhs()
     455    if is_SymbolicVariable(dvar):
     456        raise ValueError, "You have to declare dependent variable as a function, eg. y=function('y',x)"
     457    # for backwards compatibility
     458    if isinstance(dvar, list):
     459        dvar, ivar = dvar
     460    elif ivar is None:
     461        ivars = de.variables()
     462        ivars = [t for t in ivars if t != dvar]
     463        if len(ivars) != 1:
     464            raise ValueError, "Unable to determine independent variable, please specify."
     465        ivar = ivars[0]
     466    ## verbatim copy from desolve - end
     467
     468    def sanitize_var(exprs):  # 'y(x) -> y(x)
     469        return exprs.replace("'"+str(dvar),str(dvar))   
     470    de0=de._maxima_()
     471    cmd = sanitize_var("desolve("+de0.str()+","+str(dvar)+")")
     472    soln=maxima(cmd).rhs()
     473    if str(soln).strip() == 'false':
     474        raise NotImplementedError, "Maxima was unable to solve this ODE."   
     475    soln=soln.sage()
    382476    if ics!=None:
    383477        d = len(ics)
    384478        for i in range(0,d-1):
    385             ic = "atvalue(diff('"+vars[1]+"("+vars[0]+"),"+str(vars[0])+","+str(i)+"),"+str(vars[0])+"="+str(ics[0])+","+str(ics[1+i])+")"
    386             maxima(ic)
    387             #print i,ic
    388     cmd = "desolve("+de._maxima_init_()+",'"+vars[1]+"("+vars[0]+"));"
    389     return maxima(cmd).rhs()._maxima_init_()
     479            soln=eval('soln.substitute(diff(dvar,ivar,i)('+str(ivar)+'=ics[0])==ics[i+1])')
     480    return soln
     481   
    390482   
    391483def desolve_system(des, vars, ics=None, ivar=None):
    392484    """
    393     Solves any size system of 1st order odes using maxima. Initials conditions
    394     are optional.
     485    Solves any size system of 1st order ODE's. Initials conditions are optional.
    395486
    396487    INPUT:
    397488        des   -- list of ODEs
     
    449540
    450541def desolve_system_strings(des,vars,ics=None):
    451542    """
    452     Solves any size system of 1st order odes using maxima. Initials conditions
    453     are optional.
     543    Solves any size system of 1st order ODE's. Initials conditions are optional.
    454544
    455545
    456546    INPUT:
     
    686776    Q1 = line([[x[0],x[1]] for x in soln], rgbcolor=(1/4,1/8,3/4))
    687777    Q2 = line([[x[0],x[2]] for x in soln], rgbcolor=(1/2,1/8,1/4))
    688778    return [Q1,Q2]
     779
     780
     781def desolve_rk4(de, dvar, ics=None, ivar=None, end_point=10, step=0.1, output='list', method='maxima',**kwds):
     782    """
     783    Solves numerically one first-order ordinary differential equation.
     784
     785    INPUT:
     786    input is similar to desolve command. The differential equation can be written in
     787    a form close to the plot_slope_field or desolve command
     788
     789    Variant 1
     790        de - right hand side, i.e. the function f(x,y) from ODE y'=f(x,y)
     791        dvar - dependent variable (symbolic variable declared by var)
     792
     793    Variant 2
     794        de - equation, including term with diff(y,x)
     795        dvar - dependent variable (declared as funciton of independent variable)
     796
     797    Other parameters
     798
     799        ivar - should be specified, if there are more variables or if the equation is autonomous
     800        ics  - initial conditions in the form [x0,y0]
     801        end_point  - the end point of the interval, if smaller than x0 we extend the solution
     802               to the left
     803        step - the length of the step
     804        output - 'list', 'plot', 'slope_field' (graph of the solution with slope field)
     805
     806
     807    OUTPUT:
     808    Returns a list of points, or plot produced by listplot, optionally with slope field.
     809
     810
     811    EXAMPLES:
     812    sage: from sage.calculus.desolvers import desolve_rk4
     813    sage: x,y=var('x y')
     814    sage: P=desolve_rk4(y*(2-y),y,ics=[0,.1],ivar=x,output='slope_field',end_point=6,thickness=3)
     815
     816    sage: x,y=var('x y')
     817    sage: P=desolve_rk4(y*(2-y),y,ics=[0,.5],ivar=x,output='plot',end_point=6,thickness=3)
     818    sage: PP=desolve_rk4(y*(2-y),y,ics=[0,.5],ivar=x,output='plot',end_point=-6,thickness=3)
     819    sage: PPP=plot_slope_field(y*(2-y),(x,PP.xmin(),P.xmax()),(y,PP.ymin(),P.ymax()))
     820    sage: Q=P+PP+PPP
     821
     822    ALGORITHM:
     823    4th order Runge-Kutta method. Wrapper for command rk in Maxima's dynamics package.
     824    Perhaps could be faster by using fast_float instead.
     825
     826    AUTHOR: Robert Marik (10-2009)
     827    """
     828    if ics is None:
     829        raise ValueError, "No initial conditions, specify with ics=[x0,y0]."
     830
     831    if ivar is None:
     832        ivars = de.variables()
     833        ivars = [t for t in ivars if t != dvar]
     834        if len(ivars) != 1:
     835            raise ValueError, "Unable to determine independent variable, please specify."
     836        ivar = ivars[0]
     837
     838    if not is_SymbolicVariable(dvar):
     839        from sage.calculus.var import var
     840        from sage.calculus.all import diff
     841        from sage.symbolic.relation import solve
     842        if is_SymbolicEquation(de):
     843            de = de.lhs() - de.rhs()
     844        dummy_dvar=var('dummy_dvar')
     845        # consider to add warning if the solution is not unique
     846        de=solve(de,diff(dvar,ivar),solution_dict=True)
     847        if len(de) != 1:
     848            raise NotImplementedError, "Sorry, cannot find explicit formula for right-hand side of the ODE."
     849        de=de[0][diff(dvar,ivar)].subs(dvar==dummy_dvar)
     850    else:
     851        dummy_dvar=dvar
     852
     853    if end_point<ics[0]:
     854        step = -abs(step)
     855
     856    de0=de._maxima_().str()
     857    cmd="rk(%s,%s,%s,[%s,%s,%s,%s])\
     858    "%(de0,str(dummy_dvar),str(ics[1]),str(ivar),str(ics[0]),end_point,step)
     859    maxima("load('dynamics)")
     860    sol=maxima(cmd).sage()
     861
     862    if output=='list':
     863        return sol
     864    from sage.plot.plot import list_plot
     865    from sage.plot.plot_field import plot_slope_field
     866    R = list_plot(sol,plotjoined=True,**kwds)
     867    if output=='plot':
     868        return R
     869    if output=='slope_field':
     870        XMIN=sol[0][0]
     871        YMIN=sol[0][1]
     872        XMAX=XMIN
     873        YMAX=YMIN
     874        for s,t in sol:
     875            if s>XMAX:XMAX=s
     876            if s<XMIN:XMIN=s
     877            if t>YMAX:YMAX=t
     878            if t<YMIN:YMIN=t
     879        return plot_slope_field(de,(ivar,XMIN,XMAX),(dummy_dvar,YMIN,YMAX))+R
     880
     881    raise ValueError, "Option output should be 'list', 'plot' or 'slope_field'."
     882
     883
     884def desolve_system_rk4(des, vars, ics=None, ivar=None, end_point=10, step=0.1):
     885    """
     886    Solves numerically system of first-order ordinary differential
     887    equations using the 4th order Runge-Kutta method. Wrapper for
     888    Maxima command rk. 
     889
     890    INPUT:
     891    input is similar to desolve_system and desolve_rk4 commands
     892        des - right hand sides of the system
     893        vars - dependent variables
     894        ivar - should be specified, if there are more variables or if the equation is autonomous\
     895               and the independent variable is missing
     896        ics  - initial conditions in the form [x0,y01,y02,y03,....]
     897        end_point  - the end point of the interval, upper bound for x
     898        step - the length of the step
     899
     900
     901    OUTPUT:
     902    Returns a list of points.
     903
     904
     905    EXAMPLES:
     906    Lotka Volterra system
     907
     908    sage: from sage.calculus.desolvers import desolve_system_rk4
     909    sage: x,y,t=var('x y t')
     910    sage: P=desolve_system_rk4([x*(1-y),-y*(1-x)],[x,y],ics=[0,0.5,2],ivar=t,end_point=20)
     911    sage: Q=[ [i,j] for i,j,k in P]
     912    sage: LP=list_plot(Q)
     913
     914    sage: Q=[ [j,k] for i,j,k in P]
     915    sage: LP=list_plot(Q)
     916   
     917    ALGORITHM:
     918    4th order Runge-Kutta method. Wrapper for command rk in Maxima's dynamics package.
     919    Perhaps could be faster by using fast_float instead.
     920
     921    AUTHOR: Robert Marik (10-2009)
     922    """
     923
     924    if ics is None:
     925        raise ValueError, "No initial conditions, specify with ics=[x0,y01,y02,...]."
     926
     927    ivars = set([])
     928
     929    for de in des:
     930        ivars = ivars.union(set(de.variables()))
     931    if ivar is None:
     932        ivars = ivars - set(vars)
     933        if len(ivars) != 1:
     934            raise ValueError, "Unable to determine independent variable, please specify."
     935        ivar = list(ivars)[0]
     936
     937    dess = [de._maxima_().str() for de in des]
     938    desstr = "[" + ",".join(dess) + "]"
     939    varss = [varsi._maxima_().str() for varsi in vars]
     940    varstr = "[" + ",".join(varss) + "]"
     941    x0=ics[0]
     942    icss = [ics[i]._maxima_().str() for i in range(1,len(ics))]
     943    icstr = "[" + ",".join(icss) + "]"
     944    if end_point<x0:
     945        step=-abs(step)
     946    cmd="rk(%s,%s,%s,[%s,%s,%s,%s])\
     947    "%(desstr,varstr,icstr,str(ivar),str(x0),end_point,step)
     948    maxima("load('dynamics)")
     949    sol=maxima(cmd)
     950
     951    return sol.sage()