# 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 maxima = Maxima() def desolve(de, dvar, ics=None, ivar=None): def desolve(de, dvar, ics=None, ivar=None, show_method=False, contrib_ode=False): """ Solves a 1st or 2nd order linear ODE via maxima. Initials conditions are not given. Solves a 1st or 2nd order linear ODE via maxima. Including IVP and BVP. In most cases returns SymbolicEquation which defines the solution implicitly. If the result is in the form y(x)=... (happens for linear eqs.), returns the right-hand side only. The possible constant solutions of separable ODE's are omitted. INPUT: de    -- an expression or equation representing the ODE for a first-order equation, specify the initial x and y for a second-order equation, specify the initial x, y, and dy/dx for a second-order boundary solution, specify initial and final x and y initial conditions gives error if the solution is not SymbolisEquation (as happens for example for clairot equation) ivar  -- (optional) the independent variable  (hereafter called x), which must be specified if there is more than one independent variable in the equation. show_method  -- (optional) if true, then Sage returns pair [solution, method], where method is the string describing method which has been used to get solution (Maxima uses the following order for first order equations: linear, separable, exact (including exact with integrating factor), homogeneous, bernoulli, generalized homogeneous) - use carefully in class, see below for the example of the equation which is separable but this property is not recognized by Maxima and equation is solved as exact. contrib_ode  -- (optional) if true, desolve allows to solve clairot, lagrange, riccati and 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 (does not contain singular solution, for example) EXAMPLES: sage: x = var('x') sage: desolve(de, y) k1*e^x + k2*e^(-x) - x sage: f = desolve(de, y, [10,2,1]); f 1/2*(y(10) + 12)*e^(x - 10) + 1/2*(e^10*y(10) + 8*e^10)*e^(-x) - x sage: f(x=10).expand().simplify() y(10) sage: diff(f,x)(x=10).expand().simplify() -x + 5*e^(-x + 10) + 7*e^(x - 10) sage: f(x=10) 2 sage: diff(f,x)(x=10) 1 sage: de = diff(y,x,2) + y == 0 sage: desolve(de, y) k1*sin(x) + k2*cos(x) sage: desolve(de, y, [0,1,pi/2,4]) 4*sin(x) + cos(x) sage: desolve(y*diff(y,x)+sin(x)==0,y) -1/2*y(x)^2 == c - cos(x) ### Produces error (difficult ODE) #sage: desolve(sqrt(y)*diff(y,x)+e^(y)+cos(x)-sin(x+y)==0,y) ### Produces error (difficult ODE) - moreover, takes a long time #sage: desolve(sqrt(y)*diff(y,x)+e^(y)+cos(x)-sin(x+y)==0,y,contrib_ode=True) sage: desolve(x*diff(y,x)^2-(1+x*y)*diff(y,x)+y==0,y,contrib_ode=True) [y(x) == c + log(x), y(x) == c*e^x] sage: desolve(x*diff(y,x)^2-(1+x*y)*diff(y,x)+y==0,y,contrib_ode=True,show_method=True) [[y(x) == c + log(x), y(x) == c*e^x], 'factor'] sage: desolve(diff(y,x)^2+x*diff(y,x)-y==0,y,contrib_ode=True) [y(x) == c^2 + c*x, y(x) == -1/4*x^2] sage: desolve(diff(y,x)^2+x*diff(y,x)-y==0,y,contrib_ode=True,show_method=True) [[y(x) == c^2 + c*x, y(x) == -1/4*x^2], 'clairault'] sage: desolve(diff(y,x)==(x+y)^2,y,contrib_ode=True) [[x == c - arctan(sqrt(t)), y(x) == -x - sqrt(t)], [x == c + arctan(sqrt(t)), y(x) == -x + sqrt(t)]] sage: desolve(diff(y,x)==(x+y)^2,y,contrib_ode=True,show_method=True) [[[x == c - arctan(sqrt(t)), y(x) == -x - sqrt(t)], [x == c + arctan(sqrt(t)), y(x) == -x + sqrt(t)]], 'lagrange'] For equations involving more variables we specify independent variable sage: a,b,c,n=var('a b c n') sage: desolve(x^2*diff(y,x)==a+b*x^n+c*x^2*y^2,y,ivar=x,contrib_ode=True) [[y(x) == 0, (b*x^(n - 2) + a/x^2)*c^2*u == 0]] 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) [[[y(x) == 0, (b*x^(n - 2) + a/x^2)*c^2*u == 0]], 'riccati'] Higher orded, not involving independent variable sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y).expand() 1/6*y(x)^3 + k1*y(x) == k2 + x sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,1,3]).expand() 1/6*y(x)^3 - 5/3*y(x) == x - 3/2 sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,1,3],show_method=True) [1/6*y(x)^3 - 5/3*y(x) == x - 3/2, 'freeofx'] ### These two examples produce error (as expected, Sage cannot solve equations from initial conditions) # sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,2]).expand() # sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,2],show_method=True) Separable equations - Sage returns solution in implicit form sage: desolve(diff(y,x)*sin(y) == cos(x),y) -cos(y(x)) == c + sin(x) sage: desolve(diff(y,x)*sin(y) == cos(x),y,show_method=True) [-cos(y(x)) == c + sin(x), 'separable'] sage: desolve(diff(y,x)*sin(y) == cos(x),y,[pi/2,1]) -cos(y(x)) == sin(x) - cos(1) - 1 sage: desolve(diff(y,x)*sin(y) == cos(x),y,[pi/2,1],show_method=True) [-cos(y(x)) == sin(x) - cos(1) - 1, 'separable'] sage: desolve(diff(y,x)*(y) == cos(x),y) 1/2*y(x)^2 == c + sin(x) Linear equation - Sage returns the expression on the right hand side only sage: desolve(diff(y,x)+(y) == cos(x),y) 1/2*((sin(x) + cos(x))*e^x + 2*c)*e^(-x) sage: desolve(diff(y,x)+(y) == cos(x),y,show_method=True) [1/2*((sin(x) + cos(x))*e^x + 2*c)*e^(-x), 'linear'] sage: desolve(diff(y,x)+(y) == cos(x),y,[0,1]) 1/2*(e^x*sin(x) + e^x*cos(x) + 1)*e^(-x) sage: desolve(diff(y,x)+(y) == cos(x),y,[0,1],show_method=True) [1/2*(e^x*sin(x) + e^x*cos(x) + 1)*e^(-x), 'linear'] Second order linear sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y) (k2*x + k1)*e^(-x) + 1/2*sin(x) sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,show_method=True) [(k2*x + k1)*e^(-x) + 1/2*sin(x), 'variationofparameters'] sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,[0,3,1]) 1/2*(7*x + 6)*e^(-x) + 1/2*sin(x) sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,[0,3,1],show_method=True) [1/2*(7*x + 6)*e^(-x) + 1/2*sin(x), 'variationofparameters'] sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,[0,3,pi/2,2]) 3*((e^(1/2*pi) - 2)*x/pi + 1)*e^(-x) + 1/2*sin(x) sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,[0,3,pi/2,2],show_method=True) [3*((e^(1/2*pi) - 2)*x/pi + 1)*e^(-x) + 1/2*sin(x), 'variationofparameters'] sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y) (k2*x + k1)*e^(-x) sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,show_method=True) [(k2*x + k1)*e^(-x), 'constcoeff'] sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,[0,3,1]) (4*x + 3)*e^(-x) sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,[0,3,1],show_method=True) [(4*x + 3)*e^(-x), 'constcoeff'] sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,[0,3,pi/2,2]) (2*(2*e^(1/2*pi) - 3)*x/pi + 3)*e^(-x) sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,[0,3,pi/2,2],show_method=True) [(2*(2*e^(1/2*pi) - 3)*x/pi + 3)*e^(-x), 'constcoeff'] This ODE with separated variables is solved as exact Explanation: factor does not split exp(x-y) in maxima into exp(x)*exp(-y) sage: desolve(diff(y,x)==exp(x-y),y,show_method=True) [-e^x + e^y(x) == c, 'exact'] This equation can be solved within Maxima but not within Sage needs assumptions assume(x>0,y>0) and works in maxima, but not in Sage # sage: assume(x>0) # sage: assume(y>0) # sage: maxima("assume(x>0,y>0)") # sage: desolve(x*diff(y,x)-x*sqrt(y^2+x^2)-y,y,show_method=True) This equation can be solved within Maxima but not within Sage # sage: desolve(x*(x+1)*diff(y,x,2)+(x+5)*diff(y,x,1)+(-4)*y,y,contrib_ode=True) AUTHOR: David Joyner (1-2006) Robert Bradshaw (10-2008) CHANGES: Robert Marik (10-2009) Remark: The method contrib_ode has not been tested with initial conditions, since I was not able to find equation which cannot be solved using ode2, can be solved using contrib_ode and the result is sufficiently simple to continue safely (one SymbolicEquation) - robert.marik """ if is_SymbolicEquation(de): de = de.lhs() - de.rhs() if len(ivars) != 1: raise ValueError, "Unable to determine independent variable, please specify." ivar = ivars[0] de0 = de._maxima_() soln = de0.ode2(dvar, ivar) def sanitize_var(exprs): return exprs.replace("'"+dvar_str+"("+ivar_str+")",dvar_str) dvar_str=dvar.operator()._maxima_().str() ivar_str=ivar._maxima_().str() de00 = de._maxima_().str() de0 = sanitize_var(de00) ode_solver="ode2" 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) # we produce string like this # ode2('diff(y,x,2)+2*'diff(y,x,1)+y-cos(x),y(x),x) soln = maxima(cmd) if str(soln).strip() == 'false': raise NotImplementedError, "Maxima was unable to solve this system." def to_eqns(lhs, exprs): eqns = [] for lhs, expr in zip(lhs, exprs): if is_SymbolicEquation(expr): eqns.append(expr) else: if lhs == dvar and len(exprs) == 2: ivar_ic = exprs[0] # figure this out... lhs = lhs.subs(lhs.default_variable()==ivar_ic) eqns.append(lhs == expr) return eqns if ics is not None: if contrib_ode: ode_solver="contrib_ode" maxima("load('contrib_ode)") 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) # we produce string like this # (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)) soln = maxima(cmd) if str(soln).strip() == 'false': raise NotImplementedError, "Maxima was unable to solve this ODE." else: raise NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True." if show_method: maxima_method=maxima("method") if (ics is not None): if not is_SymbolicEquation(soln.sage()): raise NotImplementedError, "Maxima was unable to use initial condition for this equation (%s)"%(maxima_method.str()) if len(ics) == 2: ics = to_eqns([ivar, dvar], ics) soln = soln.ic1(ics[0], ics[1]) tempic=(ivar==ics[0])._maxima_().str() tempic=tempic+","+(dvar==ics[1])._maxima_().str() 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) cmd=sanitize_var(cmd) # we produce string like this # (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)) soln=maxima(cmd) if len(ics) == 3: ics = to_eqns([ivar, dvar, dvar.derivative(ivar)], ics) soln = soln.ic2(*ics) #fixed ic2 command from Maxima - we have to ensure that %k1, %k2 do not depend on variables, should be removed when fixed in Maxima 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))") tempic=(ivar==ics[0])._maxima_().str() tempic=tempic+","+(dvar==ics[1])._maxima_().str() tempic=tempic+",'diff("+dvar_str+","+ivar_str+")="+(ics[2])._maxima_().str() 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) cmd=sanitize_var(cmd) # we produce string like this # (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)) soln=maxima(cmd) if str(soln).strip() == 'false': raise NotImplementedError, "Maxima was unable to solve this IVP. Remove the initial contition to get the general solution." if len(ics) == 4: ics = to_eqns([ivar, dvar, ivar, dvar], ics) soln = soln.bc(*ics) if soln.lhs() == dvar: #fixed bc2 command from Maxima - we have to ensure that %k1, %k2 do not depend on variables, should be removed when fixed in Maxima 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))") 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()) cmd="(TEMP:%s,substitute(%s=%s(%s),TEMP))"%(cmd,dvar_str,dvar_str,ivar_str) cmd=sanitize_var(cmd) # we produce string like this # (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)) soln=maxima(cmd) if str(soln).strip() == 'false': raise NotImplementedError, "Maxima was unable to solve this BVP. Remove the initial contition to get the general solution." soln=soln.sage() if is_SymbolicEquation(soln) and soln.lhs() == dvar: # Remark: Here we do not check that the right hand side does not depend on dvar. # This probably will not hapen for soutions obtained via ode2, anyway. soln = soln.rhs() return soln.sage() if show_method: return [soln,maxima_method.str()] else: return soln #def desolve_laplace2(de,vars,ics=None): ##     """