# Ticket #10682: trac_10682-reviewer.patch

File trac_10682-reviewer.patch, 55.8 KB (added by Jean-Pierre Flori, 11 years ago)

Reviewer patch, formatting issues

• ## sage/calculus/calculus.py

# HG changeset patch
# User Jean-Pierre Flori <jean-pierre.flor@ssi.gouv.fr>
# Date 1331030079 -3600
# Node ID 4ebd0aae52ae9b89870a211578c815127fd98b5b
# Parent  0a84cd9c2c97b45d3bf8e2e8a10d337f7139b98f
#10682: Reviewer patch

diff --git a/sage/calculus/calculus.py b/sage/calculus/calculus.py
 a sage: CC(f) 1.12762596520638 + 1.17520119364380*I sage: ComplexField(200)(f) 1.1276259652063807852262251614026720125478471180986674836290 + 1.1752011936438014568823818505956008151557179813340958702296*I 1.1276259652063807852262251614026720125478471180986674836290 + 1.1752011936438014568823818505956008151557179813340958702296*I sage: ComplexField(100)(f) 1.1276259652063807852262251614 + 1.1752011936438014568823818506*I sage: f = sum(1/var('n%s'%i)^i for i in range(10)) sage: f 1/n1 + 1/n2^2 + 1/n3^3 + 1/n4^4 + 1/n5^5 + 1/n6^6 + 1/n7^7 + 1/n8^8 + 1/n9^9 + 1 1/n1 + 1/n2^2 + 1/n3^3 + 1/n4^4 + 1/n5^5 + 1/n6^6 + 1/n7^7 + 1/n8^8 + 1/n9^9 + 1 Note that after calling var, the variables are immediately available for use:: We can, of course, substitute:: sage: f(n9=9,n7=n6) 1/n1 + 1/n2^2 + 1/n3^3 + 1/n4^4 + 1/n5^5 + 1/n6^6 + 1/n6^7 + 1/n8^8 + 387420490/387420489 1/n1 + 1/n2^2 + 1/n3^3 + 1/n4^4 + 1/n5^5 + 1/n6^6 + 1/n6^7 + 1/n8^8 + 387420490/387420489 TESTS: - 5 - integral is probably divergent or slowly convergent - 6 - the input is invalid; this includes the case of desired_relative_error being too small to be achieved - 6 - the input is invalid; this includes the case of desired_relative_error being too small to be achieved ALIAS: nintegrate is the same as nintegral else: raise TypeError, err #This is just a work around until there is a response to #http://www.math.utexas.edu/pipermail/maxima/2008/012975.html # Maxima returns an unevaluated expression when the underlying SLATEC # library fails. See: # http://www.math.utexas.edu/pipermail/maxima/2008/012975.html # and Trac ticket #10682 for more information. if 'quad_qags' in str(v): raise ValueError, "Maxima (via quadpack) cannot compute the integral" sage: sin(pi/7).minpoly() x^6 - 7/4*x^4 + 7/8*x^2 - 7/64 sage: minpoly(exp(I*pi/17)) x^16 - x^15 + x^14 - x^13 + x^12 - x^11 + x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 x^16 - x^15 + x^14 - x^13 + x^12 - x^11 + x^10 - x^9 + x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1 Here we verify it gives the same result as the abstract number field. algorithm='numerical':: sage: cos(pi/33).minpoly(algorithm='algebraic') x^10 + 1/2*x^9 - 5/2*x^8 - 5/4*x^7 + 17/8*x^6 + 17/16*x^5 - 43/64*x^4 - 43/128*x^3 + 3/64*x^2 + 3/128*x + 1/1024 x^10 + 1/2*x^9 - 5/2*x^8 - 5/4*x^7 + 17/8*x^6 + 17/16*x^5 - 43/64*x^4 - 43/128*x^3 + 3/64*x^2 + 3/128*x + 1/1024 sage: cos(pi/33).minpoly(algorithm='numerical') Traceback (most recent call last): ... .. math:: F(s) = \frac{1}{2\pi i} \int_{\gamma-i\infty}^{\gamma + i\infty} e^{st} F(s) dt, F(s) = \frac{1}{2\pi i} \int_{\gamma-i\infty}^{\gamma + i\infty} e^{st} F(s) dt, where \gamma is chosen so that the contour path of integration is in the region of convergence of F(s). sage: inverse_laplace(L, s, t) t |--> t*cos(t) sage: inverse_laplace(1/(s^3+1), s, t) 1/3*(sqrt(3)*sin(1/2*sqrt(3)*t) - cos(1/2*sqrt(3)*t))*e^(1/2*t) + 1/3*e^(-t) 1/3*(sqrt(3)*sin(1/2*sqrt(3)*t) - cos(1/2*sqrt(3)*t))*e^(1/2*t) + 1/3*e^(-t) No explicit inverse Laplace transform, so one is returned formally as a function ilt:: sage: diff(u(x+h), x) D[0](u)(h + x) sage: taylor(u(x+h),h,0,4) 1/24*h^4*D[0, 0, 0, 0](u)(x) + 1/6*h^3*D[0, 0, 0](u)(x) + 1/2*h^2*D[0, 0](u)(x) + h*D[0](u)(x) + u(x) 1/24*h^4*D[0, 0, 0, 0](u)(x) + 1/6*h^3*D[0, 0, 0](u)(x) + 1/2*h^2*D[0, 0](u)(x) + h*D[0](u)(x) + u(x) We compute a Laplace transform::
• ## sage/calculus/desolvers.py

diff --git a/sage/calculus/desolvers.py b/sage/calculus/desolvers.py
 a - eulers_method_2x2_plot - Plots the sequence of points obtained from Euler's method. AUTHORS: AUTHORS: - David Joyner (3-2006) - Initial version of functions def desolve(de, dvar, ics=None, ivar=None, show_method=False, contrib_ode=False): r""" Solves a 1st or 2nd order linear ODE via maxima. Including IVP and BVP. *Use* desolve?  *if the output in truncated in notebook.* INPUT: - de - an expression or equation representing the ODE - dvar - the dependent variable (hereafter called y) - ics - (optional) the initial or boundary conditions - for a first-order equation, specify the initial x and y - for a second-order equation, specify the initial x, y, and dy/dx, i.e. write [x_0, y(x_0), y'(x_0)] - for a second-order boundary solution, specify initial and final x and y boundary conditions, i.e. write [x_0, y(x_0), x_1, y(x_1)]. - gives an error if the solution is not SymbolicEquation (as happens for final x and y boundary conditions, i.e. write [x_0, y(x_0), x_1, y(x_1)]. - gives an error if the solution is not SymbolicEquation (as happens for example for Clairaut 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 contain singular solution, for example) OUTPUT: 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 sage: y = function('y', x) sage: desolve(diff(y,x) + y - 1, y) (c + e^x)*e^(-x) :: sage: f = desolve(diff(y,x) + y - 1, y, ics=[10,2]); f :: sage: plot(f) We can also solve second-order differential equations.:: sage: x = var('x') sage: y = function('y', x) sage: de = diff(y,x,2) - y == x sage: desolve(de, y) k1*e^x + k2*e^(-x) - x :: :: sage: desolve(de, y, [0,1,pi/2,4]) sage: desolve(de, y, [0,1,pi/2,4]) 4*sin(x) + cos(x) :: 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'] For equations involving more variables we specify independent variable:: sage: a,b,c,n=var('a b c n') condition at x=0, since this point is singlar point of the equation. Anyway, if the solution should be bounded at x=0, then k2=0.:: sage: desolve(x^2*diff(y,x,x)+x*diff(y,x)+(x^2-4)*y==0,y) k1*bessel_j(2, x) + k2*bessel_y(2, x) Difficult ODE produces error:: sage: desolve(sqrt(y)*diff(y,x)+e^(y)+cos(x)-sin(x+y)==0,y) # not tested Traceback (click to the left for traceback) ... NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True." Difficult ODE produces error - moreover, takes a long time :: sage: desolve(sqrt(y)*diff(y,x)+e^(y)+cos(x)-sin(x+y)==0,y,contrib_ode=True) # not tested [[y(x) == c + log(x), y(x) == c*e^x], 'factor'] :: 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'] Traceback (click to the left for traceback) ... NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True." :: sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,2],show_method=True) # not tested 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'] TESTS: Trac #9961 fixed (allow assumptions on the dependent variable in desolve):: sage: y=function('y',x); assume(x>0); assume(y>0) sage: sage.calculus.calculus.maxima('domain:real')  # needed since Maxima 5.26.0 to get the answer as below real sage: desolve(x*diff(y,x)-x*sqrt(y^2+x^2)-y == 0, y, contrib_ode=True) [x - arcsinh(y(x)/x) == c] Trac #10682 updated Maxima to 5.26, and it started to show a different solution in the complex domain for the ODE above:: Trac #10682 updated Maxima to 5.26, and it started to show a different solution in the complex domain for the ODE above:: sage: sage.calculus.calculus.maxima('domain:complex')  # back to the dafault, complex, domain sage: sage.calculus.calculus.maxima('domain:complex')  # back to the default, complex, domain complex sage: desolve(x*diff(y,x)-x*sqrt(y^2+x^2)-y == 0, y, contrib_ode=True) [1/2*(2*x^2*sqrt(x^(-2)) - 2*x*sqrt(x^(-2))*arcsinh(y(x)/sqrt(x^2)) - 2*x*sqrt(x^(-2))*arcsinh(y(x)^2/(sqrt(y(x)^2)*x)) + log(4*(2*x^2*sqrt((x^2*y(x)^2 + y(x)^4)/x^2)*sqrt(x^(-2)) + x^2 + 2*y(x)^2)/x^2))/(x*sqrt(x^(-2))) == c] [1/2*(2*x^2*sqrt(x^(-2)) - 2*x*sqrt(x^(-2))*arcsinh(y(x)/sqrt(x^2)) - 2*x*sqrt(x^(-2))*arcsinh(y(x)^2/(sqrt(y(x)^2)*x)) + log(4*(2*x^2*sqrt((x^2*y(x)^2 + y(x)^4)/x^2)*sqrt(x^(-2)) + x^2 + 2*y(x)^2)/x^2))/(x*sqrt(x^(-2))) == c] Trac #6479 fixed:: x :: sage: desolve( diff(y,x,x) == 0, y, [0,1,1]) x + 1 sage: x=var('x'); f=function('f',x); k=var('k'); assume(k>0) sage: desolve(diff(f,x,2)/f==k,f,ivar=x) k1*e^(sqrt(k)*x) + k2*e^(-sqrt(k)*x) AUTHORS: - David Joyner (1-2006) - Robert Bradshaw (10-2008) raise ValueError, "Unable to determine independent variable, please specify." ivar = ivars[0] def sanitize_var(exprs): return exprs.replace("'"+dvar_str+"("+ivar_str+")",dvar_str) return exprs.replace("'"+dvar_str+"("+ivar_str+")",dvar_str) de00 = de._maxima_() P = de00.parent() dvar_str=P(dvar.operator()).str() de00 = de00.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) 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 = P(cmd) if contrib_ode: ode_solver="contrib_ode" P("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) 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 = P(cmd) if str(soln).strip() == 'false': raise NotImplementedError, "Maxima was unable to solve this ODE." 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." # (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=P(cmd) if len(ics) == 3: #fixed ic2 command from Maxima - we have to ensure that %k1, %k2 do not depend on variables, should be removed when fixed in Maxima #fixed ic2 command from Maxima - we have to ensure that %k1, %k2 do not depend on variables, should be removed when fixed in Maxima P("ic2_sage(soln,xa,ya,dya):=block([programmode:true,backsubst:true,singsolve:true,temp,%k2,%k1,TEMP_k], \ noteqn(xa), noteqn(ya), noteqn(dya), boundtest('%k1,%k1), boundtest('%k2,%k2), \ temp: lhs(soln) - rhs(soln), \ # (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=P(cmd) if str(soln).strip() == 'false': raise NotImplementedError, "Maxima was unable to solve this IVP. Remove the initial condition to get the general solution." raise NotImplementedError, "Maxima was unable to solve this IVP. Remove the initial condition to get the general solution." if len(ics) == 4: #fixed bc2 command from Maxima - we have to ensure that %k1, %k2 do not depend on variables, should be removed when fixed in Maxima #fixed bc2 command from Maxima - we have to ensure that %k1, %k2 do not depend on variables, should be removed when fixed in Maxima P("bc2_sage(soln,xa,ya,xb,yb):=block([programmode:true,backsubst:true,singsolve:true,temp,%k1,%k2,TEMP_k], \ 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]), \ 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,P(ivar==ics[0]).str(),dvar_str,P(ics[1]).str(),P(ivar==ics[2]).str(),dvar_str,P(ics[3]).str()) cmd="(TEMP:%s,substitute(%s=%s(%s),TEMP))"%(cmd,dvar_str,dvar_str,ivar_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=P(cmd) if str(soln).strip() == 'false': raise NotImplementedError, "Maxima was unable to solve this BVP. Remove the initial condition to get the general solution." raise NotImplementedError, "Maxima was unable to solve this BVP. Remove the initial condition to get the general solution." soln=soln.sage() if is_SymbolicEquation(soln) and soln.lhs() == dvar: soln = soln.rhs() if show_method: return [soln,maxima_method.str()] else: else: return soln #    #return maxima.eval(cmd) #    return de0.desolve(vars[0]).rhs() def desolve_laplace(de, dvar, ics=None, ivar=None): """ Solves an ODE using laplace transforms. Initials conditions are optional. INPUT: - de - a lambda expression representing the ODE (eg, de = diff(y,x,2) == diff(y,x)+sin(x)) Solution of the ODE as symbolic expression EXAMPLES:: sage: u=function('u',x) sage: eq = diff(u,x) - exp(-x) - u == 0 sage: desolve_laplace(eq,u) 1/2*(2*u(0) + 1)*e^x - 1/2*e^(-x) We can use initial conditions:: sage: desolve_laplace(eq,u,ics=[0,3]) -1/2*e^(-x) + 7/2*e^x The initial conditions do not persist in the system (as they persisted The initial conditions do not persist in the system (as they persisted in previous versions):: sage: desolve_laplace(eq,u) sage: eq = diff(f,x) + f == 0 sage: desolve_laplace(eq,f,[0,1]) e^(-x) :: sage: x = var('x') sage: f = function('f', x) sage: de = diff(f,x,x) - 2*diff(f,x) + f sage: soln(t=3) e^(-3) + 1 AUTHORS: AUTHORS: - David Joyner (1-2006,8-2007) - Robert Marik (10-2009) """ #This is the original code from David Joyner (inputs and outputs strings) ## verbatim copy from desolve - end def sanitize_var(exprs):  # 'y(x) -> y(x) return exprs.replace("'"+str(dvar),str(dvar)) return exprs.replace("'"+str(dvar),str(dvar)) de0=de._maxima_() P = de0.parent() cmd = sanitize_var("desolve("+de0.str()+","+str(dvar)+")") soln=P(cmd).rhs() if str(soln).strip() == 'false': raise NotImplementedError, "Maxima was unable to solve this ODE." raise NotImplementedError, "Maxima was unable to solve this ODE." soln=soln.sage() if ics!=None: d = len(ics) for i in range(0,d-1): soln=eval('soln.substitute(diff(dvar,ivar,i)('+str(ivar)+'=ics[0])==ics[i+1])') soln=eval('soln.substitute(diff(dvar,ivar,i)('+str(ivar)+'=ics[0])==ics[i+1])') return soln def desolve_system(des, vars, ics=None, ivar=None): """ Solves any size system of 1st order ODE's. Initials conditions are optional. Onedimensional systems are passed to :meth:desolve_laplace. Onedimensional systems are passed to :meth:desolve_laplace. INPUT: - des - list of ODEs - vars - list of dependent variables - ics - (optional) list of initial values for ivar and vars - ivar - (optional) the independent variable, which must be specified if there is more than one independent variable in the equation. sage: desolve_system([de1, de2], [x,y]) [x(t) == (x(0) - 1)*cos(t) - (y(0) - 1)*sin(t) + 1, y(t) == (x(0) - 1)*sin(t) + (y(0) - 1)*cos(t) + 1] Now we give some initial conditions:: sage: sol = desolve_system([de1, de2], [x,y], ics=[0,1,2]); sol [x(t) == -sin(t) + 1, y(t) == cos(t) + 1] :: sage: solnx, solny = sol[0].rhs(), sol[1].rhs() sage: plot([solnx,solny],(0,1))  # not tested sage: parametric_plot((solnx,solny),(0,1))  # not tested TESTS: Trac #9823 fixed:: sage: t = var('t') sage: x = function('x', t) sage: de1 = diff(x,t) + 1 == 0 sage: desolve_system([de1], [x]) -t + x(0) AUTHORS: AUTHORS: - Robert Bradshaw (10-2008) """ for dvar, ic in zip(dvars, ics[:1]): dvar.atvalue(ivar==ivar_ic, dvar) return soln def desolve_system_strings(des,vars,ics=None): r""" sage: s = var('s') sage: function('x', s) x(s) :: sage: function('y', s) y(s) :: sage: de1 = lambda z: diff(z[0],s) + z[1] - 1 sage: de2 = lambda z: diff(z[1],s) - z[0] + 1 sage: des = [de1([x(s),y(s)]),de2([x(s),y(s)])] sage: vars = ["s","x","y"] sage: desolve_system_strings(des,vars) ["(1-'y(0))*sin(s)+('x(0)-1)*cos(s)+1", "('x(0)-1)*sin(s)+('y(0)-1)*cos(s)+1"] ["(1-'y(0))*sin(s)+('x(0)-1)*cos(s)+1", "('x(0)-1)*sin(s)+('y(0)-1)*cos(s)+1"] :: sage: ics = [0,1,-1] sage: soln = desolve_system_strings(des,vars,ics); soln ['2*sin(s)+1', '1-2*cos(s)'] :: sage: solnx, solny = map(SR, soln) sage: RR(solnx(s=3)) 1.28224001611973 :: sage: P1 = plot([solnx,solny],(0,1)) sage: P2 = parametric_plot((solnx,solny),(0,1)) Now type show(P1), show(P2) to view these. AUTHORS: AUTHORS: - David Joyner (3-2006, 8-2007) """ last column. *For pedagogical purposes only.* EXAMPLES:: sage: from sage.calculus.desolvers import eulers_method sage: x,y = PolynomialRing(QQ,2,"xy").gens() sage: eulers_method(5*x+y-5,0,1,1/2,1) sage: P2 = line(pts) sage: (P1+P2).show() AUTHORS: AUTHORS: - David Joyner """ next (last) column. *For pedagogical purposes only.* EXAMPLES:: sage: from sage.calculus.desolvers import eulers_method_2x2 sage: t, x, y = PolynomialRing(QQ,3,"txy").gens() sage: f = x+y+t; g = x-y sage: eulers_method_2x2(f,g, 0, 0, 0, 1/3, 1,algorithm="none") [[0, 0, 0], [1/3, 0, 0], [2/3, 1/9, 0], [1, 10/27, 1/27], [4/3, 68/81, 4/27]] [[0, 0, 0], [1/3, 0, 0], [2/3, 1/9, 0], [1, 10/27, 1/27], [4/3, 68/81, 4/27]] :: 3/4                 0.63                   -0.0078               -0.031                 0.11 1                 0.63                     0.020                0.079                0.071 To numerically approximate y(1), where y''+ty'+y=0, y(0)=1, y'(0)=0:: To numerically approximate y(1), where y''+ty'+y=0, y(0)=1, y'(0)=0:: sage: t,x,y=PolynomialRing(RR,3,"txy").gens() sage: f = y; g = -x-y*t 3/4                 0.88                     -0.15                -0.62                -0.10 1                 0.75                     -0.17                -0.68               -0.015 AUTHORS: AUTHORS: - David Joyner """ x(a)=x_0, y' = g(t,x,y), y(a) = y_0. *For pedagogical purposes only.* EXAMPLES:: sage: from sage.calculus.desolvers import eulers_method_2x2_plot def desolve_rk4_determine_bounds(ics,end_points=None): """ Used to determine bounds for numerical integration. - If end_points is None, the interval for integration is from ics[0] to ics[0]+10 - If end_points is a or [a], the interval for integration is from min(ics[0],a) to max(ics[0],a) - If end_points is a or [a], the interval for integration is from min(ics[0],a) to max(ics[0],a) - If end_points is [a,b], the interval for integration is from min(ics[0],a) to max(ics[0],b) sage: from sage.calculus.desolvers import desolve_rk4_determine_bounds sage: desolve_rk4_determine_bounds([0,2],1) (0, 1) :: sage: desolve_rk4_determine_bounds([0,2]) sage: desolve_rk4_determine_bounds([0,2]) (0, 10) :: sage: desolve_rk4_determine_bounds([0,2],[-2]) (-2, 0) :: sage: desolve_rk4_determine_bounds([0,2],[-2,4]) (-2, 4) """ if end_points is None: if end_points is None: return((ics[0],ics[0]+10)) if not isinstance(end_points,list): end_points=[end_points] return (min(ics[0],end_points[0]),max(ics[0],end_points[0])) else: return (min(ics[0],end_points[0]),max(ics[0],end_points[1])) def desolve_rk4(de, dvar, ics=None, ivar=None, end_points=None, step=0.1, output='list', **kwds): """ equation. See also ode_solver. INPUT: input is similar to desolve command. The differential equation can be input is similar to desolve command. The differential equation can be written in a form close to the plot_slope_field or desolve command - Variant 1 (function in two variables) - de - right hand side, i.e. the function f(x,y) from ODE y'=f(x,y) - de - right hand side, i.e. the function f(x,y) from ODE y'=f(x,y) - dvar - dependent variable (symbolic variable declared by var) - Variant 2 (symbolic equation) - de - equation, including term with diff(y,x) - dvar - dependent variable (declared as funciton of independent variable) - dvar - dependent variable (declared as a function of independent variables) - Other parameters - Other parameters - ivar - should be specified, if there are more variables or if the equation is autonomous - ivar - should be specified, if there are more variables or if the equation is autonomous - ics - initial conditions in the form [x0,y0] - end_points - the end points of the interval - if end_points is a or [a], we integrate on between min(ics[0],a) and max(ics[0],a) - if end_points is a or [a], we integrate on between min(ics[0],a) and max(ics[0],a) - if end_points is None, we use end_points=ics[0]+10 - if end_points is [a,b] we integrate on between min(ics[0],a) and max(ics[0],b) - if end_points is [a,b] we integrate on between min(ics[0],a) and max(ics[0],b) - step - (optional, default:0.1) the length of the step (positive number) - step - (optional, default:0.1) the length of the step (positive number) - output - (optional, default: 'list') one of 'list', 'plot', 'slope_field' (graph of the solution with slope field) OUTPUT: OUTPUT: Returns a list of points, or plot produced by list_plot, optionally with slope field. Variant 1 for input - we can pass ODE in the form used by desolve function In this example we integrate bakwards, since end_points < ics[0]:: end_points < ics[0]:: sage: y=function('y',x) sage: desolve_rk4(diff(y,x)+y*(y-1) == x-2,y,ics=[1,1],step=0.5, end_points=0) sage: y=function('y',x) sage: desolve_rk4(diff(y,x)+y*(y-1) == x-2,y,ics=[1,1],step=0.5, end_points=0) [[0.0, 8.90425710896], [0.5, 1.90932794536], [1, 1]] Here we show how to plot simple pictures. For more advanced sage: x,y=var('x y') sage: P=desolve_rk4(y*(2-y),y,ics=[0,.1],ivar=x,output='slope_field',end_points=[-4,6],thickness=3) ALGORITHM: 4th order Runge-Kutta method. Wrapper for command rk in Maxima's dynamics package.  Perhaps could be faster by using fast_float instead. AUTHORS: AUTHORS: - Robert Marik (10-2009) """ if ics is None: if ics is None: raise ValueError, "No initial conditions, specify with ics=[x0,y0]." if ivar is None: YMIN=sol[0][1] XMAX=XMIN YMAX=YMIN for s,t in sol: for s,t in sol: if s>XMAX:XMAX=s if sYMAX:YMAX=t r""" Solves numerically system of first-order ordinary differential equations using the 4th order Runge-Kutta method. Wrapper for Maxima command rk. See also ode_solver. Maxima command rk. See also ode_solver. INPUT: input is similar to desolve_system and desolve_rk4 commands - des - right hand sides of the system - vars - dependent variables missing - ics - initial conditions in the form [x0,y01,y02,y03,....] - end_points - the end points of the interval - if end_points is a or [a], we integrate on between min(ics[0],a) and max(ics[0],a) - if end_points is a or [a], we integrate on between min(ics[0],a) and max(ics[0],a) - if end_points is None, we use end_points=ics[0]+10 - if end_points is [a,b] we integrate on between min(ics[0],a) and max(ics[0],b) - if end_points is [a,b] we integrate on between min(ics[0],a) and max(ics[0],b) - step -- (optional, default: 0.1) the length of the step sage: P=desolve_system_rk4([x*(1-y),-y*(1-x)],[x,y],ics=[0,0.5,2],ivar=t,end_points=20) sage: Q=[ [i,j] for i,j,k in P] sage: LP=list_plot(Q) sage: Q=[ [j,k] for i,j,k in P] sage: LP=list_plot(Q) ALGORITHM: 4th order Runge-Kutta method. Wrapper for command rk in Maxima's dynamics package.  Perhaps could be faster by using fast_float instead. AUTHOR: AUTHOR: - Robert Marik (10-2009) """ if ics is None: if ics is None: raise ValueError, "No initial conditions, specify with ics=[x0,y01,y02,...]." ivars = set([]) x0=ics[0] icss = [ics[i]._maxima_().str() for i in range(1,len(ics))] icstr = "[" + ",".join(icss) + "]" step=abs(step) step=abs(step) maxima("load('dynamics)") lower_bound,upper_bound=desolve_rk4_determine_bounds(ics,end_points) , rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0 , mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0): r""" Solves numerically a system of first-order ordinary differential equations Solves numerically a system of first-order ordinary differential equations using odeint from scipy.integrate module. INPUT: - des  -- right hand sides of the system - ics  -- initial conditions - ics  -- initial conditions - times -- a sequence of time points in which the solution must be found - ivar -- independent variable, optional. - compute_jac -- boolean. If True, the Jacobian of des is computed and - compute_jac -- boolean. If True, the Jacobian of des is computed and used during the integration of Stiff Systems. Default value is False. Other Parameters (taken from the documentation of odeint function from Other Parameters (taken from the documentation of odeint function from scipy.integrate module) - rtol, atol : float - hmin : float, (0: solver-determined) The minimum absolute step size allowed. - ixpr : boolean. - ixpr : boolean. Whether to generate extra printing at method switches. - mxstep : integer, (0: solver-determined) - mxords : integer, (0: solver-determined) Maximum order to be allowed for the stiff (BDF) method. OUTPUT: Returns a list with the solution of the system at each time in times. EXAMPLES: Lotka Volterra Equations:: sage: from sage.calculus.desolvers import desolve_odeint sage: sol=desolve_odeint(f,[0.5,2],srange(0,10,0.1),[x,y]) sage: p=line(zip(sol[:,0],sol[:,1])) sage: p.show() Lorenz Equations:: sage: x,y,z=var('x,y,z') sage: rho=28 sage: beta=8/3 sage: # The Lorenz equations sage: lorenz=[sigma*(y-x),x*(rho-z)-y,x*y-beta*z] sage: lorenz=[sigma*(y-x),x*(rho-z)-y,x*y-beta*z] sage: # Time and initial conditions sage: times=srange(0,50.05,0.05) sage: ics=[0,1,1] sage: sol=desolve_odeint(lorenz,ics,times,[x,y,z],rtol=1e-13,atol=1e-14) One-dimensional Stiff system:: sage: y= var('y') sage: sol=desolve_odeint(f,ic,t,y,rtol=1e-9,atol=1e-10,compute_jac=True) sage: p=points(zip(t,sol)) sage: p.show() Another Stiff system with some optional parameters with no Another Stiff system with some optional parameters with no default value:: sage: y1,y2,y3=var('y1,y2,y3') sage: t=srange(0,10,0.01) sage: v=[y1,y2,y3] sage: sol=desolve_odeint(f,ci,t,v,rtol=1e-3,atol=1e-4,h0=0.1,hmax=1,hmin=1e-4,mxstep=1000,mxords=17) AUTHOR: - Oriol Castejon (05-2010) """ def Dfun(y,t): v = list(y[:]) v.append(t) return [[element(*v) for element in row] for row in J] return [[element(*v) for element in row] for row in J] sol=odeint(func, ics, times, args=args, Dfun=Dfun, rtol=rtol, atol=atol, tcrit=tcrit, h0=h0, hmax=hmax, hmin=hmin, ixpr=ixpr, mxstep=mxstep, mxhnil=mxhnil, mxordn=mxordn, mxords=mxords, printmessg=printmessg)
• ## sage/symbolic/relation.py

diff --git a/sage/symbolic/relation.py b/sage/symbolic/relation.py
 a 0 == -10*a + b - 144 Subtract two symbolic equations:: sage: var('a,b') (a, b) (a, b) sage: m = 144 == 20 * a + b sage: n = 136 == 10 * a + b sage: m - n sage: assume(x>0, y < 2) sage: assumptions() [x > 0, y < 2] sage: (y < 2).forget() sage: (y < 2).forget() sage: assumptions() [x > 0] sage: forget() sage: eq = (x == 2) sage: eq._maple_init_() 'x = 2' Comparison:: sage: x = var('x') sage: f = x^5 + a sage: solve(f==0,x) [x == (-a)^(1/5)*e^(2/5*I*pi), x == (-a)^(1/5)*e^(4/5*I*pi), x == (-a)^(1/5)*e^(-4/5*I*pi), x == (-a)^(1/5)*e^(-2/5*I*pi), x == (-a)^(1/5)] [x == (-a)^(1/5)*e^(2/5*I*pi), x == (-a)^(1/5)*e^(4/5*I*pi), x == (-a)^(1/5)*e^(-4/5*I*pi), x == (-a)^(1/5)*e^(-2/5*I*pi), x == (-a)^(1/5)] You can also do arithmetic with inequalities, as illustrated below:: sage: eqn = x^3 + 2/3 >= x sage: loads(dumps(eqn)) x^3 + 2/3 >= x sage: loads(dumps(eqn)) == eqn sage: loads(dumps(eqn)) == eqn True AUTHORS: Used internally by the symbolic solve command to convert the output of Maxima's solve command to a list of solutions in Sage's symbolic package. EXAMPLES: We derive the (monic) quadratic formula:: sage: var('x,a,b') (x, a, b) sage: solve(x^2 + a*x + b == 0, x) [x == -1/2*a - 1/2*sqrt(a^2 - 4*b), x == -1/2*a + 1/2*sqrt(a^2 - 4*b)] Behind the scenes when the above is evaluated the function :func:string_to_list_of_solutions is called with input the string s below:: sage: s = '[x=-(sqrt(a^2-4*b)+a)/2,x=(sqrt(a^2-4*b)-a)/2]' sage: sage.symbolic.relation.string_to_list_of_solutions(s) [x == -1/2*a - 1/2*sqrt(a^2 - 4*b), x == -1/2*a + 1/2*sqrt(a^2 - 4*b)] def solve(f, *args, **kwds): r""" Algebraically solve an equation or system of equations (over the complex numbers) for given variables. Inequalities and systems of inequalities are also supported. complex numbers) for given variables. Inequalities and systems of inequalities are also supported. INPUT: -  f - equation or system of equations (given by a list or tuple) -  *args - variables to solve for. -  solution_dict - bool (default: False); if True or non-zero, -  solution_dict - bool (default: False); if True or non-zero, return a list of dictionaries containing the solutions. If there are no solutions, return an empty list (rather than a list containing an empty dictionary). Likewise, if there's only a single solution, return a list containing one dictionary with that solution. EXAMPLES:: sage: x, y = var('x, y') sage: solve([x+y==6, x-y==4], x, y) [[x == 5, y == 1]] sage: solve([x^2+y^2 == 1, y^2 == x^3 + x + 1], x, y) [[x == -1/2*I*sqrt(3) - 1/2, y == -sqrt(-1/2*I*sqrt(3) + 3/2)], [x == -1/2*I*sqrt(3) - 1/2, y == sqrt(-1/2*I*sqrt(3) + 3/2)], [x == 1/2*I*sqrt(3) - 1/2, y == -sqrt(1/2*I*sqrt(3) + 3/2)], [x == 1/2*I*sqrt(3) - 1/2, y == sqrt(1/2*I*sqrt(3) + 3/2)], [x == 0, y == -1], [[x == -1/2*I*sqrt(3) - 1/2, y == -sqrt(-1/2*I*sqrt(3) + 3/2)], [x == -1/2*I*sqrt(3) - 1/2, y == sqrt(-1/2*I*sqrt(3) + 3/2)], [x == 1/2*I*sqrt(3) - 1/2, y == -sqrt(1/2*I*sqrt(3) + 3/2)], [x == 1/2*I*sqrt(3) - 1/2, y == sqrt(1/2*I*sqrt(3) + 3/2)], [x == 0, y == -1], [x == 0, y == 1]] sage: solve([sqrt(x) + sqrt(y) == 5, x + y == 10], x, y) [[x == -5/2*I*sqrt(5) + 5, y == 5/2*I*sqrt(5) + 5], [x == 5/2*I*sqrt(5) + 5, y == -5/2*I*sqrt(5) + 5]] [[x == -5/2*I*sqrt(5) + 5, y == 5/2*I*sqrt(5) + 5], [x == 5/2*I*sqrt(5) + 5, y == -5/2*I*sqrt(5) + 5]] sage: solutions=solve([x^2+y^2 == 1, y^2 == x^3 + x + 1], x, y, solution_dict=True) sage: for solution in solutions: print solution[x].n(digits=3), ",", solution[y].n(digits=3) -0.500 - 0.866*I , -1.27 + 0.341*I 0.000 , -1.00 0.000 , 1.00 Whenever possible, answers will be symbolic, but with systems of equations, at times approximations will be given, due to the Whenever possible, answers will be symbolic, but with systems of equations, at times approximations will be given, due to the underlying algorithm in Maxima:: sage: sols = solve([x^3==y,y^2==x],[x,y]); sols[-1], sols[0] ([x == 0, y == 0], [x == (0.309016994375 + 0.951056516295*I),  y == (-0.809016994375 - 0.587785252292*I)]) ([x == 0, y == 0], [x == (0.309016994375 + 0.951056516295*I), y == (-0.809016994375 - 0.587785252292*I)]) sage: sols[0][0].rhs().pyobject().parent() Complex Double Field for symbolic expressions, which defaults to exact answers only:: sage: solve([y^6==y],y) [y == e^(2/5*I*pi), y == e^(4/5*I*pi), y == e^(-4/5*I*pi), y == e^(-2/5*I*pi), y == 1, y == 0] [y == e^(2/5*I*pi), y == e^(4/5*I*pi), y == e^(-4/5*I*pi), y == e^(-2/5*I*pi), y == 1, y == 0] sage: solve( [y^6 == y], y)==solve( y^6 == y, y) True TypeError: 5 is not a valid variable. If we ask for dictionaries containing the solutions, we get them:: sage: solve([x^2-1],x,solution_dict=True) [{x: -1}, {x: 1}] sage: solve([x^2-4*x+4],x,solution_dict=True) x: 2, y: 4 x: -2, y: 4 If there is a parameter in the answer, that will show up as a new variable.  In the following example, r1 is a real free variable (because of the r):: If there is a parameter in the answer, that will show up as a new variable. In the following example, r1 is a real free variable (because of the r):: sage: solve([x+y == 3, 2*x+2*y == 6],x,y) [[x == -r1 + 3, y == r1]] solutions are returned. E.g., note that the first 3==3 evaluates to True, not to a symbolic equation. :: sage: solve([3==3, 1.00000000000000*x^3 == 0], x) [x == 0] sage: solve([1.00000000000000*x^3 == 0], x) [x == 0] Here, the first equation evaluates to False, so there are no solutions:: sage: solve([1==3, 1.00000000000000*x^3 == 0], x) [] Completely symbolic solutions are supported:: sage: var('s,j,b,m,g') (s, j, b, m, g) sage: sys = [ m*(1-s) - b*s*j, b*s*j-g*j ]; Use use_grobner if no solution is obtained from to_poly_solve:: sage: x,y=var('x y'); c1(x,y)=(x-5)^2+y^2-16; c2(x,y)=(y-3)^2+x^2-9 sage: solve([c1(x,y),c2(x,y)],[x,y]) [[x == -9/68*sqrt(55) + 135/68, y == -15/68*sqrt(5)*sqrt(11) + 123/68], [x == 9/68*sqrt(55) + 135/68, y == 15/68*sqrt(5)*sqrt(11) + 123/68]] sage: solve([c1(x,y),c2(x,y)],[x,y]) [[x == -9/68*sqrt(55) + 135/68, y == -15/68*sqrt(5)*sqrt(11) + 123/68], [x == 9/68*sqrt(55) + 135/68, y == 15/68*sqrt(5)*sqrt(11) + 123/68]] TESTS:: sage: solve([sin(x)==x,y^2==x],x,y) variables = args else: variables = tuple(args[0]) for v in variables: if not is_SymbolicVariable(v): raise TypeError, "%s is not a valid variable."%v Return all solutions to an equation or list of equations modulo the given integer modulus. Each equation must involve only polynomials in 1 or many variables. By default the solutions are returned as n-tuples, where n is the number of variables appearing anywhere in the given equations. The variables are in alphabetical order. INPUT: -  eqns - equation or list of equations -  modulus - an integer -  solution_dict - bool (default: False); if True or non-zero, return a list of dictionaries containing the solutions. If there are no solutions, return an empty list (rather than a list containing an empty dictionary). Likewise, if there's only a single solution, return a list containing one dictionary with that solution. EXAMPLES:: sage: var('x,y') (x, y) sage: solve_mod([x^2 + 2 == x, x^2 + y == y^2], 14) [(4, 2), (4, 6), (4, 9), (4, 13)] sage: solve_mod([x^2 == 1, 4*x  == 11], 15) [(14,)] Fermat's equation modulo 3 with exponent 5:: sage: var('x,y,z') (x, y, z) sage: solve_mod([x^5 + y^5 == z^5], 3) [(0, 0, 0), (0, 1, 1), (0, 2, 2), (1, 0, 1), (1, 1, 2), (1, 2, 0), (2, 0, 2), (2, 1, 0), (2, 2, 1)] [(0, 0, 0), (0, 1, 1), (0, 2, 2), (1, 0, 1), (1, 1, 2), (1, 2, 0), (2, 0, 2), (2, 1, 0), (2, 2, 1)] We can solve with respect to a bigger modulus if it consists only of small prime factors:: We can solve with respect to a bigger modulus if it consists only of small prime factors:: sage: [d] = solve_mod([5*x + y == 3, 2*x - 3*y == 9], 3*5*7*11*19*23*29, solution_dict = True) sage: d[x] is large:: sage: sorted(solve_mod([x^2 == 41], 10^20)) [(4538602480526452429,), (11445932736758703821,), (38554067263241296179,), (45461397519473547571,), (54538602480526452429,), (61445932736758703821,), (88554067263241296179,), (95461397519473547571,)] [(4538602480526452429,), (11445932736758703821,), (38554067263241296179,), (45461397519473547571,), (54538602480526452429,), (61445932736758703821,), (88554067263241296179,), (95461397519473547571,)] We solve a simple equation modulo 2:: sage: x,y = var('x,y') sage: solve_mod([x == y], 2) [(0, 0), (1, 1)] [(8, 5, 6)] Confirm that modulus 1 now behaves as it should:: sage: x, y = var("x y") sage: solve_mod([x==1], 1) [(0,)] sage: solve_mod([2*x^2+x*y, -x*y+2*y^2+x-2*y, -2*x^2+2*x*y-y^2-x-y], 1) [(0, 0)] """ from sage.rings.all import Integer, Integers, crt_basis raise ValueError, "the modulus must be a positive integer" vars = list(set(sum([list(e.variables()) for e in eqns], []))) vars.sort(key=repr) if modulus == 1: # degenerate case ans = [tuple(Integers(1)(0) for v in vars)] return ans factors = modulus.factor() crt_basis = vector(Integers(modulus), crt_basis([p**i for p,i in factors])) solutions = [] has_solution = True for p,i in factors: solution =_solve_mod_prime_power(eqns, p, i, vars) else: has_solution = False break ans = [] if has_solution: Internal help function for solve_mod, does little checking since it expects solve_mod to do that Return all solutions to an equation or list of equations modulo p^m. Return all solutions to an equation or list of equations modulo p^m. Each equation must involve only polynomials in 1 or many variables. The solutions are returned as n-tuples, where n is the number of variables in vars. INPUT: -  eqns - equation or list of equations -  p - a prime -  i - an integer > 0 -  vars - a list of variables to solve for EXAMPLES:: sage: var('x,y') (x, y) sage: solve_mod([x^2 + 2 == x, x^2 + y == y^2], 14) [(4, 2), (4, 6), (4, 9), (4, 13)] sage: solve_mod([x^2 == 1, 4*x  == 11], 15) [(14,)] Fermat's equation modulo 3 with exponent 5:: sage: var('x,y,z') (x, y, z) sage: solve_mod([x^5 + y^5 == z^5], 3) [(0, 0, 0), (0, 1, 1), (0, 2, 2), (1, 0, 1), (1, 1, 2), (1, 2, 0), (2, 0, 2), (2, 1, 0), (2, 2, 1)] [(0, 0, 0), (0, 1, 1), (0, 2, 2), (1, 0, 1), (1, 1, 2), (1, 2, 0), (2, 0, 2), (2, 1, 0), (2, 2, 1)] We solve a simple equation modulo 2:: sage: x,y = var('x,y') sage: solve_mod([x == y], 2) [(0, 0), (1, 1)] .. warning:: Currently this constructs possible solutions by building up from the smallest prime factor of the modulus.  The interface is good, but the algorithm is horrible if the modulus isn't the interface. TESTS: Confirm we can reproduce the first few terms of OEIS A187719:: sage: from sage.symbolic.relation import _solve_mod_prime_power sage: [sorted(_solve_mod_prime_power([x^2==41], 10, i, [x]))[0][0] for i in [1..13]] [1, 21, 71, 1179, 2429, 47571, 1296179, 8703821, 26452429, 526452429, 13241296179, 19473547571, 2263241296179] [1, 21, 71, 1179, 2429, 47571, 1296179, 8703821, 26452429, 526452429, 13241296179, 19473547571, 2263241296179] """ from sage.rings.all import Integers, PolynomialRing from sage.modules.all import vector from sage.misc.all import cartesian_product_iterator mrunning = 1 ans = [] for mi in xrange(m): def solve_ineq_univar(ineq): """ Function solves rational inequality in one variable. Function solves rational inequality in one variable. INPUT: - ineq - inequality in one variable OUTPUT: - list -- output is list of solutions as a list of simple inequalities output [A,B,C] means (A or B or C) each A, B, C is again a list and if A=[a,b], then A means (a and b). The list is empty if there is no output [A,B,C] means (A or B or C) each A, B, C is again a list and if A=[a,b], then A means (a and b). The list is empty if there is no solution. EXAMPLES:: sage: from sage.symbolic.relation import solve_ineq_univar sage: solve_ineq_univar(x-1/x>0) [[x > -1, x < 0], [x > 1]] sage: solve_ineq_univar(x^2-1/x>0) [[x < 0], [x > 1]] sage: solve_ineq_univar((x^3-1)*x<=0) [[x >= 0, x <= 1]] ALGORITHM: Calls Maxima command solve_rat_ineq AUTHORS: - Robert Marik (01-2010) """ ineqvar = ineq.variables() ineq0 = ineq._maxima_() ineq0.parent().eval("if solve_rat_ineq_loaded#true then (solve_rat_ineq_loaded:true,load(\"solve_rat_ineq.mac\")) ") sol = ineq0.solve_rat_ineq().sage() if repr(sol)=="all": from sage.rings.infinity import Infinity sol = [ineqvar[0]4],[x,y]) [[y + 4 < x, x < -y + 9, y < (5/2)]] [[y + 4 < x, x < -y + 9, y < (5/2)]] sage: solve_ineq_fourier([x+y<9,x-y>4],[y,x]) [[y < min(x - 4, -x + 9)]] sage: solve_ineq_fourier([x^2>=0]) [[x < +Infinity]] [[y < x, 0 < y]] ALGORITHM: Calls Maxima command fourier_elim AUTHORS: - Robert Marik (01-2010) """ if vars is None: if not atom(x) and op(x)=\"or\" then args(x) \ else [x]") sol = sol.or_to_list().sage() if repr(sol) == "[emptyset]": if repr(sol) == "[emptyset]": sol = [] if repr(sol) == "[universalset]": from sage.rings.infinity import Infinity from sage.rings.infinity import Infinity sol = [[i3) sage: solve_ineq(x^2-1>3) [[x < -2], [x > 2]] sage: solve_ineq(1/(x-1)<=8) [[x < 1], [x >= (9/8)]] System of inequalities with automatically detected inequalities:: sage: y=var('y') [[x < y, y < -x + 3, x < (3/2)]] ALGORITHM: Calls solve_ineq_fourier if inequalities are list and solve_ineq_univar of the inequality is symbolic expression. See the description of these commands for more details related to the there is no solution. AUTHORS: - Robert Marik (01-2010) """ if isinstance(ineq,list):