Ticket #10682: trac_10682reviewer.patch
File trac_10682reviewer.patch, 55.8 KB (added by , 10 years ago) 


sage/calculus/calculus.py
# HG changeset patch # User JeanPierre Flori <jeanpierre.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 b 247 247 sage: CC(f) 248 248 1.12762596520638 + 1.17520119364380*I 249 249 sage: ComplexField(200)(f) 250 1.1276259652063807852262251614026720125478471180986674836290 + 1.1752011936438014568823818505956008151557179813340958702296*I 250 1.1276259652063807852262251614026720125478471180986674836290 251 + 1.1752011936438014568823818505956008151557179813340958702296*I 251 252 sage: ComplexField(100)(f) 252 253 1.1276259652063807852262251614 + 1.1752011936438014568823818506*I 253 254 … … 256 257 257 258 sage: f = sum(1/var('n%s'%i)^i for i in range(10)) 258 259 sage: f 259 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 260 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 261 + 1 260 262 261 263 Note that after calling var, the variables are immediately 262 264 available for use:: … … 267 269 We can, of course, substitute:: 268 270 269 271 sage: f(n9=9,n7=n6) 270 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 272 1/n1 + 1/n2^2 + 1/n3^3 + 1/n4^4 + 1/n5^5 + 1/n6^6 + 1/n6^7 + 1/n8^8 273 + 387420490/387420489 271 274 272 275 TESTS: 273 276 … … 632 635  ``5``  integral is probably divergent or slowly 633 636 convergent 634 637 635  ``6``  the input is invalid; this includes the case of desired_relative_error636 being too small to be achieved638  ``6``  the input is invalid; this includes the case of 639 desired_relative_error being too small to be achieved 637 640 638 641 ALIAS: nintegrate is the same as nintegral 639 642 … … 716 719 else: 717 720 raise TypeError, err 718 721 719 #This is just a work around until there is a response to 720 #http://www.math.utexas.edu/pipermail/maxima/2008/012975.html 722 # Maxima returns an unevaluated expression when the underlying SLATEC 723 # library fails. See: 724 # http://www.math.utexas.edu/pipermail/maxima/2008/012975.html 725 # and Trac ticket #10682 for more information. 721 726 if 'quad_qags' in str(v): 722 727 raise ValueError, "Maxima (via quadpack) cannot compute the integral" 723 728 … … 806 811 sage: sin(pi/7).minpoly() 807 812 x^6  7/4*x^4 + 7/8*x^2  7/64 808 813 sage: minpoly(exp(I*pi/17)) 809 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 814 x^16  x^15 + x^14  x^13 + x^12  x^11 + x^10  x^9 + x^8  x^7 + x^6 815  x^5 + x^4  x^3 + x^2  x + 1 810 816 811 817 Here we verify it gives the same result as the abstract number 812 818 field. … … 884 890 algorithm='numerical':: 885 891 886 892 sage: cos(pi/33).minpoly(algorithm='algebraic') 887 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 893 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 894  43/128*x^3 + 3/64*x^2 + 3/128*x + 1/1024 888 895 sage: cos(pi/33).minpoly(algorithm='numerical') 889 896 Traceback (most recent call last): 890 897 ... … … 1306 1313 1307 1314 .. math:: 1308 1315 1309 1316 F(s) = \frac{1}{2\pi i} \int_{\gammai\infty}^{\gamma + i\infty} e^{st} F(s) dt, 1310 1317 1311 1318 where `\gamma` is chosen so that the contour path of 1312 1319 integration is in the region of convergence of `F(s)`. … … 1327 1334 sage: inverse_laplace(L, s, t) 1328 1335 t > t*cos(t) 1329 1336 sage: inverse_laplace(1/(s^3+1), s, t) 1330 1/3*(sqrt(3)*sin(1/2*sqrt(3)*t)  cos(1/2*sqrt(3)*t))*e^(1/2*t) + 1/3*e^(t) 1337 1/3*(sqrt(3)*sin(1/2*sqrt(3)*t)  cos(1/2*sqrt(3)*t))*e^(1/2*t) 1338 + 1/3*e^(t) 1331 1339 1332 1340 No explicit inverse Laplace transform, so one is returned formally 1333 1341 as a function ``ilt``:: … … 1369 1377 sage: diff(u(x+h), x) 1370 1378 D[0](u)(h + x) 1371 1379 sage: taylor(u(x+h),h,0,4) 1372 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) 1380 1/24*h^4*D[0, 0, 0, 0](u)(x) + 1/6*h^3*D[0, 0, 0](u)(x) 1381 + 1/2*h^2*D[0, 0](u)(x) + h*D[0](u)(x) + u(x) 1373 1382 1374 1383 We compute a Laplace transform:: 1375 1384 
sage/calculus/desolvers.py
diff git a/sage/calculus/desolvers.py b/sage/calculus/desolvers.py
a b 35 35  ``eulers_method_2x2_plot``  Plots the sequence of points obtained 36 36 from Euler's method. 37 37 38 AUTHORS: 38 AUTHORS: 39 39 40 40  David Joyner (32006)  Initial version of functions 41 41 … … 68 68 def desolve(de, dvar, ics=None, ivar=None, show_method=False, contrib_ode=False): 69 69 r""" 70 70 Solves a 1st or 2nd order linear ODE via maxima. Including IVP and BVP. 71 71 72 72 *Use* ``desolve? <tab>`` *if the output in truncated in notebook.* 73 73 74 74 INPUT: 75 75 76 76  ``de``  an expression or equation representing the ODE 77 77 78 78  ``dvar``  the dependent variable (hereafter called ``y``) 79 79 80 80  ``ics``  (optional) the initial or boundary conditions 81 81 82 82  for a firstorder equation, specify the initial ``x`` and ``y`` 83 83 84 84  for a secondorder equation, specify the initial ``x``, ``y``, 85 85 and ``dy/dx``, i.e. write `[x_0, y(x_0), y'(x_0)]` 86 86 87 87  for a secondorder boundary solution, specify initial and 88 final ``x`` and ``y`` boundary conditions, i.e. write `[x_0, y(x_0), x_1, y(x_1)]`. 89 90  gives an error if the solution is not SymbolicEquation (as happens for 88 final ``x`` and ``y`` boundary conditions, i.e. write 89 `[x_0, y(x_0), x_1, y(x_1)]`. 90 91  gives an error if the solution is not SymbolicEquation (as happens for 91 92 example for Clairaut equation) 92 93 93 94  ``ivar``  (optional) the independent variable (hereafter called 94 95 x), which must be specified if there is more than one 95 96 independent variable in the equation. 96 97 97 98  ``show_method``  (optional) if true, then Sage returns pair 98 99 ``[solution, method]``, where method is the string describing 99 100 method which has been used to get solution (Maxima uses the … … 111 112 contain singular solution, for example) 112 113 113 114 OUTPUT: 114 115 115 116 In most cases returns SymbolicEquation which defines the solution 116 117 implicitly. If the result is in the form y(x)=... (happens for 117 118 linear eqs.), returns the righthand side only. The possible … … 124 125 sage: y = function('y', x) 125 126 sage: desolve(diff(y,x) + y  1, y) 126 127 (c + e^x)*e^(x) 127 128 128 129 :: 129 130 130 131 sage: f = desolve(diff(y,x) + y  1, y, ics=[10,2]); f … … 133 134 :: 134 135 135 136 sage: plot(f) 136 137 137 138 We can also solve secondorder differential equations.:: 138 139 139 140 sage: x = var('x') 140 141 sage: y = function('y', x) 141 142 sage: de = diff(y,x,2)  y == x 142 143 sage: desolve(de, y) 143 144 k1*e^x + k2*e^(x)  x 144 145 145 146 146 147 :: 147 148 … … 166 167 167 168 :: 168 169 169 sage: desolve(de, y, [0,1,pi/2,4]) 170 sage: desolve(de, y, [0,1,pi/2,4]) 170 171 4*sin(x) + cos(x) 171 172 172 173 :: … … 178 179 179 180 sage: desolve(diff(y,x)^2+x*diff(y,x)y==0,y,contrib_ode=True,show_method=True) 180 181 [[y(x) == c^2 + c*x, y(x) == 1/4*x^2], 'clairault'] 181 182 182 183 For equations involving more variables we specify independent variable:: 183 184 184 185 sage: a,b,c,n=var('a b c n') … … 248 249 condition at x=0, since this point is singlar point of the 249 250 equation. Anyway, if the solution should be bounded at x=0, then 250 251 k2=0.:: 251 252 252 253 sage: desolve(x^2*diff(y,x,x)+x*diff(y,x)+(x^24)*y==0,y) 253 254 k1*bessel_j(2, x) + k2*bessel_y(2, x) 254 255 255 256 Difficult ODE produces error:: 256 257 257 258 sage: desolve(sqrt(y)*diff(y,x)+e^(y)+cos(x)sin(x+y)==0,y) # not tested 258 259 Traceback (click to the left for traceback) 259 260 ... 260 261 NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True." 261 262 262 263 Difficult ODE produces error  moreover, takes a long time :: 263 264 264 265 sage: desolve(sqrt(y)*diff(y,x)+e^(y)+cos(x)sin(x+y)==0,y,contrib_ode=True) # not tested … … 269 270 [[y(x) == c + log(x), y(x) == c*e^x], 'factor'] 270 271 271 272 :: 272 273 273 274 sage: desolve(diff(y,x)==(x+y)^2,y,contrib_ode=True,show_method=True) 274 275 [[[x == c  arctan(sqrt(t)), y(x) == x  sqrt(t)], [x == c + arctan(sqrt(t)), y(x) == x + sqrt(t)]], 'lagrange'] 275 276 … … 281 282 Traceback (click to the left for traceback) 282 283 ... 283 284 NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True." 284 285 285 286 :: 286 287 287 288 sage: desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,2],show_method=True) # not tested … … 293 294 294 295 sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y) 295 296 (k2*x + k1)*e^(x) + 1/2*sin(x) 296 297 297 298 :: 298 299 299 300 sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,show_method=True) 300 301 [(k2*x + k1)*e^(x) + 1/2*sin(x), 'variationofparameters'] 301 302 302 303 :: 303 304 304 305 sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,[0,3,1]) 305 306 1/2*(7*x + 6)*e^(x) + 1/2*sin(x) 306 307 307 308 :: 308 309 309 310 sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,[0,3,1],show_method=True) 310 311 [1/2*(7*x + 6)*e^(x) + 1/2*sin(x), 'variationofparameters'] 311 312 312 313 :: 313 314 314 315 sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,[0,3,pi/2,2]) 315 316 3*((e^(1/2*pi)  2)*x/pi + 1)*e^(x) + 1/2*sin(x) 316 317 317 318 :: 318 319 319 320 sage: desolve(diff(y,x,2)+2*diff(y,x)+y == cos(x),y,[0,3,pi/2,2],show_method=True) 320 321 [3*((e^(1/2*pi)  2)*x/pi + 1)*e^(x) + 1/2*sin(x), 'variationofparameters'] 321 322 322 323 :: 323 324 324 325 sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y) 325 326 (k2*x + k1)*e^(x) 326 327 327 328 :: 328 329 329 330 sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,show_method=True) 330 331 [(k2*x + k1)*e^(x), 'constcoeff'] 331 332 332 333 :: 333 334 334 335 sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,[0,3,1]) 335 336 (4*x + 3)*e^(x) 336 337 337 338 :: 338 339 339 340 sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,[0,3,1],show_method=True) 340 341 [(4*x + 3)*e^(x), 'constcoeff'] 341 342 342 343 :: 343 344 344 345 sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,[0,3,pi/2,2]) 345 346 (2*(2*e^(1/2*pi)  3)*x/pi + 3)*e^(x) 346 347 347 348 :: 348 349 349 350 sage: desolve(diff(y,x,2)+2*diff(y,x)+y == 0,y,[0,3,pi/2,2],show_method=True) 350 351 [(2*(2*e^(1/2*pi)  3)*x/pi + 3)*e^(x), 'constcoeff'] 351 352 352 353 TESTS: 353 354 354 355 Trac #9961 fixed (allow assumptions on the dependent variable in desolve):: 355 356 356 357 sage: y=function('y',x); assume(x>0); assume(y>0) 357 358 sage: sage.calculus.calculus.maxima('domain:real') # needed since Maxima 5.26.0 to get the answer as below 358 359 real 359 360 sage: desolve(x*diff(y,x)x*sqrt(y^2+x^2)y == 0, y, contrib_ode=True) 360 361 [x  arcsinh(y(x)/x) == c] 361 362 362 Trac #10682 updated Maxima to 5.26, and it started to show a different solution in the complex domain for the ODE above:: 363 Trac #10682 updated Maxima to 5.26, and it started to show a different 364 solution in the complex domain for the ODE above:: 363 365 364 sage: sage.calculus.calculus.maxima('domain:complex') # back to the d afault, complex, domain366 sage: sage.calculus.calculus.maxima('domain:complex') # back to the default, complex, domain 365 367 complex 366 368 sage: desolve(x*diff(y,x)x*sqrt(y^2+x^2)y == 0, y, contrib_ode=True) 367 [1/2*(2*x^2*sqrt(x^(2))  2*x*sqrt(x^(2))*arcsinh(y(x)/sqrt(x^2)) 368  2*x*sqrt(x^(2))*arcsinh(y(x)^2/(sqrt(y(x)^2)*x)) 369 + 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] 369 [1/2*(2*x^2*sqrt(x^(2))  2*x*sqrt(x^(2))*arcsinh(y(x)/sqrt(x^2)) 370  2*x*sqrt(x^(2))*arcsinh(y(x)^2/(sqrt(y(x)^2)*x)) 371 + log(4*(2*x^2*sqrt((x^2*y(x)^2 + y(x)^4)/x^2)*sqrt(x^(2)) 372 + x^2 + 2*y(x)^2)/x^2))/(x*sqrt(x^(2))) == c] 370 373 371 374 Trac #6479 fixed:: 372 375 … … 376 379 x 377 380 378 381 :: 379 382 380 383 sage: desolve( diff(y,x,x) == 0, y, [0,1,1]) 381 384 x + 1 382 385 … … 394 397 sage: x=var('x'); f=function('f',x); k=var('k'); assume(k>0) 395 398 sage: desolve(diff(f,x,2)/f==k,f,ivar=x) 396 399 k1*e^(sqrt(k)*x) + k2*e^(sqrt(k)*x) 397 398 400 401 399 402 AUTHORS: 400 403 401 404  David Joyner (12006) 402 405 403 406  Robert Bradshaw (102008) … … 419 422 raise ValueError, "Unable to determine independent variable, please specify." 420 423 ivar = ivars[0] 421 424 def sanitize_var(exprs): 422 return exprs.replace("'"+dvar_str+"("+ivar_str+")",dvar_str) 425 return exprs.replace("'"+dvar_str+"("+ivar_str+")",dvar_str) 423 426 de00 = de._maxima_() 424 427 P = de00.parent() 425 428 dvar_str=P(dvar.operator()).str() … … 427 430 de00 = de00.str() 428 431 de0 = sanitize_var(de00) 429 432 ode_solver="ode2" 430 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) 433 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) 431 434 # we produce string like this 432 435 # ode2('diff(y,x,2)+2*'diff(y,x,1)+ycos(x),y(x),x) 433 436 soln = P(cmd) … … 436 439 if contrib_ode: 437 440 ode_solver="contrib_ode" 438 441 P("load('contrib_ode)") 439 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) 442 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) 440 443 # we produce string like this 441 444 # (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)) 442 445 soln = P(cmd) 443 446 if str(soln).strip() == 'false': 444 raise NotImplementedError, "Maxima was unable to solve this ODE." 447 raise NotImplementedError, "Maxima was unable to solve this ODE." 445 448 else: 446 449 raise NotImplementedError, "Maxima was unable to solve this ODE. Consider to set option contrib_ode to True." 447 450 … … 462 465 # (TEMP:ic2(ode2('diff(y,x,2)+2*'diff(y,x,1)+ycos(x),y,x),x=0,y=3,'diff(y,x)=1),substitute(y=y(x),TEMP)) 463 466 soln=P(cmd) 464 467 if len(ics) == 3: 465 #fixed ic2 command from Maxima  we have to ensure that %k1, %k2 do not depend on variables, should be removed when fixed in Maxima 468 #fixed ic2 command from Maxima  we have to ensure that %k1, %k2 do not depend on variables, should be removed when fixed in Maxima 466 469 P("ic2_sage(soln,xa,ya,dya):=block([programmode:true,backsubst:true,singsolve:true,temp,%k2,%k1,TEMP_k], \ 467 470 noteqn(xa), noteqn(ya), noteqn(dya), boundtest('%k1,%k1), boundtest('%k2,%k2), \ 468 471 temp: lhs(soln)  rhs(soln), \ … … 479 482 # (TEMP:ic2(ode2('diff(y,x,2)+2*'diff(y,x,1)+ycos(x),y,x),x=0,y=3,'diff(y,x)=1),substitute(y=y(x),TEMP)) 480 483 soln=P(cmd) 481 484 if str(soln).strip() == 'false': 482 raise NotImplementedError, "Maxima was unable to solve this IVP. Remove the initial condition to get the general solution." 485 raise NotImplementedError, "Maxima was unable to solve this IVP. Remove the initial condition to get the general solution." 483 486 if len(ics) == 4: 484 #fixed bc2 command from Maxima  we have to ensure that %k1, %k2 do not depend on variables, should be removed when fixed in Maxima 487 #fixed bc2 command from Maxima  we have to ensure that %k1, %k2 do not depend on variables, should be removed when fixed in Maxima 485 488 P("bc2_sage(soln,xa,ya,xb,yb):=block([programmode:true,backsubst:true,singsolve:true,temp,%k1,%k2,TEMP_k], \ 486 489 noteqn(xa), noteqn(ya), noteqn(xb), noteqn(yb), boundtest('%k1,%k1), boundtest('%k2,%k2), \ 487 490 TEMP_k:solve([subst([xa,ya],soln), subst([xb,yb],soln)], [%k1,%k2]), \ … … 489 492 temp: maplist(lambda([zz], subst(zz,soln)),TEMP_k), \ 490 493 if length(temp)=1 then return(first(temp)) else return(temp))") 491 494 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()) 492 cmd="(TEMP:%s,substitute(%s=%s(%s),TEMP))"%(cmd,dvar_str,dvar_str,ivar_str) 495 cmd="(TEMP:%s,substitute(%s=%s(%s),TEMP))"%(cmd,dvar_str,dvar_str,ivar_str) 493 496 cmd=sanitize_var(cmd) 494 497 # we produce string like this 495 498 # (TEMP:bc2(ode2('diff(y,x,2)+2*'diff(y,x,1)+ycos(x),y,x),x=0,y=3,x=%pi/2,y=2),substitute(y=y(x),TEMP)) 496 499 soln=P(cmd) 497 500 if str(soln).strip() == 'false': 498 raise NotImplementedError, "Maxima was unable to solve this BVP. Remove the initial condition to get the general solution." 501 raise NotImplementedError, "Maxima was unable to solve this BVP. Remove the initial condition to get the general solution." 499 502 500 503 soln=soln.sage() 501 504 if is_SymbolicEquation(soln) and soln.lhs() == dvar: … … 504 507 soln = soln.rhs() 505 508 if show_method: 506 509 return [soln,maxima_method.str()] 507 else: 510 else: 508 511 return soln 509 512 510 513 … … 549 552 # #return maxima.eval(cmd) 550 553 # return de0.desolve(vars[0]).rhs() 551 554 552 555 553 556 def desolve_laplace(de, dvar, ics=None, ivar=None): 554 557 """ 555 558 Solves an ODE using laplace transforms. Initials conditions are optional. 556 559 557 560 INPUT: 558 561 559 562  ``de``  a lambda expression representing the ODE (eg, de = 560 563 diff(y,x,2) == diff(y,x)+sin(x)) 561 564 … … 573 576 Solution of the ODE as symbolic expression 574 577 575 578 EXAMPLES:: 576 579 577 580 sage: u=function('u',x) 578 581 sage: eq = diff(u,x)  exp(x)  u == 0 579 582 sage: desolve_laplace(eq,u) 580 583 1/2*(2*u(0) + 1)*e^x  1/2*e^(x) 581 584 582 585 We can use initial conditions:: 583 586 584 587 sage: desolve_laplace(eq,u,ics=[0,3]) 585 588 1/2*e^(x) + 7/2*e^x 586 589 587 The initial conditions do not persist in the system (as they persisted 590 The initial conditions do not persist in the system (as they persisted 588 591 in previous versions):: 589 592 590 593 sage: desolve_laplace(eq,u) … … 596 599 sage: eq = diff(f,x) + f == 0 597 600 sage: desolve_laplace(eq,f,[0,1]) 598 601 e^(x) 599 602 600 603 :: 601 604 602 605 sage: x = var('x') 603 606 sage: f = function('f', x) 604 607 sage: de = diff(f,x,x)  2*diff(f,x) + f … … 624 627 625 628 sage: soln(t=3) 626 629 e^(3) + 1 627 628 AUTHORS: 630 631 AUTHORS: 629 632 630 633  David Joyner (12006,82007) 631 634 632 635  Robert Marik (102009) 633 636 """ 634 637 #This is the original code from David Joyner (inputs and outputs strings) … … 659 662 ## verbatim copy from desolve  end 660 663 661 664 def sanitize_var(exprs): # 'y(x) > y(x) 662 return exprs.replace("'"+str(dvar),str(dvar)) 665 return exprs.replace("'"+str(dvar),str(dvar)) 663 666 de0=de._maxima_() 664 667 P = de0.parent() 665 668 cmd = sanitize_var("desolve("+de0.str()+","+str(dvar)+")") 666 669 soln=P(cmd).rhs() 667 670 if str(soln).strip() == 'false': 668 raise NotImplementedError, "Maxima was unable to solve this ODE." 671 raise NotImplementedError, "Maxima was unable to solve this ODE." 669 672 soln=soln.sage() 670 673 if ics!=None: 671 674 d = len(ics) 672 675 for i in range(0,d1): 673 soln=eval('soln.substitute(diff(dvar,ivar,i)('+str(ivar)+'=ics[0])==ics[i+1])') 676 soln=eval('soln.substitute(diff(dvar,ivar,i)('+str(ivar)+'=ics[0])==ics[i+1])') 674 677 return soln 675 676 678 679 677 680 def desolve_system(des, vars, ics=None, ivar=None): 678 681 """ 679 682 Solves any size system of 1st order ODE's. Initials conditions are optional. 680 681 Onedimensional systems are passed to :meth:`desolve_laplace`. 683 684 Onedimensional systems are passed to :meth:`desolve_laplace`. 682 685 683 686 INPUT: 684 687 685 688  ``des``  list of ODEs 686 689 687 690  ``vars``  list of dependent variables 688 691 689 692  ``ics``  (optional) list of initial values for ivar and vars 690 693 691 694  ``ivar``  (optional) the independent variable, which must be 692 695 specified if there is more than one independent variable in the 693 696 equation. … … 702 705 sage: desolve_system([de1, de2], [x,y]) 703 706 [x(t) == (x(0)  1)*cos(t)  (y(0)  1)*sin(t) + 1, 704 707 y(t) == (x(0)  1)*sin(t) + (y(0)  1)*cos(t) + 1] 705 708 706 709 Now we give some initial conditions:: 707 710 708 711 sage: sol = desolve_system([de1, de2], [x,y], ics=[0,1,2]); sol 709 712 [x(t) == sin(t) + 1, y(t) == cos(t) + 1] 710 713 711 714 :: 712 715 713 716 sage: solnx, solny = sol[0].rhs(), sol[1].rhs() 714 717 sage: plot([solnx,solny],(0,1)) # not tested 715 718 sage: parametric_plot((solnx,solny),(0,1)) # not tested … … 717 720 TESTS: 718 721 719 722 Trac #9823 fixed:: 720 723 721 724 sage: t = var('t') 722 725 sage: x = function('x', t) 723 726 sage: de1 = diff(x,t) + 1 == 0 724 727 sage: desolve_system([de1], [x]) 725 728 t + x(0) 726 729 727 AUTHORS: 730 AUTHORS: 728 731 729 732  Robert Bradshaw (102008) 730 733 """ … … 756 759 for dvar, ic in zip(dvars, ics[:1]): 757 760 dvar.atvalue(ivar==ivar_ic, dvar) 758 761 return soln 759 762 760 763 761 764 def desolve_system_strings(des,vars,ics=None): 762 765 r""" … … 788 791 sage: s = var('s') 789 792 sage: function('x', s) 790 793 x(s) 791 794 792 795 :: 793 796 794 797 sage: function('y', s) 795 798 y(s) 796 799 797 800 :: 798 801 799 802 sage: de1 = lambda z: diff(z[0],s) + z[1]  1 800 803 sage: de2 = lambda z: diff(z[1],s)  z[0] + 1 801 804 sage: des = [de1([x(s),y(s)]),de2([x(s),y(s)])] 802 805 sage: vars = ["s","x","y"] 803 806 sage: desolve_system_strings(des,vars) 804 ["(1'y(0))*sin(s)+('x(0)1)*cos(s)+1", "('x(0)1)*sin(s)+('y(0)1)*cos(s)+1"] 805 807 ["(1'y(0))*sin(s)+('x(0)1)*cos(s)+1", 808 "('x(0)1)*sin(s)+('y(0)1)*cos(s)+1"] 809 806 810 :: 807 811 808 812 sage: ics = [0,1,1] 809 813 sage: soln = desolve_system_strings(des,vars,ics); soln 810 814 ['2*sin(s)+1', '12*cos(s)'] 811 815 812 816 :: 813 817 814 818 sage: solnx, solny = map(SR, soln) 815 819 sage: RR(solnx(s=3)) 816 820 1.28224001611973 817 821 818 822 :: 819 823 820 824 sage: P1 = plot([solnx,solny],(0,1)) 821 825 sage: P2 = parametric_plot((solnx,solny),(0,1)) 822 826 823 827 Now type show(P1), show(P2) to view these. 824 825 828 826 AUTHORS: 829 830 AUTHORS: 827 831 828 832  David Joyner (32006, 82007) 829 833 """ … … 856 860 last column. 857 861 858 862 *For pedagogical purposes only.* 859 863 860 864 EXAMPLES:: 861 865 862 866 sage: from sage.calculus.desolvers import eulers_method 863 867 sage: x,y = PolynomialRing(QQ,2,"xy").gens() 864 868 sage: eulers_method(5*x+y5,0,1,1/2,1) … … 912 916 sage: P2 = line(pts) 913 917 sage: (P1+P2).show() 914 918 915 AUTHORS: 919 AUTHORS: 916 920 917 921  David Joyner 918 922 """ … … 948 952 next (last) column. 949 953 950 954 *For pedagogical purposes only.* 951 955 952 956 EXAMPLES:: 953 957 954 958 sage: from sage.calculus.desolvers import eulers_method_2x2 955 959 sage: t, x, y = PolynomialRing(QQ,3,"txy").gens() 956 960 sage: f = x+y+t; g = xy 957 961 sage: eulers_method_2x2(f,g, 0, 0, 0, 1/3, 1,algorithm="none") 958 [[0, 0, 0], [1/3, 0, 0], [2/3, 1/9, 0], [1, 10/27, 1/27], [4/3, 68/81, 4/27]] 962 [[0, 0, 0], [1/3, 0, 0], [2/3, 1/9, 0], [1, 10/27, 1/27], 963 [4/3, 68/81, 4/27]] 959 964 960 965 :: 961 966 … … 994 999 3/4 0.63 0.0078 0.031 0.11 995 1000 1 0.63 0.020 0.079 0.071 996 1001 997 To numerically approximate y(1), where `y''+ty'+y=0`, `y(0)=1`, `y'(0)=0`:: 1002 To numerically approximate y(1), where `y''+ty'+y=0`, `y(0)=1`, `y'(0)=0`:: 998 1003 999 1004 sage: t,x,y=PolynomialRing(RR,3,"txy").gens() 1000 1005 sage: f = y; g = xy*t … … 1006 1011 3/4 0.88 0.15 0.62 0.10 1007 1012 1 0.75 0.17 0.68 0.015 1008 1013 1009 AUTHORS: 1014 AUTHORS: 1010 1015 1011 1016  David Joyner 1012 1017 """ … … 1036 1041 `x(a)=x_0`, `y' = g(t,x,y)`, `y(a) = y_0`. 1037 1042 1038 1043 *For pedagogical purposes only.* 1039 1044 1040 1045 EXAMPLES:: 1041 1046 1042 1047 sage: from sage.calculus.desolvers import eulers_method_2x2_plot … … 1066 1071 def desolve_rk4_determine_bounds(ics,end_points=None): 1067 1072 """ 1068 1073 Used to determine bounds for numerical integration. 1069 1074 1070 1075  If end_points is None, the interval for integration is from ics[0] 1071 1076 to ics[0]+10 1072 1077 1073  If end_points is a or [a], the interval for integration is from min(ics[0],a)1074 to max(ics[0],a)1078  If end_points is a or [a], the interval for integration is from 1079 min(ics[0],a) to max(ics[0],a) 1075 1080 1076 1081  If end_points is [a,b], the interval for integration is from min(ics[0],a) 1077 1082 to max(ics[0],b) … … 1081 1086 sage: from sage.calculus.desolvers import desolve_rk4_determine_bounds 1082 1087 sage: desolve_rk4_determine_bounds([0,2],1) 1083 1088 (0, 1) 1084 1089 1085 1090 :: 1086 1091 1087 sage: desolve_rk4_determine_bounds([0,2]) 1092 sage: desolve_rk4_determine_bounds([0,2]) 1088 1093 (0, 10) 1089 1094 1090 1095 :: 1091 1096 1092 1097 sage: desolve_rk4_determine_bounds([0,2],[2]) 1093 1098 (2, 0) 1094 1099 1095 1100 :: 1096 1101 1097 1102 sage: desolve_rk4_determine_bounds([0,2],[2,4]) 1098 1103 (2, 4) 1099 1104 1100 1105 """ 1101 if end_points is None: 1106 if end_points is None: 1102 1107 return((ics[0],ics[0]+10)) 1103 1108 if not isinstance(end_points,list): 1104 1109 end_points=[end_points] … … 1106 1111 return (min(ics[0],end_points[0]),max(ics[0],end_points[0])) 1107 1112 else: 1108 1113 return (min(ics[0],end_points[0]),max(ics[0],end_points[1])) 1109 1114 1110 1115 1111 1116 def desolve_rk4(de, dvar, ics=None, ivar=None, end_points=None, step=0.1, output='list', **kwds): 1112 1117 """ … … 1114 1119 equation. See also ``ode_solver``. 1115 1120 1116 1121 INPUT: 1117 1118 input is similar to ``desolve`` command. The differential equation can be 1122 1123 input is similar to ``desolve`` command. The differential equation can be 1119 1124 written in a form close to the plot_slope_field or desolve command 1120 1125 1121 1126  Variant 1 (function in two variables) 1122 1123  ``de``  right hand side, i.e. the function `f(x,y)` from ODE `y'=f(x,y)` 1124 1127 1128  ``de``  right hand side, i.e. the function `f(x,y)` from ODE 1129 `y'=f(x,y)` 1130 1125 1131  ``dvar``  dependent variable (symbolic variable declared by var) 1126 1132 1127 1133  Variant 2 (symbolic equation) 1128 1134 1129 1135  ``de``  equation, including term with ``diff(y,x)`` 1130 1136 1131  ``dvar```  dependent variable (declared as funciton of independent variable) 1137  ``dvar```  dependent variable (declared as a function of independent 1138 variables) 1132 1139 1133  Other parameters 1140  Other parameters 1134 1141 1135  ``ivar``  should be specified, if there are more variables or if the equation is autonomous 1142  ``ivar``  should be specified, if there are more variables or if the 1143 equation is autonomous 1136 1144 1137 1145  ``ics``  initial conditions in the form [x0,y0] 1138 1146 1139 1147  ``end_points``  the end points of the interval 1140 1148 1141  if end_points is a or [a], we integrate on between min(ics[0],a) and max(ics[0],a) 1149  if end_points is a or [a], we integrate on between min(ics[0],a) and 1150 max(ics[0],a) 1151 1142 1152  if end_points is None, we use end_points=ics[0]+10 1143 1153 1144  if end_points is [a,b] we integrate on between min(ics[0],a) and max(ics[0],b) 1154  if end_points is [a,b] we integrate on between min(ics[0],a) and 1155 max(ics[0],b) 1145 1156 1146  ``step``  (optional, default:0.1) the length of the step (positive number) 1147 1157  ``step``  (optional, default:0.1) the length of the step (positive 1158 number) 1159 1148 1160  ``output``  (optional, default: 'list') one of 'list', 1149 1161 'plot', 'slope_field' (graph of the solution with slope field) 1150 1162 1151 OUTPUT: 1163 OUTPUT: 1152 1164 1153 1165 Returns a list of points, or plot produced by list_plot, 1154 1166 optionally with slope field. … … 1166 1178 1167 1179 Variant 1 for input  we can pass ODE in the form used by 1168 1180 desolve function In this example we integrate bakwards, since 1169 ``end_points < ics[0]``:: 1181 ``end_points < ics[0]``:: 1170 1182 1171 sage: y=function('y',x) 1172 sage: desolve_rk4(diff(y,x)+y*(y1) == x2,y,ics=[1,1],step=0.5, end_points=0) 1183 sage: y=function('y',x) 1184 sage: desolve_rk4(diff(y,x)+y*(y1) == x2,y,ics=[1,1],step=0.5, end_points=0) 1173 1185 [[0.0, 8.90425710896], [0.5, 1.90932794536], [1, 1]] 1174 1186 1175 1187 Here we show how to plot simple pictures. For more advanced … … 1178 1190 1179 1191 sage: x,y=var('x y') 1180 1192 sage: P=desolve_rk4(y*(2y),y,ics=[0,.1],ivar=x,output='slope_field',end_points=[4,6],thickness=3) 1181 1193 1182 1194 ALGORITHM: 1183 1195 1184 1196 4th order RungeKutta method. Wrapper for command ``rk`` in 1185 1197 Maxima's dynamics package. Perhaps could be faster by using 1186 1198 fast_float instead. 1187 1199 1188 AUTHORS: 1200 AUTHORS: 1189 1201 1190 1202  Robert Marik (102009) 1191 1203 """ 1192 if ics is None: 1204 if ics is None: 1193 1205 raise ValueError, "No initial conditions, specify with ics=[x0,y0]." 1194 1206 1195 1207 if ivar is None: … … 1246 1258 YMIN=sol[0][1] 1247 1259 XMAX=XMIN 1248 1260 YMAX=YMIN 1249 for s,t in sol: 1261 for s,t in sol: 1250 1262 if s>XMAX:XMAX=s 1251 1263 if s<XMIN:XMIN=s 1252 1264 if t>YMAX:YMAX=t … … 1260 1272 r""" 1261 1273 Solves numerically system of firstorder ordinary differential 1262 1274 equations using the 4th order RungeKutta method. Wrapper for 1263 Maxima command ``rk``. See also ``ode_solver``. 1275 Maxima command ``rk``. See also ``ode_solver``. 1264 1276 1265 1277 INPUT: 1266 1278 1267 1279 input is similar to desolve_system and desolve_rk4 commands 1268 1280 1269 1281  ``des``  right hand sides of the system 1270 1282 1271 1283  ``vars``  dependent variables … … 1275 1287 missing 1276 1288 1277 1289  ``ics``  initial conditions in the form [x0,y01,y02,y03,....] 1278 1290 1279 1291  ``end_points``  the end points of the interval 1280 1281  if end_points is a or [a], we integrate on between min(ics[0],a) and max(ics[0],a) 1292 1293  if end_points is a or [a], we integrate on between min(ics[0],a) and 1294 max(ics[0],a) 1295 1282 1296  if end_points is None, we use end_points=ics[0]+10 1283 1297 1284  if end_points is [a,b] we integrate on between min(ics[0],a) and max(ics[0],b) 1298  if end_points is [a,b] we integrate on between min(ics[0],a) and 1299 max(ics[0],b) 1285 1300 1286 1301  ``step``  (optional, default: 0.1) the length of the step 1287 1302 … … 1300 1315 sage: P=desolve_system_rk4([x*(1y),y*(1x)],[x,y],ics=[0,0.5,2],ivar=t,end_points=20) 1301 1316 sage: Q=[ [i,j] for i,j,k in P] 1302 1317 sage: LP=list_plot(Q) 1303 1318 1304 1319 sage: Q=[ [j,k] for i,j,k in P] 1305 1320 sage: LP=list_plot(Q) 1306 1321 1307 1322 ALGORITHM: 1308 1323 1309 1324 4th order RungeKutta method. Wrapper for command ``rk`` in Maxima's 1310 1325 dynamics package. Perhaps could be faster by using ``fast_float`` 1311 1326 instead. 1312 1327 1313 AUTHOR: 1328 AUTHOR: 1314 1329 1315 1330  Robert Marik (102009) 1316 1331 """ 1317 1332 1318 if ics is None: 1333 if ics is None: 1319 1334 raise ValueError, "No initial conditions, specify with ics=[x0,y01,y02,...]." 1320 1335 1321 1336 ivars = set([]) … … 1335 1350 x0=ics[0] 1336 1351 icss = [ics[i]._maxima_().str() for i in range(1,len(ics))] 1337 1352 icstr = "[" + ",".join(icss) + "]" 1338 step=abs(step) 1353 step=abs(step) 1339 1354 1340 1355 maxima("load('dynamics)") 1341 1356 lower_bound,upper_bound=desolve_rk4_determine_bounds(ics,end_points) … … 1361 1376 , rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0 1362 1377 , mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0): 1363 1378 r""" 1364 Solves numerically a system of firstorder ordinary differential equations 1379 Solves numerically a system of firstorder ordinary differential equations 1365 1380 using ``odeint`` from scipy.integrate module. 1366 1381 1367 1382 INPUT: 1368 1383 1369 1384  ``des``  right hand sides of the system 1370 1385 1371  ``ics``  initial conditions 1386  ``ics``  initial conditions 1372 1387 1373 1388  ``times``  a sequence of time points in which the solution must be found 1374 1389 … … 1377 1392 1378 1393  ``ivar``  independent variable, optional. 1379 1394 1380  ``compute_jac``  boolean. If True, the Jacobian of des is computed and 1395  ``compute_jac``  boolean. If True, the Jacobian of des is computed and 1381 1396 used during the integration of Stiff Systems. Default value is False. 1382 1397 1383 Other Parameters (taken from the documentation of odeint function from 1398 Other Parameters (taken from the documentation of odeint function from 1384 1399 scipy.integrate module) 1385 1400 1386 1401  ``rtol``, ``atol`` : float … … 1410 1425  ``hmin`` : float, (0: solverdetermined) 1411 1426 The minimum absolute step size allowed. 1412 1427 1413  ``ixpr`` : boolean. 1428  ``ixpr`` : boolean. 1414 1429 Whether to generate extra printing at method switches. 1415 1430 1416 1431  ``mxstep`` : integer, (0: solverdetermined) … … 1425 1440 1426 1441  ``mxords`` : integer, (0: solverdetermined) 1427 1442 Maximum order to be allowed for the stiff (BDF) method. 1428 1443 1429 1444 OUTPUT: 1430 1445 1431 1446 Returns a list with the solution of the system at each time in times. 1432 1447 1433 1448 EXAMPLES: 1434 1449 1435 1450 Lotka Volterra Equations:: 1436 1451 1437 1452 sage: from sage.calculus.desolvers import desolve_odeint … … 1440 1455 sage: sol=desolve_odeint(f,[0.5,2],srange(0,10,0.1),[x,y]) 1441 1456 sage: p=line(zip(sol[:,0],sol[:,1])) 1442 1457 sage: p.show() 1443 1458 1444 1459 Lorenz Equations:: 1445 1460 1446 1461 sage: x,y,z=var('x,y,z') … … 1449 1464 sage: rho=28 1450 1465 sage: beta=8/3 1451 1466 sage: # The Lorenz equations 1452 sage: lorenz=[sigma*(yx),x*(rhoz)y,x*ybeta*z] 1467 sage: lorenz=[sigma*(yx),x*(rhoz)y,x*ybeta*z] 1453 1468 sage: # Time and initial conditions 1454 1469 sage: times=srange(0,50.05,0.05) 1455 1470 sage: ics=[0,1,1] 1456 1471 sage: sol=desolve_odeint(lorenz,ics,times,[x,y,z],rtol=1e13,atol=1e14) 1457 1472 1458 1473 Onedimensional Stiff system:: 1459 1474 1460 1475 sage: y= var('y') … … 1465 1480 sage: sol=desolve_odeint(f,ic,t,y,rtol=1e9,atol=1e10,compute_jac=True) 1466 1481 sage: p=points(zip(t,sol)) 1467 1482 sage: p.show() 1468 1469 Another Stiff system with some optional parameters with no 1483 1484 Another Stiff system with some optional parameters with no 1470 1485 default value:: 1471 1486 1472 1487 sage: y1,y2,y3=var('y1,y2,y3') … … 1478 1493 sage: t=srange(0,10,0.01) 1479 1494 sage: v=[y1,y2,y3] 1480 1495 sage: sol=desolve_odeint(f,ci,t,v,rtol=1e3,atol=1e4,h0=0.1,hmax=1,hmin=1e4,mxstep=1000,mxords=17) 1481 1496 1482 1497 AUTHOR: 1483 1498 1484 1499  Oriol Castejon (052010) 1485 1500 """ 1486 1501 … … 1547 1562 def Dfun(y,t): 1548 1563 v = list(y[:]) 1549 1564 v.append(t) 1550 return [[element(*v) for element in row] for row in J] 1551 1565 return [[element(*v) for element in row] for row in J] 1566 1552 1567 sol=odeint(func, ics, times, args=args, Dfun=Dfun, rtol=rtol, atol=atol, 1553 1568 tcrit=tcrit, h0=h0, hmax=hmax, hmin=hmin, ixpr=ixpr, mxstep=mxstep, 1554 1569 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 b 62 62 0 == 10*a + b  144 63 63 64 64 Subtract two symbolic equations:: 65 65 66 66 sage: var('a,b') 67 (a, b) 67 (a, b) 68 68 sage: m = 144 == 20 * a + b 69 69 sage: n = 136 == 10 * a + b 70 70 sage: m  n … … 188 188 sage: assume(x>0, y < 2) 189 189 sage: assumptions() 190 190 [x > 0, y < 2] 191 sage: (y < 2).forget() 191 sage: (y < 2).forget() 192 192 sage: assumptions() 193 193 [x > 0] 194 194 sage: forget() … … 222 222 sage: eq = (x == 2) 223 223 sage: eq._maple_init_() 224 224 'x = 2' 225 225 226 226 Comparison:: 227 227 228 228 sage: x = var('x') … … 262 262 263 263 sage: f = x^5 + a 264 264 sage: solve(f==0,x) 265 [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)] 265 [x == (a)^(1/5)*e^(2/5*I*pi), x == (a)^(1/5)*e^(4/5*I*pi), 266 x == (a)^(1/5)*e^(4/5*I*pi), x == (a)^(1/5)*e^(2/5*I*pi), 267 x == (a)^(1/5)] 266 268 267 269 You can also do arithmetic with inequalities, as illustrated 268 270 below:: … … 297 299 sage: eqn = x^3 + 2/3 >= x 298 300 sage: loads(dumps(eqn)) 299 301 x^3 + 2/3 >= x 300 sage: loads(dumps(eqn)) == eqn 302 sage: loads(dumps(eqn)) == eqn 301 303 True 302 304 303 305 AUTHORS: … … 433 435 Used internally by the symbolic solve command to convert the output 434 436 of Maxima's solve command to a list of solutions in Sage's symbolic 435 437 package. 436 438 437 439 EXAMPLES: 438 440 439 441 We derive the (monic) quadratic formula:: 440 442 441 443 sage: var('x,a,b') 442 444 (x, a, b) 443 445 sage: solve(x^2 + a*x + b == 0, x) 444 446 [x == 1/2*a  1/2*sqrt(a^2  4*b), x == 1/2*a + 1/2*sqrt(a^2  4*b)] 445 447 446 448 Behind the scenes when the above is evaluated the function 447 449 :func:`string_to_list_of_solutions` is called with input the 448 450 string `s` below:: 449 451 450 452 sage: s = '[x=(sqrt(a^24*b)+a)/2,x=(sqrt(a^24*b)a)/2]' 451 453 sage: sage.symbolic.relation.string_to_list_of_solutions(s) 452 454 [x == 1/2*a  1/2*sqrt(a^2  4*b), x == 1/2*a + 1/2*sqrt(a^2  4*b)] … … 464 466 def solve(f, *args, **kwds): 465 467 r""" 466 468 Algebraically solve an equation or system of equations (over the 467 complex numbers) for given variables. Inequalities and systems 468 of inequalities are also supported. 469 469 complex numbers) for given variables. Inequalities and systems 470 of inequalities are also supported. 471 470 472 INPUT: 471 473 472 474  ``f``  equation or system of equations (given by a 473 475 list or tuple) 474 476 475 477  ``*args``  variables to solve for. 476 477  ``solution_dict``  bool (default: False); if True or nonzero, 478 479  ``solution_dict``  bool (default: False); if True or nonzero, 478 480 return a list of dictionaries containing the solutions. If there 479 481 are no solutions, return an empty list (rather than a list containing 480 482 an empty dictionary). Likewise, if there's only a single solution, 481 483 return a list containing one dictionary with that solution. 482 484 483 485 EXAMPLES:: 484 486 485 487 sage: x, y = var('x, y') 486 488 sage: solve([x+y==6, xy==4], x, y) 487 489 [[x == 5, y == 1]] 488 490 sage: solve([x^2+y^2 == 1, y^2 == x^3 + x + 1], x, y) 489 [[x == 1/2*I*sqrt(3)  1/2, y == sqrt(1/2*I*sqrt(3) + 3/2)], 490 [x == 1/2*I*sqrt(3)  1/2, y == sqrt(1/2*I*sqrt(3) + 3/2)], 491 [x == 1/2*I*sqrt(3)  1/2, y == sqrt(1/2*I*sqrt(3) + 3/2)], 492 [x == 1/2*I*sqrt(3)  1/2, y == sqrt(1/2*I*sqrt(3) + 3/2)], 493 [x == 0, y == 1], 491 [[x == 1/2*I*sqrt(3)  1/2, y == sqrt(1/2*I*sqrt(3) + 3/2)], 492 [x == 1/2*I*sqrt(3)  1/2, y == sqrt(1/2*I*sqrt(3) + 3/2)], 493 [x == 1/2*I*sqrt(3)  1/2, y == sqrt(1/2*I*sqrt(3) + 3/2)], 494 [x == 1/2*I*sqrt(3)  1/2, y == sqrt(1/2*I*sqrt(3) + 3/2)], 495 [x == 0, y == 1], 494 496 [x == 0, y == 1]] 495 497 sage: solve([sqrt(x) + sqrt(y) == 5, x + y == 10], x, y) 496 [[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]] 498 [[x == 5/2*I*sqrt(5) + 5, y == 5/2*I*sqrt(5) + 5], 499 [x == 5/2*I*sqrt(5) + 5, y == 5/2*I*sqrt(5) + 5]] 497 500 sage: solutions=solve([x^2+y^2 == 1, y^2 == x^3 + x + 1], x, y, solution_dict=True) 498 501 sage: for solution in solutions: print solution[x].n(digits=3), ",", solution[y].n(digits=3) 499 502 0.500  0.866*I , 1.27 + 0.341*I … … 503 506 0.000 , 1.00 504 507 0.000 , 1.00 505 508 506 Whenever possible, answers will be symbolic, but with systems of 507 equations, at times approximations will be given, due to the 509 Whenever possible, answers will be symbolic, but with systems of 510 equations, at times approximations will be given, due to the 508 511 underlying algorithm in Maxima:: 509 512 510 513 sage: sols = solve([x^3==y,y^2==x],[x,y]); sols[1], sols[0] 511 ([x == 0, y == 0], [x == (0.309016994375 + 0.951056516295*I), y == (0.809016994375  0.587785252292*I)]) 514 ([x == 0, y == 0], 515 [x == (0.309016994375 + 0.951056516295*I), 516 y == (0.809016994375  0.587785252292*I)]) 512 517 sage: sols[0][0].rhs().pyobject().parent() 513 518 Complex Double Field 514 519 … … 516 521 for symbolic expressions, which defaults to exact answers only:: 517 522 518 523 sage: solve([y^6==y],y) 519 [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] 524 [y == e^(2/5*I*pi), y == e^(4/5*I*pi), y == e^(4/5*I*pi), 525 y == e^(2/5*I*pi), y == 1, y == 0] 520 526 sage: solve( [y^6 == y], y)==solve( y^6 == y, y) 521 527 True 522 528 … … 540 546 TypeError: 5 is not a valid variable. 541 547 542 548 If we ask for dictionaries containing the solutions, we get them:: 543 549 544 550 sage: solve([x^21],x,solution_dict=True) 545 551 [{x: 1}, {x: 1}] 546 552 sage: solve([x^24*x+4],x,solution_dict=True) … … 550 556 x: 2, y: 4 551 557 x: 2, y: 4 552 558 553 If there is a parameter in the answer, that will show up as 554 a new variable. In the following example, ``r1`` is a real free555 variable(because of the ``r``)::559 If there is a parameter in the answer, that will show up as a new variable. 560 In the following example, ``r1`` is a real free variable 561 (because of the ``r``):: 556 562 557 563 sage: solve([x+y == 3, 2*x+2*y == 6],x,y) 558 564 [[x == r1 + 3, y == r1]] … … 574 580 solutions are returned. E.g., note that the first 575 581 ``3==3`` evaluates to ``True``, not to a 576 582 symbolic equation. 577 583 578 584 :: 579 585 580 586 sage: solve([3==3, 1.00000000000000*x^3 == 0], x) 581 587 [x == 0] 582 588 sage: solve([1.00000000000000*x^3 == 0], x) 583 589 [x == 0] 584 590 585 591 Here, the first equation evaluates to ``False``, so 586 592 there are no solutions:: 587 593 588 594 sage: solve([1==3, 1.00000000000000*x^3 == 0], x) 589 595 [] 590 596 591 597 Completely symbolic solutions are supported:: 592 598 593 599 sage: var('s,j,b,m,g') 594 600 (s, j, b, m, g) 595 601 sage: sys = [ m*(1s)  b*s*j, b*s*jg*j ]; … … 608 614 Use use_grobner if no solution is obtained from to_poly_solve:: 609 615 610 616 sage: x,y=var('x y'); c1(x,y)=(x5)^2+y^216; c2(x,y)=(y3)^2+x^29 611 sage: solve([c1(x,y),c2(x,y)],[x,y]) 612 [[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]] 613 617 sage: solve([c1(x,y),c2(x,y)],[x,y]) 618 [[x == 9/68*sqrt(55) + 135/68, y == 15/68*sqrt(5)*sqrt(11) + 123/68], 619 [x == 9/68*sqrt(55) + 135/68, y == 15/68*sqrt(5)*sqrt(11) + 123/68]] 620 614 621 TESTS:: 615 622 616 623 sage: solve([sin(x)==x,y^2==x],x,y) … … 672 679 variables = args 673 680 else: 674 681 variables = tuple(args[0]) 675 682 676 683 for v in variables: 677 684 if not is_SymbolicVariable(v): 678 685 raise TypeError, "%s is not a valid variable."%v … … 732 739 Return all solutions to an equation or list of equations modulo the 733 740 given integer modulus. Each equation must involve only polynomials 734 741 in 1 or many variables. 735 742 736 743 By default the solutions are returned as `n`tuples, where `n` 737 744 is the number of variables appearing anywhere in the given 738 745 equations. The variables are in alphabetical order. 739 746 740 747 INPUT: 741 742 748 749 743 750  ``eqns``  equation or list of equations 744 751 745 752  ``modulus``  an integer 746 753 747 754  ``solution_dict``  bool (default: False); if True or nonzero, 748 755 return a list of dictionaries containing the solutions. If there 749 756 are no solutions, return an empty list (rather than a list containing 750 757 an empty dictionary). Likewise, if there's only a single solution, 751 758 return a list containing one dictionary with that solution. 752 759 753 760 754 761 EXAMPLES:: 755 762 756 763 sage: var('x,y') 757 764 (x, y) 758 765 sage: solve_mod([x^2 + 2 == x, x^2 + y == y^2], 14) 759 766 [(4, 2), (4, 6), (4, 9), (4, 13)] 760 767 sage: solve_mod([x^2 == 1, 4*x == 11], 15) 761 768 [(14,)] 762 769 763 770 Fermat's equation modulo 3 with exponent 5:: 764 771 765 772 sage: var('x,y,z') 766 773 (x, y, z) 767 774 sage: solve_mod([x^5 + y^5 == z^5], 3) 768 [(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)] 775 [(0, 0, 0), (0, 1, 1), (0, 2, 2), (1, 0, 1), (1, 1, 2), (1, 2, 0), 776 (2, 0, 2), (2, 1, 0), (2, 2, 1)] 769 777 770 We can solve with respect to a bigger modulus if it consists only of small prime factors:: 778 We can solve with respect to a bigger modulus if it consists only of small 779 prime factors:: 771 780 772 781 sage: [d] = solve_mod([5*x + y == 3, 2*x  3*y == 9], 3*5*7*11*19*23*29, solution_dict = True) 773 782 sage: d[x] … … 780 789 is large:: 781 790 782 791 sage: sorted(solve_mod([x^2 == 41], 10^20)) 783 [(4538602480526452429,), (11445932736758703821,), (38554067263241296179,), 784 (45461397519473547571,), (54538602480526452429,), (61445932736758703821,), 785 (88554067263241296179,), (95461397519473547571,)] 792 [(4538602480526452429,), (11445932736758703821,), 793 (38554067263241296179,), (45461397519473547571,), 794 (54538602480526452429,), (61445932736758703821,), 795 (88554067263241296179,), (95461397519473547571,)] 786 796 787 797 We solve a simple equation modulo 2:: 788 798 789 799 sage: x,y = var('x,y') 790 800 sage: solve_mod([x == y], 2) 791 801 [(0, 0), (1, 1)] … … 821 831 [(8, 5, 6)] 822 832 823 833 Confirm that modulus 1 now behaves as it should:: 824 834 825 835 sage: x, y = var("x y") 826 836 sage: solve_mod([x==1], 1) 827 837 [(0,)] 828 838 sage: solve_mod([2*x^2+x*y, x*y+2*y^2+x2*y, 2*x^2+2*x*yy^2xy], 1) 829 839 [(0, 0)] 830 840 831 841 832 842 """ 833 843 from sage.rings.all import Integer, Integers, crt_basis … … 844 854 raise ValueError, "the modulus must be a positive integer" 845 855 vars = list(set(sum([list(e.variables()) for e in eqns], []))) 846 856 vars.sort(key=repr) 847 857 848 858 if modulus == 1: # degenerate case 849 859 ans = [tuple(Integers(1)(0) for v in vars)] 850 860 return ans … … 852 862 factors = modulus.factor() 853 863 crt_basis = vector(Integers(modulus), crt_basis([p**i for p,i in factors])) 854 864 solutions = [] 855 865 856 866 has_solution = True 857 867 for p,i in factors: 858 868 solution =_solve_mod_prime_power(eqns, p, i, vars) … … 861 871 else: 862 872 has_solution = False 863 873 break 864 874 865 875 866 876 ans = [] 867 877 if has_solution: … … 882 892 Internal help function for solve_mod, does little checking since it expects 883 893 solve_mod to do that 884 894 885 Return all solutions to an equation or list of equations modulo p^m. 895 Return all solutions to an equation or list of equations modulo p^m. 886 896 Each equation must involve only polynomials 887 897 in 1 or many variables. 888 898 889 899 The solutions are returned as `n`tuples, where `n` 890 900 is the number of variables in vars. 891 901 892 902 INPUT: 893 894 903 904 895 905  ``eqns``  equation or list of equations 896 906 897 907  ``p``  a prime 898 908 899 909  ``i``  an integer > 0 900 910 901 911  ``vars``  a list of variables to solve for 902 903 912 913 904 914 EXAMPLES:: 905 915 906 916 sage: var('x,y') 907 917 (x, y) 908 918 sage: solve_mod([x^2 + 2 == x, x^2 + y == y^2], 14) 909 919 [(4, 2), (4, 6), (4, 9), (4, 13)] 910 920 sage: solve_mod([x^2 == 1, 4*x == 11], 15) 911 921 [(14,)] 912 922 913 923 Fermat's equation modulo 3 with exponent 5:: 914 924 915 925 sage: var('x,y,z') 916 926 (x, y, z) 917 927 sage: solve_mod([x^5 + y^5 == z^5], 3) 918 [(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)] 928 [(0, 0, 0), (0, 1, 1), (0, 2, 2), (1, 0, 1), (1, 1, 2), (1, 2, 0), 929 (2, 0, 2), (2, 1, 0), (2, 2, 1)] 919 930 920 931 We solve a simple equation modulo 2:: 921 932 922 933 sage: x,y = var('x,y') 923 934 sage: solve_mod([x == y], 2) 924 935 [(0, 0), (1, 1)] 925 936 926 937 927 938 .. warning:: 928 939 929 940 Currently this constructs possible solutions by building up 930 941 from the smallest prime factor of the modulus. The interface 931 942 is good, but the algorithm is horrible if the modulus isn't the … … 937 948 interface. 938 949 939 950 TESTS: 940 951 941 952 Confirm we can reproduce the first few terms of OEIS A187719:: 942 953 943 954 sage: from sage.symbolic.relation import _solve_mod_prime_power 944 955 sage: [sorted(_solve_mod_prime_power([x^2==41], 10, i, [x]))[0][0] for i in [1..13]] 945 [1, 21, 71, 1179, 2429, 47571, 1296179, 8703821, 26452429, 526452429, 946 13241296179, 19473547571, 2263241296179]956 [1, 21, 71, 1179, 2429, 47571, 1296179, 8703821, 26452429, 526452429, 957 13241296179, 19473547571, 2263241296179] 947 958 948 959 """ 949 960 from sage.rings.all import Integers, PolynomialRing 950 961 from sage.modules.all import vector 951 962 from sage.misc.all import cartesian_product_iterator 952 963 953 964 mrunning = 1 954 965 ans = [] 955 966 for mi in xrange(m): … … 970 981 971 982 def solve_ineq_univar(ineq): 972 983 """ 973 Function solves rational inequality in one variable. 974 984 Function solves rational inequality in one variable. 985 975 986 INPUT: 976 987 977 988  ``ineq``  inequality in one variable 978 989 979 990 OUTPUT: 980 991 981 992  ``list``  output is list of solutions as a list of simple inequalities 982 output [A,B,C] means (A or B or C) each A, B, C is again a list and 983 if A=[a,b], then A means (a and b). The list is empty if there is no 993 output [A,B,C] means (A or B or C) each A, B, C is again a list and 994 if A=[a,b], then A means (a and b). The list is empty if there is no 984 995 solution. 985 996 986 997 EXAMPLES:: 987 998 988 999 sage: from sage.symbolic.relation import solve_ineq_univar 989 1000 sage: solve_ineq_univar(x1/x>0) 990 1001 [[x > 1, x < 0], [x > 1]] 991 1002 992 1003 sage: solve_ineq_univar(x^21/x>0) 993 1004 [[x < 0], [x > 1]] 994 1005 995 1006 sage: solve_ineq_univar((x^31)*x<=0) 996 1007 [[x >= 0, x <= 1]] 997 1008 998 1009 ALGORITHM: 999 1010 1000 1011 Calls Maxima command solve_rat_ineq 1001 1012 1002 1013 AUTHORS: 1003 1014 1004 1015  Robert Marik (012010) 1005 1016 """ 1006 1017 ineqvar = ineq.variables() … … 1009 1020 ineq0 = ineq._maxima_() 1010 1021 ineq0.parent().eval("if solve_rat_ineq_loaded#true then (solve_rat_ineq_loaded:true,load(\"solve_rat_ineq.mac\")) ") 1011 1022 sol = ineq0.solve_rat_ineq().sage() 1012 if repr(sol)=="all": 1013 from sage.rings.infinity import Infinity 1014 sol = [ineqvar[0]<Infinity] 1023 if repr(sol)=="all": 1024 from sage.rings.infinity import Infinity 1025 sol = [ineqvar[0]<Infinity] 1015 1026 return sol 1016 1027 1017 1028 def solve_ineq_fourier(ineq,vars=None): 1018 1029 """ 1019 Solves sytem of inequalities using Maxima and fourier elimination 1030 Solves sytem of inequalities using Maxima and fourier elimination 1020 1031 1021 1032 Can be used for system of linear inequalities and for some types 1022 1033 of nonlinear inequalities. For examples see the section EXAMPLES 1023 1034 below and http://maxima.cvs.sourceforge.net/viewvc/maxima/maxima/share/contrib/fourier_elim/rtest_fourier_elim.mac 1024 1025 1035 1036 1026 1037 INPUT: 1027 1028  ``ineq``  list with system of inequalities1029 1038 1030  ``vars``  optionally list with variables for fourier elimination. 1039  ``ineq``  list with system of inequalities 1040 1041  ``vars``  optionally list with variables for fourier elimination. 1031 1042 1032 1043 OUTPUT: 1033 1034  ``list``  output is list of solutions as a list of simple inequalities 1035 output [A,B,C] means (A or B or C) each A, B, C is again a list and 1036 if A=[a,b], then A means (a and b). The list is empty if there is no 1044 1045  ``list``  output is list of solutions as a list of simple inequalities 1046 output [A,B,C] means (A or B or C) each A, B, C is again a list and 1047 if A=[a,b], then A means (a and b). The list is empty if there is no 1037 1048 solution. 1038 1049 1039 1050 EXAMPLES:: 1040 1051 1041 1052 sage: from sage.symbolic.relation import solve_ineq_fourier 1042 sage: y=var('y') 1053 sage: y=var('y') 1043 1054 sage: solve_ineq_fourier([x+y<9,xy>4],[x,y]) 1044 [[y + 4 < x, x < y + 9, y < (5/2)]] 1055 [[y + 4 < x, x < y + 9, y < (5/2)]] 1045 1056 sage: solve_ineq_fourier([x+y<9,xy>4],[y,x]) 1046 1057 [[y < min(x  4, x + 9)]] 1047 1058 1048 1059 sage: solve_ineq_fourier([x^2>=0]) 1049 1060 [[x < +Infinity]] 1050 1061 … … 1061 1072 [[y < x, 0 < y]] 1062 1073 1063 1074 ALGORITHM: 1064 1075 1065 1076 Calls Maxima command fourier_elim 1066 1077 1067 1078 AUTHORS: 1068 1079 1069 1080  Robert Marik (012010) 1070 1081 """ 1071 1082 if vars is None: … … 1080 1091 if not atom(x) and op(x)=\"or\" then args(x) \ 1081 1092 else [x]") 1082 1093 sol = sol.or_to_list().sage() 1083 if repr(sol) == "[emptyset]": 1094 if repr(sol) == "[emptyset]": 1084 1095 sol = [] 1085 1096 if repr(sol) == "[universalset]": 1086 from sage.rings.infinity import Infinity 1097 from sage.rings.infinity import Infinity 1087 1098 sol = [[i<Infinity for i in vars]] 1088 1099 return sol 1089 1100 1090 1101 def solve_ineq(ineq, vars=None): 1091 1102 """ 1092 Solves inequalities and systems of inequalities using Maxima. 1093 Switches between rational inequalities 1094 (sage.symbolic.relation.solve_ineq_rational) 1095 and fourier elimination (sage.symbolic.relation.solve_ineq_fouried). 1096 See the documentation of these functions for more details. 1097 1103 Solves inequalities and systems of inequalities using Maxima. 1104 Switches between rational inequalities 1105 (sage.symbolic.relation.solve_ineq_rational) 1106 and fourier elimination (sage.symbolic.relation.solve_ineq_fouried). 1107 See the documentation of these functions for more details. 1108 1098 1109 INPUT: 1099 1110 1100 1111  ``ineq``  one inequality or a list of inequalities … … 1110 1121 for some types of nonlinear inequalities. See 1111 1122 http://maxima.cvs.sourceforge.net/viewvc/maxima/maxima/share/contrib/fourier_elim/rtest_fourier_elim.mac 1112 1123 for a big gallery of problems covered by this algorithm. 1113 1124 1114 1125  ``vars``  optional parameter with list of variables. This list 1115 1126 is used only if fourier elimination is used. If omitted or if 1116 1127 rational inequality is solved, then variables are determined 1117 1128 automatically. 1118 1129 1119 OUTPUT: 1130 OUTPUT: 1120 1131 1121 1132  ``list``  output is list of solutions as a list of simple inequalities 1122 output [A,B,C] means (A or B or C) each A, B, C is again a list and 1133 output [A,B,C] means (A or B or C) each A, B, C is again a list and 1123 1134 if A=[a,b], then A means (a and b). 1124 1135 1125 1136 EXAMPLES:: 1126 1137 1127 1138 sage: from sage.symbolic.relation import solve_ineq 1128 1139 1129 1140 Inequalities in one variable. The variable is detected automatically:: 1130 1141 1131 sage: solve_ineq(x^21>3) 1142 sage: solve_ineq(x^21>3) 1132 1143 [[x < 2], [x > 2]] 1133 1144 1134 1145 sage: solve_ineq(1/(x1)<=8) 1135 1146 [[x < 1], [x >= (9/8)]] 1136 1147 1137 1148 System of inequalities with automatically detected inequalities:: 1138 1149 1139 1150 sage: y=var('y') … … 1150 1161 [[x < y, y < x + 3, x < (3/2)]] 1151 1162 1152 1163 ALGORITHM: 1153 1164 1154 1165 Calls solve_ineq_fourier if inequalities are list and 1155 1166 solve_ineq_univar of the inequality is symbolic expression. See 1156 1167 the description of these commands for more details related to the … … 1158 1169 there is no solution. 1159 1170 1160 1171 AUTHORS: 1161 1172 1162 1173  Robert Marik (012010) 1163 1174 """ 1164 1175 if isinstance(ineq,list):