Ticket #6479: trac_6479_marik_revised.patch

File trac_6479_marik_revised.patch, 16.0 KB (added by robert.marik, 12 years ago)

Apply only this patch

  • sage/calculus/desolvers.py

    # HG changeset patch
    # User Robert Marik <marik@mendelu.cz>
    # Date 1254921466 -7200
    # Node ID 227631aed73eb8d55656213f88e4496384063baf
    # Parent  e334f15125a3faecdecfefc6079a7ba4536955e8
    trac #6479 plus some enhancements
    
    diff -r e334f15125a3 -r 227631aed73e sage/calculus/desolvers.py
    a b  
    4141
    4242maxima = Maxima()
    4343
    44 def desolve(de, dvar, ics=None, ivar=None):
     44def desolve(de, dvar, ics=None, ivar=None, show_method=False, contrib_ode=False):
    4545    """
    46     Solves a 1st or 2nd order linear ODE via maxima. Initials conditions
    47     are not given.
     46    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.
    4851
    4952    INPUT:
    5053        de    -- an expression or equation representing the ODE
     
    5356                    for a first-order equation, specify the initial x and y
    5457                    for a second-order equation, specify the initial x, y, and dy/dx
    5558                    for a second-order boundary solution, specify initial and final x and y
     59                    initial conditions gives error if the solution is not SymbolisEquation
     60                    (as happens for example for clairot equation)
    5661        ivar  -- (optional) the independent variable  (hereafter called x),
    5762                    which must be specified if there is more than one
    5863                    independent variable in the equation.
     64        show_method  -- (optional) if true, then Sage returns pair [solution, method], where method
     65                           is the string describing method which has been used to get solution
     66                           (Maxima uses the following order for first order equations: linear, separable,
     67                           exact (including exact with integrating factor), homogeneous, bernoulli,
     68                           generalized homogeneous) - use carefully in class, see below for the example
     69                           of the equation which is separable but this property is not recognized
     70                           by Maxima and equation is solved as exact.
     71        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
     73                           (does not contain singular solution, for example)
     74                     
    5975
    6076    EXAMPLES:
    6177        sage: x = var('x')
     
    7389        sage: desolve(de, y)
    7490        k1*e^x + k2*e^(-x) - x
    7591        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()
     92        -x + 5*e^(-x + 10) + 7*e^(x - 10)
     93        sage: f(x=10)
     94        2
     95        sage: diff(f,x)(x=10)
    8096        1
    8197
     98        sage: de = diff(y,x,2) + y == 0
     99        sage: desolve(de, y)
     100        k1*sin(x) + k2*cos(x)
     101        sage: desolve(de, y, [0,1,pi/2,4])
     102        4*sin(x) + cos(x)
     103
     104        sage: desolve(y*diff(y,x)+sin(x)==0,y)
     105        -1/2*y(x)^2 == c - cos(x)
     106
     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]
     119        sage: desolve(diff(y,x)^2+x*diff(y,x)-y==0,y,contrib_ode=True,show_method=True)
     120        [[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
     126    For equations involving more variables we specify independent variable
     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        sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y).expand()
     136        1/6*y(x)^3 + k1*y(x) == k2 + x
     137        sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,1,3]).expand()
     138        1/6*y(x)^3 - 5/3*y(x) == x - 3/2
     139        sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,1,3],show_method=True)
     140        [1/6*y(x)^3 - 5/3*y(x) == x - 3/2, 'freeofx']
     141
     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
     147    Separable equations - Sage returns solution in implicit form
     148
     149        sage: desolve(diff(y,x)*sin(y) == cos(x),y)
     150        -cos(y(x)) == c + sin(x)
     151        sage: desolve(diff(y,x)*sin(y) == cos(x),y,show_method=True)
     152        [-cos(y(x)) == c + sin(x), 'separable']
     153        sage: desolve(diff(y,x)*sin(y) == cos(x),y,[pi/2,1])
     154        -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)
     159
     160    Linear equation - Sage returns the expression on the right hand side only
     161        sage: desolve(diff(y,x)+(y) == cos(x),y)
     162        1/2*((sin(x) + cos(x))*e^x + 2*c)*e^(-x)
     163        sage: desolve(diff(y,x)+(y) == cos(x),y,show_method=True)
     164        [1/2*((sin(x) + cos(x))*e^x + 2*c)*e^(-x), 'linear']
     165        sage: desolve(diff(y,x)+(y) == cos(x),y,[0,1])
     166        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']
     169
     170    Second order linear
     171        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y)
     172        (k2*x + k1)*e^(-x) + 1/2*sin(x)
     173        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,show_method=True)
     174        [(k2*x + k1)*e^(-x) + 1/2*sin(x), 'variationofparameters']
     175        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,[0,3,1])
     176        1/2*(7*x + 6)*e^(-x) + 1/2*sin(x)
     177        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,[0,3,1],show_method=True)
     178        [1/2*(7*x + 6)*e^(-x) + 1/2*sin(x), 'variationofparameters']
     179        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,[0,3,pi/2,2])
     180        3*((e^(1/2*pi) - 2)*x/pi + 1)*e^(-x) + 1/2*sin(x)
     181        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,[0,3,pi/2,2],show_method=True)
     182        [3*((e^(1/2*pi) - 2)*x/pi + 1)*e^(-x) + 1/2*sin(x), 'variationofparameters']
     183
     184        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y)
     185        (k2*x + k1)*e^(-x)
     186        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,show_method=True)
     187        [(k2*x + k1)*e^(-x), 'constcoeff']
     188        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,[0,3,1])
     189        (4*x + 3)*e^(-x)
     190        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,[0,3,1],show_method=True)
     191        [(4*x + 3)*e^(-x), 'constcoeff']
     192        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,[0,3,pi/2,2])
     193        (2*(2*e^(1/2*pi) - 3)*x/pi + 3)*e^(-x)
     194        sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,[0,3,pi/2,2],show_method=True)
     195        [(2*(2*e^(1/2*pi) - 3)*x/pi + 3)*e^(-x), 'constcoeff']
     196
     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
     202       
     203    This equation can be solved within Maxima but not within Sage
     204    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)   
     209
     210    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)
     212
    82213    AUTHOR: David Joyner (1-2006)
    83214            Robert Bradshaw (10-2008)
     215    CHANGES: Robert Marik (10-2009)
     216
     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           
    84222    """
    85223    if is_SymbolicEquation(de):
    86224        de = de.lhs() - de.rhs()
     
    93231        if len(ivars) != 1:
    94232            raise ValueError, "Unable to determine independent variable, please specify."
    95233        ivar = ivars[0]
    96     de0 = de._maxima_()
    97     soln = de0.ode2(dvar, ivar)
     234    def sanitize_var(exprs):
     235        return exprs.replace("'"+dvar_str+"("+ivar_str+")",dvar_str)   
     236    dvar_str=dvar.operator()._maxima_().str()
     237    ivar_str=ivar._maxima_().str()
     238    de00 = de._maxima_().str()
     239    de0 = sanitize_var(de00)
     240    ode_solver="ode2"
     241    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)   
     242    # we produce string like this
     243    # ode2('diff(y,x,2)+2*'diff(y,x,1)+y-cos(x),y(x),x)
     244    soln = maxima(cmd)
     245
    98246    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:
     247        if contrib_ode:
     248            ode_solver="contrib_ode"
     249            maxima("load('contrib_ode)")
     250            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)   
     251            # we produce string like this
     252            # (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))
     253            soln = maxima(cmd)
     254            if str(soln).strip() == 'false':
     255                raise NotImplementedError, "Maxima was unable to solve this ODE."   
     256        else:
     257            raise NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True."
     258
     259    if show_method:
     260        maxima_method=maxima("method")
     261
     262    if (ics is not None):
     263        if not is_SymbolicEquation(soln.sage()):
     264             raise NotImplementedError, "Maxima was unable to use initial condition for this equation (%s)"%(maxima_method.str())   
    112265        if len(ics) == 2:
    113             ics = to_eqns([ivar, dvar], ics)
    114             soln = soln.ic1(ics[0], ics[1])
     266            tempic=(ivar==ics[0])._maxima_().str()
     267            tempic=tempic+","+(dvar==ics[1])._maxima_().str()
     268            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)
     269            cmd=sanitize_var(cmd)
     270            # we produce string like this
     271            # (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))
     272            soln=maxima(cmd)
    115273        if len(ics) == 3:
    116             ics = to_eqns([ivar, dvar, dvar.derivative(ivar)], ics)
    117             soln = soln.ic2(*ics)
     274            #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))")
     276            tempic=(ivar==ics[0])._maxima_().str()
     277            tempic=tempic+","+(dvar==ics[1])._maxima_().str()
     278            tempic=tempic+",'diff("+dvar_str+","+ivar_str+")="+(ics[2])._maxima_().str()
     279            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)
     280            cmd=sanitize_var(cmd)
     281            # we produce string like this
     282            # (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))
     283            soln=maxima(cmd)
     284            if str(soln).strip() == 'false':
     285                raise NotImplementedError, "Maxima was unable to solve this IVP. Remove the initial contition to get the general solution."   
    118286        if len(ics) == 4:
    119             ics = to_eqns([ivar, dvar, ivar, dvar], ics)
    120             soln = soln.bc(*ics)
    121     if soln.lhs() == dvar:
     287            #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))")
     289            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())
     290            cmd="(TEMP:%s,substitute(%s=%s(%s),TEMP))"%(cmd,dvar_str,dvar_str,ivar_str)           
     291            cmd=sanitize_var(cmd)
     292            # we produce string like this
     293            # (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))
     294            soln=maxima(cmd)
     295            if str(soln).strip() == 'false':
     296                raise NotImplementedError, "Maxima was unable to solve this BVP. Remove the initial contition to get the general solution."   
     297
     298    soln=soln.sage()
     299    if is_SymbolicEquation(soln) and soln.lhs() == dvar:
     300        # Remark: Here we do not check that the right hand side does not depend on dvar.
     301        # This probably will not hapen for soutions obtained via ode2, anyway.
    122302        soln = soln.rhs()
    123     return soln.sage()
     303    if show_method:
     304        return [soln,maxima_method.str()]
     305    else:   
     306        return soln
    124307
    125308#def desolve_laplace2(de,vars,ics=None):
    126309##     """