Ticket #14801: trac_14801_piecewiese.patch
File trac_14801_piecewiese.patch, 141.0 KB (added by , 6 years ago) 


sage/functions/all.py
# HG changeset patch # User Volker Braun <vbraun@stp.dias.ie> # Date 1372346122 14400 # Thu Jun 27 11:15:22 2013 0400 # Node ID 530f6f96f708d3565c153752a7c6d04ed07c8bd0 # Parent 3f607223351f015a7fef512ba512497ed782e0b9 Piecewise functions using the symbolic expression framework diff git a/sage/functions/all.py b/sage/functions/all.py
a b 1 from piecewise import piecewise, Piecewise 1 from sage.misc.lazy_import import lazy_import 2 3 lazy_import('sage.functions.piecewise', 'Piecewise') # deprecated 4 lazy_import('sage.functions.piecewise', 'piecewise') 2 5 3 6 from trig import ( sin, cos, sec, csc, cot, tan, 4 7 asin, acos, atan, 
new file sage/functions/piecewiseold.py
diff git a/sage/functions/piecewiseold.py b/sage/functions/piecewiseold.py new file mode 100644
 + 1 r""" 2 Piecewisedefined Functions 3 4 Sage implements a very simple class of piecewisedefined functions. 5 Functions may be any type of symbolic expression. Infinite 6 intervals are not supported. The endpoints of each interval must 7 line up. 8 9 TODO: 10 11  Implement max/min location and values, 12 13  Need: parent object  ring of piecewise functions 14 15  This class should derive from an elementtype class, and should 16 define ``_add_``, ``_mul_``, etc. That will automatically take care 17 of left multiplication and proper coercion. The coercion mentioned 18 below for scalar mult on right is bad, since it only allows ints and 19 rationals. The right way is to use an element class and only define 20 ``_mul_``, and have a parent, so anything gets coerced properly. 21 22 AUTHORS: 23 24  David Joyner (200604): initial version 25 26  David Joyner (200609): added __eq__, extend_by_zero_to, unextend, 27 convolution, trapezoid, trapezoid_integral_approximation, 28 riemann_sum, riemann_sum_integral_approximation, tangent_line fixed 29 bugs in __mul__, __add__ 30 31  David Joyner (200703): adding Hann filter for FS, added general FS 32 filter methods for computing and plotting, added options to plotting 33 of FS (eg, specifying rgb values are now allowed). Fixed bug in 34 documentation reported by Pablo De Napoli. 35 36  David Joyner (200709): bug fixes due to behaviour of 37 SymbolicArithmetic 38 39  David Joyner (200804): fixed docstring bugs reported by J Morrow; added 40 support for Laplace transform of functions with infinite support. 41 42  David Joyner (200807): fixed a left multiplication bug reported by 43 C. Boncelet (by defining __rmul__ = __mul__). 44 45  Paul Butler (200901): added indefinite integration and default_variable 46 47 TESTS:: 48 49 sage: R.<x> = QQ[] 50 sage: f = Piecewise([[(0,1),1*x^0]]) 51 sage: 2*f 52 Piecewise defined function with 1 parts, [[(0, 1), 2]] 53 """ 54 55 #***************************************************************************** 56 # Copyright (C) 2006 William Stein <wstein@gmail.com> 57 # 2006 David Joyner <wdjoyner@gmail.com> 58 # 59 # Distributed under the terms of the GNU General Public License (GPL) 60 # 61 # This code is distributed in the hope that it will be useful, 62 # but WITHOUT ANY WARRANTY; without even the implied warranty of 63 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 64 # General Public License for more details. 65 # 66 # The full text of the GPL is available at: 67 # 68 # http://www.gnu.org/licenses/ 69 #***************************************************************************** 70 71 from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing 72 from sage.misc.sage_eval import sage_eval 73 from sage.rings.all import QQ, RR, Integer, Rational, infinity 74 from sage.calculus.functional import derivative 75 from sage.symbolic.expression import is_Expression 76 from sage.symbolic.assumptions import assume, forget 77 78 from sage.calculus.calculus import SR, maxima 79 from sage.calculus.all import var 80 81 def piecewise(list_of_pairs, var=None): 82 """ 83 Returns a piecewise function from a list of (interval, function) 84 pairs. 85 86 ``list_of_pairs`` is a list of pairs (I, fcn), where 87 fcn is a Sage function (such as a polynomial over RR, or functions 88 using the lambda notation), and I is an interval such as I = (1,3). 89 Two consecutive intervals must share a common endpoint. 90 91 If the optional ``var`` is specified, then any symbolic expressions 92 in the list will be converted to symbolic functions using 93 ``fcn.function(var)``. (This says which variable is considered to 94 be "piecewise".) 95 96 We assume that these definitions are consistent (ie, no checking is 97 done). 98 99 EXAMPLES:: 100 101 sage: f1(x) = 1 102 sage: f2(x) = 2 103 sage: f = Piecewise([[(0,pi/2),f1],[(pi/2,pi),f2]]) 104 sage: f(1) 105 1 106 sage: f(3) 107 2 108 sage: f = Piecewise([[(0,1),x], [(1,2),x^2]], x); f 109 Piecewise defined function with 2 parts, [[(0, 1), x > x], [(1, 2), x > x^2]] 110 sage: f(0.9) 111 0.900000000000000 112 sage: f(1.1) 113 1.21000000000000 114 """ 115 return PiecewisePolynomial(list_of_pairs, var=var) 116 117 Piecewise = piecewise 118 119 class PiecewisePolynomial: 120 """ 121 Returns a piecewise function from a list of (interval, function) 122 pairs. 123 124 EXAMPLES:: 125 126 sage: f1(x) = 1 127 sage: f2(x) = 2 128 sage: f = Piecewise([[(0,pi/2),f1],[(pi/2,pi),f2]]) 129 sage: f(1) 130 1 131 sage: f(3) 132 2 133 """ 134 def __init__(self, list_of_pairs, var=None): 135 r""" 136 ``list_of_pairs`` is a list of pairs (I, fcn), where 137 fcn is a Sage function (such as a polynomial over RR, or functions 138 using the lambda notation), and I is an interval such as I = (1,3). 139 Two consecutive intervals must share a common endpoint. 140 141 If the optional ``var`` is specified, then any symbolic 142 expressions in the list will be converted to symbolic 143 functions using ``fcn.function(var)``. (This says which 144 variable is considered to be "piecewise".) 145 146 We assume that these definitions are consistent (ie, no checking is 147 done). 148 149 EXAMPLES:: 150 151 sage: f1(x) = 1 152 sage: f2(x) = 1  x 153 sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]) 154 sage: f.list() 155 [[(0, 1), x > 1], [(1, 2), x > x + 1]] 156 sage: f.length() 157 2 158 """ 159 self._length = len(list_of_pairs) 160 self._intervals = [x[0] for x in list_of_pairs] 161 functions = [x[1] for x in list_of_pairs] 162 if var is not None: 163 for i in range(len(functions)): 164 if is_Expression(functions[i]): 165 functions[i] = functions[i].function(var) 166 self._functions = functions 167 # We regenerate self._list in case self._functions was modified 168 # above. This also protects us in case somebody mutates a list 169 # after they use it as an argument to piecewise(). 170 self._list = [[self._intervals[i], self._functions[i]] for i in range(self._length)] 171 172 def list(self): 173 """ 174 Returns the pieces of this function as a list of functions. 175 176 EXAMPLES:: 177 178 sage: f1(x) = 1 179 sage: f2(x) = 1  x 180 sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]) 181 sage: f.list() 182 [[(0, 1), x > 1], [(1, 2), x > x + 1]] 183 """ 184 return self._list 185 186 def length(self): 187 """ 188 Returns the number of pieces of this function. 189 190 EXAMPLES:: 191 192 sage: f1(x) = 1 193 sage: f2(x) = 1  x 194 sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]) 195 sage: f.length() 196 2 197 """ 198 return self._length 199 200 def __repr__(self): 201 """ 202 EXAMPLES:: 203 204 sage: f1(x) = 1 205 sage: f2(x) = 1  x 206 sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]); f 207 Piecewise defined function with 2 parts, [[(0, 1), x > 1], [(1, 2), x > x + 1]] 208 """ 209 return 'Piecewise defined function with %s parts, %s'%( 210 self.length(),self.list()) 211 212 def _latex_(self): 213 r""" 214 EXAMPLES:: 215 216 sage: f1(x) = 1 217 sage: f2(x) = 1  x 218 sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]) 219 sage: latex(f) 220 \begin{cases} 221 x \ {\mapsto}\ 1 &\text{on $(0, 1)$}\cr 222 x \ {\mapsto}\ x + 1 &\text{on $(1, 2)$}\cr 223 \end{cases} 224 225 :: 226 227 sage: f(x) = sin(x*pi/2) 228 sage: g(x) = 1(x1)^2 229 sage: h(x) = x 230 sage: P = Piecewise([[(0,1), f], [(1,3),g], [(3,5), h]]) 231 sage: latex(P) 232 \begin{cases} 233 x \ {\mapsto}\ \sin\left(\frac{1}{2} \, \pi x\right) &\text{on $(0, 1)$}\cr 234 x \ {\mapsto}\ {\left(x  1\right)}^{2} + 1 &\text{on $(1, 3)$}\cr 235 x \ {\mapsto}\ x &\text{on $(3, 5)$}\cr 236 \end{cases} 237 """ 238 from sage.misc.latex import latex 239 tex = ['\\begin{cases}\n'] 240 for (left, right), f in self.list(): 241 tex.append('%s &\\text{on $(%s, %s)$}\\cr\n' % (latex(f), left, right)) 242 tex.append(r'\end{cases}') 243 return ''.join(tex) 244 245 def intervals(self): 246 """ 247 A piecewise nonpolynomial example. 248 249 EXAMPLES:: 250 251 sage: f1(x) = 1 252 sage: f2(x) = 1x 253 sage: f3(x) = exp(x) 254 sage: f4(x) = sin(2*x) 255 sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]]) 256 sage: f.intervals() 257 [(0, 1), (1, 2), (2, 3), (3, 10)] 258 """ 259 return self._intervals 260 261 def domain(self): 262 """ 263 Returns the domain of the function. 264 265 EXAMPLES:: 266 267 sage: f1(x) = 1 268 sage: f2(x) = 1x 269 sage: f3(x) = exp(x) 270 sage: f4(x) = sin(2*x) 271 sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]]) 272 sage: f.domain() 273 (0, 10) 274 """ 275 endpoints = sum(self.intervals(), ()) 276 return (min(endpoints), max(endpoints)) 277 278 def functions(self): 279 """ 280 Returns the list of functions (the "pieces"). 281 282 EXAMPLES:: 283 284 sage: f1(x) = 1 285 sage: f2(x) = 1x 286 sage: f3(x) = exp(x) 287 sage: f4(x) = sin(2*x) 288 sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]]) 289 sage: f.functions() 290 [x > 1, x > x + 1, x > e^x, x > sin(2*x)] 291 """ 292 return self._functions 293 294 def extend_by_zero_to(self,xmin=1000,xmax=1000): 295 """ 296 This function simply returns the piecewise defined function which 297 is extended by 0 so it is defined on all of (xmin,xmax). This is 298 needed to add two piecewise functions in a reasonable way. 299 300 EXAMPLES:: 301 302 sage: f1(x) = 1 303 sage: f2(x) = 1  x 304 sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]) 305 sage: f.extend_by_zero_to(1, 3) 306 Piecewise defined function with 4 parts, [[(1, 0), 0], [(0, 1), x > 1], [(1, 2), x > x + 1], [(2, 3), 0]] 307 """ 308 zero = QQ['x'](0) 309 list_of_pairs = self.list() 310 a, b = self.domain() 311 if xmin < a: 312 list_of_pairs = [[(xmin, a), zero]] + list_of_pairs 313 if xmax > b: 314 list_of_pairs = list_of_pairs + [[(b, xmax), zero]] 315 return Piecewise(list_of_pairs) 316 317 def unextend(self): 318 """ 319 This removes any parts in the front or back of the function which 320 is zero (the inverse to extend_by_zero_to). 321 322 EXAMPLES:: 323 324 sage: R.<x> = QQ[] 325 sage: f = Piecewise([[(3,1),1+2+x],[(1,1),1x^2]]) 326 sage: e = f.extend_by_zero_to(10,10); e 327 Piecewise defined function with 4 parts, [[(10, 3), 0], [(3, 1), x + 3], [(1, 1), x^2 + 1], [(1, 10), 0]] 328 sage: d = e.unextend(); d 329 Piecewise defined function with 2 parts, [[(3, 1), x + 3], [(1, 1), x^2 + 1]] 330 sage: d==f 331 True 332 """ 333 list_of_pairs = self.list() 334 funcs = self.functions() 335 if funcs[0] == 0: 336 list_of_pairs = list_of_pairs[1:] 337 if funcs[1] == 0: 338 list_of_pairs = list_of_pairs[:1] 339 return Piecewise(list_of_pairs) 340 341 def _riemann_sum_helper(self, N, func, initial=0): 342 """ 343 A helper function for computing Riemann sums. 344 345 INPUT: 346 347 348  ``N``  the number of subdivisions 349 350  ``func``  a function to apply to the endpoints of 351 each subdivision 352 353  ``initial``  the starting value 354 355 356 EXAMPLES:: 357 358 sage: f1(x) = x^2 ## example 1 359 sage: f2(x) = 5x^2 360 sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]) 361 sage: f._riemann_sum_helper(6, lambda x0, x1: (x1x0)*f(x1)) 362 19/6 363 """ 364 a,b = self.domain() 365 rsum = initial 366 h = (ba)/N 367 for i in range(N): 368 x0 = a+i*h 369 x1 = a+(i+1)*h 370 rsum += func(x0, x1) 371 return rsum 372 373 def riemann_sum_integral_approximation(self,N,mode=None): 374 """ 375 Returns the piecewise line function defined by the Riemann sums in 376 numerical integration based on a subdivision into N subintervals. 377 378 Set mode="midpoint" for the height of the rectangles to be 379 determined by the midpoint of the subinterval; set mode="right" for 380 the height of the rectangles to be determined by the righthand 381 endpoint of the subinterval; the default is mode="left" (the height 382 of the rectangles to be determined by the lefthand endpoint of 383 the subinterval). 384 385 EXAMPLES:: 386 387 sage: f1(x) = x^2 ## example 1 388 sage: f2(x) = 5x^2 389 sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]) 390 sage: f.riemann_sum_integral_approximation(6) 391 17/6 392 sage: f.riemann_sum_integral_approximation(6,mode="right") 393 19/6 394 sage: f.riemann_sum_integral_approximation(6,mode="midpoint") 395 3 396 sage: f.integral(definite=True) 397 3 398 """ 399 if mode is None: 400 return self._riemann_sum_helper(N, lambda x0, x1: (x1x0)*self(x0)) 401 elif mode == "right": 402 return self._riemann_sum_helper(N, lambda x0, x1: (x1x0)*self(x1)) 403 elif mode == "midpoint": 404 return self._riemann_sum_helper(N, lambda x0, x1: (x1x0)*self((x0+x1)/2)) 405 else: 406 raise ValueError, "invalid mode" 407 408 def riemann_sum(self,N,mode=None): 409 """ 410 Returns the piecewise line function defined by the Riemann sums in 411 numerical integration based on a subdivision into N subintervals. 412 Set mode="midpoint" for the height of the rectangles to be 413 determined by the midpoint of the subinterval; set mode="right" for 414 the height of the rectangles to be determined by the righthand 415 endpoint of the subinterval; the default is mode="left" (the height 416 of the rectangles to be determined by the lefthand endpoint of 417 the subinterval). 418 419 EXAMPLES:: 420 421 sage: f1(x) = x^2 422 sage: f2(x) = 5x^2 423 sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]) 424 sage: f.riemann_sum(6,mode="midpoint") 425 Piecewise defined function with 6 parts, [[(0, 1/3), 1/36], [(1/3, 2/3), 1/4], [(2/3, 1), 25/36], [(1, 4/3), 131/36], [(4/3, 5/3), 11/4], [(5/3, 2), 59/36]] 426 427 :: 428 429 sage: f = Piecewise([[(1,1),(1x^2).function(x)]]) 430 sage: rsf = f.riemann_sum(7) 431 sage: P = f.plot(rgbcolor=(0.7,0.1,0.5), plot_points=40) 432 sage: Q = rsf.plot(rgbcolor=(0.7,0.6,0.6), plot_points=40) 433 sage: L = add([line([[a,0],[a,f(x=a)]],rgbcolor=(0.7,0.6,0.6)) for (a,b),f in rsf.list()]) 434 sage: P + Q + L 435 436 :: 437 438 sage: f = Piecewise([[(1,1),(1/2+xx^3)]], x) ## example 3 439 sage: rsf = f.riemann_sum(8) 440 sage: P = f.plot(rgbcolor=(0.7,0.1,0.5), plot_points=40) 441 sage: Q = rsf.plot(rgbcolor=(0.7,0.6,0.6), plot_points=40) 442 sage: L = add([line([[a,0],[a,f(x=a)]],rgbcolor=(0.7,0.6,0.6)) for (a,b),f in rsf.list()]) 443 sage: P + Q + L 444 """ 445 if mode is None: 446 rsum = self._riemann_sum_helper(N, lambda x0,x1: [[(x0,x1),SR(self(x0))]], 447 initial=[]) 448 elif mode == "right": 449 rsum = self._riemann_sum_helper(N, lambda x0,x1: [[(x0,x1),SR(self(x1))]], 450 initial=[]) 451 elif mode == "midpoint": 452 rsum = self._riemann_sum_helper(N, lambda x0,x1: [[(x0,x1),SR(self((x0+x1)/2))]], 453 initial=[]) 454 else: 455 raise ValueError, "invalid mode" 456 return Piecewise(rsum) 457 458 def trapezoid(self,N): 459 """ 460 Returns the piecewise line function defined by the trapezoid rule 461 for numerical integration based on a subdivision into N 462 subintervals. 463 464 EXAMPLES:: 465 466 sage: R.<x> = QQ[] 467 sage: f1 = x^2 468 sage: f2 = 5x^2 469 sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]) 470 sage: f.trapezoid(4) 471 Piecewise defined function with 4 parts, [[(0, 1/2), 1/2*x], [(1/2, 1), 9/2*x  2], [(1, 3/2), 1/2*x + 2], [(3/2, 2), 7/2*x + 8]] 472 473 :: 474 475 sage: R.<x> = QQ[] 476 sage: f = Piecewise([[(1,1),1x^2]]) 477 sage: tf = f.trapezoid(4) 478 sage: P = f.plot(rgbcolor=(0.7,0.1,0.5), plot_points=40) 479 sage: Q = tf.plot(rgbcolor=(0.7,0.6,0.6), plot_points=40) 480 sage: L = add([line([[a,0],[a,f(a)]],rgbcolor=(0.7,0.6,0.6)) for (a,b),f in tf.list()]) 481 sage: P+Q+L 482 483 :: 484 485 sage: R.<x> = QQ[] 486 sage: f = Piecewise([[(1,1),1/2+xx^3]]) ## example 3 487 sage: tf = f.trapezoid(6) 488 sage: P = f.plot(rgbcolor=(0.7,0.1,0.5), plot_points=40) 489 sage: Q = tf.plot(rgbcolor=(0.7,0.6,0.6), plot_points=40) 490 sage: L = add([line([[a,0],[a,f(a)]],rgbcolor=(0.7,0.6,0.6)) for (a,b),f in tf.list()]) 491 sage: P+Q+L 492 493 TESTS: 494 495 Use variables other than x (:trac:`13836`):: 496 497 sage: R.<y> = QQ[] 498 sage: f1 = y^2 499 sage: f2 = 5y^2 500 sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]) 501 sage: f.trapezoid(4) 502 Piecewise defined function with 4 parts, [[(0, 1/2), 1/2*y], [(1/2, 1), 9/2*y  2], [(1, 3/2), 1/2*y + 2], [(3/2, 2), 7/2*y + 8]] 503 504 """ 505 x = QQ[self.default_variable()].gen() 506 def f(x0, x1): 507 f0, f1 = self(x0), self(x1) 508 return [[(x0,x1),f0+(f1f0)*(x1x0)**(1)*(xx0)]] 509 rsum = self._riemann_sum_helper(N, f, initial=[]) 510 return Piecewise(rsum) 511 512 def trapezoid_integral_approximation(self,N): 513 """ 514 Returns the approximation given by the trapezoid rule for numerical 515 integration based on a subdivision into N subintervals. 516 517 EXAMPLES:: 518 519 sage: f1(x) = x^2 ## example 1 520 sage: f2(x) = 1(1x)^2 521 sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]) 522 sage: P = f.plot(rgbcolor=(0.7,0.1,0.5), plot_points=40) 523 sage: tf = f.trapezoid(6) 524 sage: Q = tf.plot(rgbcolor=(0.7,0.6,0.6), plot_points=40) 525 sage: ta = f.trapezoid_integral_approximation(6) 526 sage: t = text('trapezoid approximation = %s'%ta, (1.5, 0.25)) 527 sage: a = f.integral(definite=True) 528 sage: tt = text('area under curve = %s'%a, (1.5, 0.5)) 529 sage: P + Q + t + tt 530 531 :: 532 533 sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]) ## example 2 534 sage: tf = f.trapezoid(4) 535 sage: ta = f.trapezoid_integral_approximation(4) 536 sage: Q = tf.plot(rgbcolor=(0.7,0.6,0.6), plot_points=40) 537 sage: t = text('trapezoid approximation = %s'%ta, (1.5, 0.25)) 538 sage: a = f.integral(definite=True) 539 sage: tt = text('area under curve = %s'%a, (1.5, 0.5)) 540 sage: P+Q+t+tt 541 """ 542 def f(x0, x1): 543 f0, f1 = self(x0), self(x1) 544 return ((f1+f0)/2)*(x1x0) 545 return self._riemann_sum_helper(N, f) 546 547 def critical_points(self): 548 """ 549 Return the critical points of this piecewise function. 550 551 .. warning:: 552 553 Uses maxima, which prints the warning to use results with 554 caution. Only works for piecewise functions whose parts are 555 polynomials with real critical not occurring on the 556 interval endpoints. 557 558 EXAMPLES:: 559 560 sage: R.<x> = QQ[] 561 sage: f1 = x^0 562 sage: f2 = 10*x  x^2 563 sage: f3 = 3*x^4  156*x^3 + 3036*x^2  26208*x 564 sage: f = Piecewise([[(0,3),f1],[(3,10),f2],[(10,20),f3]]) 565 sage: expected = [5, 12, 13, 14] 566 sage: all(abs(ea) < 0.001 for e,a in zip(expected, f.critical_points())) 567 True 568 569 TESTS: 570 571 Use variables other than x (:trac:`13836`):: 572 573 sage: R.<y> = QQ[] 574 sage: f1 = y^0 575 sage: f2 = 10*y  y^2 576 sage: f3 = 3*y^4  156*y^3 + 3036*y^2  26208*y 577 sage: f = Piecewise([[(0,3),f1],[(3,10),f2],[(10,20),f3]]) 578 sage: expected = [5, 12, 13, 14] 579 sage: all(abs(ea) < 0.001 for e,a in zip(expected, f.critical_points())) 580 True 581 """ 582 from sage.calculus.calculus import maxima 583 x = QQ[self.default_variable()].gen() 584 crit_pts = [] 585 for (a,b), f in self.list(): 586 for root in maxima.allroots(SR(f).diff(x)==0): 587 root = float(root.rhs()) 588 if a < root < b: 589 crit_pts.append(root) 590 return crit_pts 591 592 def base_ring(self): 593 """ 594 Returns the base ring of the function pieces. This 595 is useful when this class is extended. 596 597 EXAMPLES:: 598 599 sage: f1(x) = 1 600 sage: f2(x) = 1x 601 sage: f3(x) = x^25 602 sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3]]) 603 sage: base_ring(f) 604 Symbolic Ring 605 606 :: 607 608 sage: R.<x> = QQ[] 609 sage: f1 = x^0 610 sage: f2 = 10*x  x^2 611 sage: f3 = 3*x^4  156*x^3 + 3036*x^2  26208*x 612 sage: f = Piecewise([[(0,3),f1],[(3,10),f2],[(10,20),f3]]) 613 sage: f.base_ring() 614 Rational Field 615 """ 616 return (self.functions()[0]).base_ring() 617 618 def end_points(self): 619 """ 620 Returns a list of all interval endpoints for this function. 621 622 EXAMPLES:: 623 624 sage: f1(x) = 1 625 sage: f2(x) = 1x 626 sage: f3(x) = x^25 627 sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3]]) 628 sage: f.end_points() 629 [0, 1, 2, 3] 630 """ 631 intervals = self.intervals() 632 return [ intervals[0][0] ] + [b for a,b in intervals] 633 634 def __call__(self,x0): 635 """ 636 Evaluates self at x0. Returns the average value of the jump if x0 637 is an interior endpoint of one of the intervals of self and the 638 usual value otherwise. 639 640 EXAMPLES:: 641 642 sage: f1(x) = 1 643 sage: f2(x) = 1x 644 sage: f3(x) = exp(x) 645 sage: f4(x) = sin(2*x) 646 sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]]) 647 sage: f(0.5) 648 1 649 sage: f(5/2) 650 e^(5/2) 651 sage: f(5/2).n() 652 12.1824939607035 653 sage: f(1) 654 1/2 655 """ 656 #x0 = QQ(x0) ## does not allow for evaluation at pi 657 n = self.length() 658 endpts = self.end_points() 659 for i in range(1,n): 660 if x0 == endpts[i]: 661 return (self.functions()[i1](x0) + self.functions()[i](x0))/2 662 if x0 == endpts[0]: 663 return self.functions()[0](x0) 664 if x0 == endpts[n]: 665 return self.functions()[n1](x0) 666 for i in range(n): 667 if endpts[i] < x0 < endpts[i+1]: 668 return self.functions()[i](x0) 669 raise ValueError,"Value not defined outside of domain." 670 671 def which_function(self,x0): 672 """ 673 Returns the function piece used to evaluate self at x0. 674 675 EXAMPLES:: 676 677 sage: f1(z) = z 678 sage: f2(x) = 1x 679 sage: f3(y) = exp(y) 680 sage: f4(t) = sin(2*t) 681 sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]]) 682 sage: f.which_function(3/2) 683 x > x + 1 684 """ 685 for (a,b), f in self.list(): 686 if a <= x0 <= b: 687 return f 688 raise ValueError,"Function not defined outside of domain." 689 690 def default_variable(self): 691 r""" 692 Return the default variable. The default variable is defined as the 693 first variable in the first piece that has a variable. If no pieces have 694 a variable (each piece is a constant value), `x` is returned. 695 696 The result is cached. 697 698 AUTHOR: Paul Butler 699 700 EXAMPLES:: 701 702 sage: f1(x) = 1 703 sage: f2(x) = 5*x 704 sage: p = Piecewise([[(0,1),f1],[(1,4),f2]]) 705 sage: p.default_variable() 706 x 707 708 sage: f1 = 3*var('y') 709 sage: p = Piecewise([[(0,1),4],[(1,4),f1]]) 710 sage: p.default_variable() 711 y 712 713 """ 714 try: 715 return self.__default_variable 716 except AttributeError: 717 pass 718 for _, fun in self._list: 719 try: 720 fun = SR(fun) 721 if fun.variables(): 722 v = fun.variables()[0] 723 self.__default_variable = v 724 return v 725 except TypeError: 726 # pass if fun is lambda function 727 pass 728 # default to x 729 v = var('x') 730 self.__default_value = v 731 return v 732 733 def integral(self, x=None, a=None, b=None, definite=False): 734 r""" 735 By default, returns the indefinite integral of the function. 736 If definite=True is given, returns the definite integral. 737 738 AUTHOR: 739 740  Paul Butler 741 742 EXAMPLES:: 743 744 sage: f1(x) = 1x 745 sage: f = Piecewise([[(0,1),1],[(1,2),f1]]) 746 sage: f.integral(definite=True) 747 1/2 748 749 :: 750 751 sage: f1(x) = 1 752 sage: f2(x) = 2 753 sage: f = Piecewise([[(0,pi/2),f1],[(pi/2,pi),f2]]) 754 sage: f.integral(definite=True) 755 1/2*pi 756 757 sage: f1(x) = 2 758 sage: f2(x) = 3  x 759 sage: f = Piecewise([[(2, 0), f1], [(0, 3), f2]]) 760 sage: f.integral() 761 Piecewise defined function with 2 parts, [[(2, 0), x > 2*x + 4], [(0, 3), x > 1/2*x^2 + 3*x + 4]] 762 763 sage: f1(y) = 1 764 sage: f2(y) = y + 3 765 sage: f3(y) = y  1 766 sage: f4(y) = y^2  1 767 sage: f5(y) = 3 768 sage: f = Piecewise([[(4,3),f1],[(3,2),f2],[(2,0),f3],[(0,2),f4],[(2,3),f5]]) 769 sage: F = f.integral(y) 770 sage: F 771 Piecewise defined function with 5 parts, [[(4, 3), y > y  4], [(3, 2), y > 1/2*y^2 + 3*y + 7/2], [(2, 0), y > 1/2*y^2  y  1/2], [(0, 2), y > 1/3*y^3  y  1/2], [(2, 3), y > 3*y  35/6]] 772 773 Ensure results are consistent with FTC:: 774 775 sage: F(3)  F(4) 776 1 777 sage: F(1)  F(3) 778 1 779 sage: F(2)  F(0) 780 2/3 781 sage: f.integral(y, 0, 2) 782 2/3 783 sage: F(3)  F(4) 784 19/6 785 sage: f.integral(y, 4, 3) 786 19/6 787 sage: f.integral(definite=True) 788 19/6 789 790 :: 791 792 sage: f1(y) = (y+3)^2 793 sage: f2(y) = y+3 794 sage: f3(y) = 3 795 sage: f = Piecewise([[(infinity, 3), f1], [(3, 0), f2], [(0, infinity), f3]]) 796 sage: f.integral() 797 Piecewise defined function with 3 parts, [[(Infinity, 3), y > 1/3*y^3 + 3*y^2 + 9*y + 9], [(3, 0), y > 1/2*y^2 + 3*y + 9/2], [(0, +Infinity), y > 3*y + 9/2]] 798 799 :: 800 801 sage: f1(x) = e^(abs(x)) 802 sage: f = Piecewise([[(infinity, infinity), f1]]) 803 sage: f.integral(definite=True) 804 2 805 sage: f.integral() 806 Piecewise defined function with 1 parts, [[(Infinity, +Infinity), x > 1/2*((sgn(x)  1)*e^(2*x)  2*e^x*sgn(x) + sgn(x) + 1)*e^(x)  1]] 807 808 :: 809 810 sage: f = Piecewise([((0, 5), cos(x))]) 811 sage: f.integral() 812 Piecewise defined function with 1 parts, [[(0, 5), x > sin(x)]] 813 814 815 TESTS: 816 817 Verify that piecewise integrals of zero work (trac #10841):: 818 819 sage: f0(x) = 0 820 sage: f = Piecewise([[(0,1),f0]]) 821 sage: f.integral(x,0,1) 822 0 823 sage: f = Piecewise([[(0,1), 0]]) 824 sage: f.integral(x,0,1) 825 0 826 sage: f = Piecewise([[(0,1), SR(0)]]) 827 sage: f.integral(x,0,1) 828 0 829 830 """ 831 if a != None and b != None: 832 F = self.integral(x) 833 return F(b)  F(a) 834 835 if a != None or b != None: 836 raise TypeError, 'only one endpoint given' 837 838 area = 0 # cumulative definite integral of parts to the left of the current interval 839 integrand_pieces = self.list() 840 integrand_pieces.sort() 841 new_pieces = [] 842 843 if x == None: 844 x = self.default_variable() 845 846 # The integral is computed by iterating over the pieces in order. 847 # The definite integral for each piece is calculated and accumulated in `area`. 848 # Thus at any time, `area` represents the definite integral of all the pieces 849 # encountered so far. The indefinite integral of each piece is also calculated, 850 # and the `area` before each piece is added to the piece. 851 # 852 # If a definite integral is requested, `area` is returned. 853 # Otherwise, a piecewise function is constructed from the indefinite integrals 854 # and returned. 855 # 856 # An exception is made if integral is called on a piecewise function 857 # that starts at infinity. In this case, we do not try to calculate the 858 # definite integral of the first piece, and the value of `area` remains 0 859 # after the first piece. 860 861 for (start, end), fun in integrand_pieces: 862 fun = SR(fun) 863 if start == infinity and not definite: 864 fun_integrated = fun.integral(x, end, x) 865 else: 866 try: 867 assume(start < x) 868 except ValueError: # Assumption is redundant 869 pass 870 fun_integrated = fun.integral(x, start, x) + area 871 forget(start < x) 872 if definite or end != infinity: 873 area += fun.integral(x, start, end) 874 new_pieces.append([(start, end), SR(fun_integrated).function(x)]) 875 876 if definite: 877 return SR(area) 878 else: 879 return Piecewise(new_pieces) 880 881 def convolution(self, other): 882 """ 883 Returns the convolution function, 884 `f*g(t)=\int_{\infty}^\infty f(u)g(tu)du`, for compactly 885 supported `f,g`. 886 887 EXAMPLES:: 888 889 sage: x = PolynomialRing(QQ,'x').gen() 890 sage: f = Piecewise([[(0,1),1*x^0]]) ## example 0 891 sage: g = f.convolution(f) 892 sage: h = f.convolution(g) 893 sage: P = f.plot(); Q = g.plot(rgbcolor=(1,1,0)); R = h.plot(rgbcolor=(0,1,1)); 894 sage: # Type show(P+Q+R) to view 895 sage: f = Piecewise([[(0,1),1*x^0],[(1,2),2*x^0],[(2,3),1*x^0]]) ## example 1 896 sage: g = f.convolution(f) 897 sage: h = f.convolution(g) 898 sage: P = f.plot(); Q = g.plot(rgbcolor=(1,1,0)); R = h.plot(rgbcolor=(0,1,1)); 899 sage: # Type show(P+Q+R) to view 900 sage: f = Piecewise([[(1,1),1]]) ## example 2 901 sage: g = Piecewise([[(0,3),x]]) 902 sage: f.convolution(g) 903 Piecewise defined function with 3 parts, [[(1, 1), 0], [(1, 2), 3/2*x], [(2, 4), 3/2*x]] 904 sage: g = Piecewise([[(0,3),1*x^0],[(3,4),2*x^0]]) 905 sage: f.convolution(g) 906 Piecewise defined function with 5 parts, [[(1, 1), x + 1], [(1, 2), 3], [(2, 3), x], [(3, 4), x + 8], [(4, 5), 2*x + 10]] 907 """ 908 f = self 909 g = other 910 M = min(min(f.end_points()),min(g.end_points())) 911 N = max(max(f.end_points()),max(g.end_points())) 912 R2 = PolynomialRing(QQ,2,names=["tt","uu"]) 913 tt,uu = R2.gens() 914 conv = 0 915 f0 = f.functions()[0] 916 g0 = g.functions()[0] 917 R1 = f0.parent() 918 xx = R1.gen() 919 var = repr(xx) 920 if len(f.intervals())==1 and len(g.intervals())==1: 921 f = f.unextend() 922 g = g.unextend() 923 a1 = f.intervals()[0][0] 924 a2 = f.intervals()[0][1] 925 b1 = g.intervals()[0][0] 926 b2 = g.intervals()[0][1] 927 i1 = repr(f0).replace(var,repr(uu)) 928 i2 = repr(g0).replace(var,"("+repr(ttuu)+")") 929 cmd1 = "integrate((%s)*(%s),%s,%s,%s)"%(i1,i2, uu, a1, ttb1) ## if a1+b1 < tt < a2+b1 930 cmd2 = "integrate((%s)*(%s),%s,%s,%s)"%(i1,i2, uu, ttb2, ttb1) ## if a1+b2 < tt < a2+b1 931 cmd3 = "integrate((%s)*(%s),%s,%s,%s)"%(i1,i2, uu, ttb2, a2) ## if a1+b2 < tt < a2+b2 932 cmd4 = "integrate((%s)*(%s),%s,%s,%s)"%(i1,i2, uu, a1, a2) ## if a2+b1 < tt < a1+b2 933 conv1 = maxima.eval(cmd1) 934 conv2 = maxima.eval(cmd2) 935 conv3 = maxima.eval(cmd3) 936 conv4 = maxima.eval(cmd4) 937 # this is a very, very, very ugly hack 938 x = PolynomialRing(QQ,'x').gen() 939 fg1 = sage_eval(conv1.replace("tt",var), {'x':x}) ## should be = R2(conv1) 940 fg2 = sage_eval(conv2.replace("tt",var), {'x':x}) ## should be = R2(conv2) 941 fg3 = sage_eval(conv3.replace("tt",var), {'x':x}) ## should be = R2(conv3) 942 fg4 = sage_eval(conv4.replace("tt",var), {'x':x}) ## should be = R2(conv4) 943 if a1b1<a2b2: 944 if a2+b1!=a1+b2: 945 h = Piecewise([[(a1+b1,a1+b2),fg1],[(a1+b2,a2+b1),fg4],[(a2+b1,a2+b2),fg3]]) 946 else: 947 h = Piecewise([[(a1+b1,a1+b2),fg1],[(a1+b2,a2+b2),fg3]]) 948 else: 949 if a1+b2!=a2+b1: 950 h = Piecewise([[(a1+b1,a2+b1),fg1],[(a2+b1,a1+b2),fg2],[(a1+b2,a2+b2),fg3]]) 951 else: 952 h = Piecewise([[(a1+b1,a2+b1),fg1],[(a2+b1,a2+b2),fg3]]) 953 return h 954 955 if len(f.intervals())>1 or len(g.intervals())>1: 956 z = Piecewise([[(3*abs(NM),3*abs(NM)),0*xx**0]]) 957 ff = f.functions() 958 gg = g.functions() 959 intvlsf = f.intervals() 960 intvlsg = g.intervals() 961 for i in range(len(ff)): 962 for j in range(len(gg)): 963 f0 = Piecewise([[intvlsf[i],ff[i]]]) 964 g0 = Piecewise([[intvlsg[j],gg[j]]]) 965 h = g0.convolution(f0) 966 z = z + h 967 return z.unextend() 968 969 def derivative(self): 970 r""" 971 Returns the derivative (as computed by maxima) 972 Piecewise(I,`(d/dx)(self_I)`), as I runs over the 973 intervals belonging to self. self must be piecewise polynomial. 974 975 EXAMPLES:: 976 977 sage: f1(x) = 1 978 sage: f2(x) = 1x 979 sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]) 980 sage: f.derivative() 981 Piecewise defined function with 2 parts, [[(0, 1), x > 0], [(1, 2), x > 1]] 982 sage: f1(x) = 1 983 sage: f2(x) = 2 984 sage: f = Piecewise([[(0,pi/2),f1],[(pi/2,pi),f2]]) 985 sage: f.derivative() 986 Piecewise defined function with 2 parts, [[(0, 1/2*pi), x > 0], [(1/2*pi, pi), x > 0]] 987 988 :: 989 990 sage: f = Piecewise([[(0,1), (x * 2)]], x) 991 sage: f.derivative() 992 Piecewise defined function with 1 parts, [[(0, 1), x > 2]] 993 """ 994 x = self.default_variable() 995 dlist = [[(a, b), derivative(f(x), x).function(x)] for (a,b),f in self.list()] 996 return Piecewise(dlist) 997 998 def tangent_line(self, pt): 999 """ 1000 Computes the linear function defining the tangent line of the 1001 piecewise function self. 1002 1003 EXAMPLES:: 1004 1005 sage: f1(x) = x^2 1006 sage: f2(x) = 5x^3+x 1007 sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]) 1008 sage: tf = f.tangent_line(0.9) ## tangent line at x=0.9 1009 sage: P = f.plot(rgbcolor=(0.7,0.1,0.5), plot_points=40) 1010 sage: Q = tf.plot(rgbcolor=(0.7,0.2,0.2), plot_points=40) 1011 sage: P + Q 1012 """ 1013 pt = QQ(pt) 1014 R = QQ[self.default_variable()] 1015 x = R.gen() 1016 der = self.derivative() 1017 tanline = (xpt)*der(pt)+self(pt) 1018 dlist = [[(a, b), tanline] for (a,b),f in self.list()] 1019 return Piecewise(dlist) 1020 1021 def plot(self, *args, **kwds): 1022 """ 1023 Returns the plot of self. 1024 1025 Keyword arguments are passed onto the plot command for each piece 1026 of the function. E.g., the plot_points keyword affects each 1027 segment of the plot. 1028 1029 EXAMPLES:: 1030 1031 sage: f1(x) = 1 1032 sage: f2(x) = 1x 1033 sage: f3(x) = exp(x) 1034 sage: f4(x) = sin(2*x) 1035 sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]]) 1036 sage: P = f.plot(rgbcolor=(0.7,0.1,0), plot_points=40) 1037 sage: P 1038 1039 Remember: to view this, type show(P) or P.save("path/myplot.png") 1040 and then open it in a graphics viewer such as GIMP. 1041 1042 TESTS: 1043 1044 We should not add each piece to the legend individually, since 1045 this creates duplicates (:trac:`12651`). This tests that only 1046 one of the graphics objects in the plot has a non``None`` 1047 ``legend_label``:: 1048 1049 sage: f1 = sin(x) 1050 sage: f2 = cos(x) 1051 sage: f = piecewise([[(1,0), f1],[(0,1), f2]]) 1052 sage: p = f.plot(legend_label='$f(x)$') 1053 sage: lines = [ 1054 ... line 1055 ... for line in p._objects 1056 ... if line.options()['legend_label'] is not None ] 1057 sage: len(lines) 1058 1 1059 """ 1060 from sage.plot.all import plot, Graphics 1061 1062 g = Graphics() 1063 1064 for i, ((a,b), f) in enumerate(self.list()): 1065 # If it's the first piece, pass all arguments. Otherwise, 1066 # filter out 'legend_label' so that we don't add each 1067 # piece to the legend separately (trac #12651). 1068 if i != 0 and 'legend_label' in kwds: 1069 del kwds['legend_label'] 1070 1071 g += plot(f, a, b, *args, **kwds) 1072 1073 return g 1074 1075 def fourier_series_cosine_coefficient(self,n,L): 1076 r""" 1077 Returns the nth Fourier series coefficient of 1078 `\cos(n\pi x/L)`, `a_n`. 1079 1080 INPUT: 1081 1082 1083  ``self``  the function f(x), defined over L x L 1084 1085  ``n``  an integer n=0 1086 1087  ``L``  (the period)/2 1088 1089 1090 OUTPUT: 1091 `a_n = \frac{1}{L}\int_{L}^L f(x)\cos(n\pi x/L)dx` 1092 1093 EXAMPLES:: 1094 1095 sage: f(x) = x^2 1096 sage: f = Piecewise([[(1,1),f]]) 1097 sage: f.fourier_series_cosine_coefficient(2,1) 1098 pi^(2) 1099 sage: f(x) = x^2 1100 sage: f = Piecewise([[(pi,pi),f]]) 1101 sage: f.fourier_series_cosine_coefficient(2,pi) 1102 1 1103 sage: f1(x) = 1 1104 sage: f2(x) = 2 1105 sage: f = Piecewise([[(pi,pi/2),f1],[(pi/2,pi),f2]]) 1106 sage: f.fourier_series_cosine_coefficient(5,pi) 1107 3/5/pi 1108 """ 1109 from sage.all import cos, pi 1110 x = var('x') 1111 result = sum([(f(x)*cos(pi*x*n/L)/L).integrate(x, a, b) 1112 for (a,b), f in self.list()]) 1113 if is_Expression(result): 1114 return result.simplify_trig() 1115 return result 1116 1117 def fourier_series_sine_coefficient(self,n,L): 1118 r""" 1119 Returns the nth Fourier series coefficient of 1120 `\sin(n\pi x/L)`, `b_n`. 1121 1122 INPUT: 1123 1124 1125  ``self``  the function f(x), defined over L x L 1126 1127  ``n``  an integer n0 1128 1129  ``L``  (the period)/2 1130 1131 1132 OUTPUT: 1133 `b_n = \frac{1}{L}\int_{L}^L f(x)\sin(n\pi x/L)dx` 1134 1135 EXAMPLES:: 1136 1137 sage: f(x) = x^2 1138 sage: f = Piecewise([[(1,1),f]]) 1139 sage: f.fourier_series_sine_coefficient(2,1) # L=1, n=2 1140 0 1141 """ 1142 from sage.all import sin, pi 1143 x = var('x') 1144 result = sum([(f(x)*sin(pi*x*n/L)/L).integrate(x, a, b) 1145 for (a,b), f in self.list()]) 1146 if is_Expression(result): 1147 return result.simplify_trig() 1148 return result 1149 1150 def _fourier_series_helper(self, N, L, scale_function): 1151 r""" 1152 A helper function for the construction of Fourier series. The 1153 argument scale_function is a function which takes in n, 1154 representing the `n^{th}` coefficient, and return an 1155 expression to scale the sine and cosine coefficients by. 1156 1157 EXAMPLES:: 1158 1159 sage: f(x) = x^2 1160 sage: f = Piecewise([[(1,1),f]]) 1161 sage: f._fourier_series_helper(3, 1, lambda n: 1) 1162 cos(2*pi*x)/pi^2  4*cos(pi*x)/pi^2 + 1/3 1163 """ 1164 from sage.all import pi, sin, cos, srange 1165 x = self.default_variable() 1166 a0 = self.fourier_series_cosine_coefficient(0,L) 1167 result = a0/2 + sum([(self.fourier_series_cosine_coefficient(n,L)*cos(n*pi*x/L) + 1168 self.fourier_series_sine_coefficient(n,L)*sin(n*pi*x/L))* 1169 scale_function(n) 1170 for n in srange(1,N)]) 1171 return result.expand() 1172 1173 1174 def fourier_series_partial_sum(self,N,L): 1175 r""" 1176 Returns the partial sum 1177 1178 .. math:: 1179 1180 f(x) \sim \frac{a_0}{2} + \sum_{n=1}^N [a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})], 1181 1182 as a string. 1183 1184 EXAMPLE:: 1185 1186 sage: f(x) = x^2 1187 sage: f = Piecewise([[(1,1),f]]) 1188 sage: f.fourier_series_partial_sum(3,1) 1189 cos(2*pi*x)/pi^2  4*cos(pi*x)/pi^2 + 1/3 1190 sage: f1(x) = 1 1191 sage: f2(x) = 2 1192 sage: f = Piecewise([[(pi,pi/2),f1],[(pi/2,pi),f2]]) 1193 sage: f.fourier_series_partial_sum(3,pi) 1194 3*cos(x)/pi  3*sin(2*x)/pi + 3*sin(x)/pi  1/4 1195 """ 1196 return self._fourier_series_helper(N, L, lambda n: 1) 1197 1198 def fourier_series_partial_sum_cesaro(self,N,L): 1199 r""" 1200 Returns the Cesaro partial sum 1201 1202 .. math:: 1203 1204 f(x) \sim \frac{a_0}{2} + \sum_{n=1}^N (1n/N)*[a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})], 1205 1206 1207 as a string. This is a "smoother" partial sum  the Gibbs 1208 phenomenon is mollified. 1209 1210 EXAMPLE:: 1211 1212 sage: f(x) = x^2 1213 sage: f = Piecewise([[(1,1),f]]) 1214 sage: f.fourier_series_partial_sum_cesaro(3,1) 1215 1/3*cos(2*pi*x)/pi^2  8/3*cos(pi*x)/pi^2 + 1/3 1216 sage: f1(x) = 1 1217 sage: f2(x) = 2 1218 sage: f = Piecewise([[(pi,pi/2),f1],[(pi/2,pi),f2]]) 1219 sage: f.fourier_series_partial_sum_cesaro(3,pi) 1220 2*cos(x)/pi  sin(2*x)/pi + 2*sin(x)/pi  1/4 1221 """ 1222 return self._fourier_series_helper(N, L, lambda n: 1n/N) 1223 1224 def fourier_series_partial_sum_hann(self,N,L): 1225 r""" 1226 Returns the Hannfiltered partial sum (named after von Hann, not 1227 Hamming) 1228 1229 .. math:: 1230 1231 f(x) \sim \frac{a_0}{2} + \sum_{n=1}^N H_N(n)*[a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})], 1232 1233 as a string, where `H_N(x) = (1+\cos(\pi x/N))/2`. This is 1234 a "smoother" partial sum  the Gibbs phenomenon is mollified. 1235 1236 EXAMPLE:: 1237 1238 sage: f(x) = x^2 1239 sage: f = Piecewise([[(1,1),f]]) 1240 sage: f.fourier_series_partial_sum_hann(3,1) 1241 1/4*cos(2*pi*x)/pi^2  3*cos(pi*x)/pi^2 + 1/3 1242 sage: f1(x) = 1 1243 sage: f2(x) = 2 1244 sage: f = Piecewise([[(pi,pi/2),f1],[(pi/2,pi),f2]]) 1245 sage: f.fourier_series_partial_sum_hann(3,pi) 1246 9/4*cos(x)/pi  3/4*sin(2*x)/pi + 9/4*sin(x)/pi  1/4 1247 """ 1248 from sage.all import cos, pi 1249 return self._fourier_series_helper(N, L, lambda n: (1+cos(pi*n/N))/2) 1250 1251 def fourier_series_partial_sum_filtered(self,N,L,F): 1252 r""" 1253 Returns the "filtered" partial sum 1254 1255 .. math:: 1256 1257 f(x) \sim \frac{a_0}{2} + \sum_{n=1}^N F_n*[a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})], 1258 1259 as a string, where `F = [F_1,F_2, ..., F_{N}]` is a list 1260 of length `N` consisting of real numbers. This can be used 1261 to plot FS solutions to the heat and wave PDEs. 1262 1263 EXAMPLE:: 1264 1265 sage: f(x) = x^2 1266 sage: f = Piecewise([[(1,1),f]]) 1267 sage: f.fourier_series_partial_sum_filtered(3,1,[1,1,1]) 1268 cos(2*pi*x)/pi^2  4*cos(pi*x)/pi^2 + 1/3 1269 sage: f1(x) = 1 1270 sage: f2(x) = 2 1271 sage: f = Piecewise([[(pi,pi/2),f1],[(pi/2,pi),f2]]) 1272 sage: f.fourier_series_partial_sum_filtered(3,pi,[1,1,1]) 1273 3*cos(x)/pi  3*sin(2*x)/pi + 3*sin(x)/pi  1/4 1274 """ 1275 return self._fourier_series_helper(N, L, lambda n: F[n]) 1276 1277 def plot_fourier_series_partial_sum(self,N,L,xmin,xmax, **kwds): 1278 r""" 1279 Plots the partial sum 1280 1281 .. math:: 1282 1283 f(x) \sim \frac{a_0}{2} + sum_{n=1}^N [a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})], 1284 1285 over xmin x xmin. 1286 1287 EXAMPLE:: 1288 1289 sage: f1(x) = 2 1290 sage: f2(x) = 1 1291 sage: f3(x) = 1 1292 sage: f4(x) = 2 1293 sage: f = Piecewise([[(pi,pi/2),f1],[(pi/2,0),f2],[(0,pi/2),f3],[(pi/2,pi),f4]]) 1294 sage: P = f.plot_fourier_series_partial_sum(3,pi,5,5) # long time 1295 sage: f1(x) = 1 1296 sage: f2(x) = 2 1297 sage: f = Piecewise([[(pi,pi/2),f1],[(pi/2,pi),f2]]) 1298 sage: P = f.plot_fourier_series_partial_sum(15,pi,5,5) # long time 1299 1300 Remember, to view this type show(P) or P.save("path/myplot.png") 1301 and then open it in a graphics viewer such as GIMP. 1302 """ 1303 from sage.plot.all import plot 1304 return plot(self.fourier_series_partial_sum(N,L), xmin, xmax, **kwds) 1305 1306 def plot_fourier_series_partial_sum_cesaro(self,N,L,xmin,xmax, **kwds): 1307 r""" 1308 Plots the partial sum 1309 1310 .. math:: 1311 1312 f(x) \sim \frac{a_0}{2} + \sum_{n=1}^N (1n/N)*[a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})], 1313 1314 1315 over xmin x xmin. This is a "smoother" partial sum  the Gibbs 1316 phenomenon is mollified. 1317 1318 EXAMPLE:: 1319 1320 sage: f1(x) = 2 1321 sage: f2(x) = 1 1322 sage: f3(x) = 1 1323 sage: f4(x) = 2 1324 sage: f = Piecewise([[(pi,pi/2),f1],[(pi/2,0),f2],[(0,pi/2),f3],[(pi/2,pi),f4]]) 1325 sage: P = f.plot_fourier_series_partial_sum_cesaro(3,pi,5,5) # long time 1326 sage: f1(x) = 1 1327 sage: f2(x) = 2 1328 sage: f = Piecewise([[(pi,pi/2),f1],[(pi/2,pi),f2]]) 1329 sage: P = f.plot_fourier_series_partial_sum_cesaro(15,pi,5,5) # long time 1330 1331 Remember, to view this type show(P) or P.save("path/myplot.png") 1332 and then open it in a graphics viewer such as GIMP. 1333 """ 1334 from sage.plot.all import plot 1335 return plot(self.fourier_series_partial_sum_cesaro(N,L), xmin, xmax, **kwds) 1336 1337 def plot_fourier_series_partial_sum_hann(self,N,L,xmin,xmax, **kwds): 1338 r""" 1339 Plots the partial sum 1340 1341 .. math:: 1342 1343 f(x) \sim \frac{a_0}{2} + \sum_{n=1}^N H_N(n)*[a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})], 1344 1345 1346 over xmin x xmin, where H_N(x) = (0.5)+(0.5)\*cos(x\*pi/N) is the 1347 Nth Hann filter. 1348 1349 EXAMPLE:: 1350 1351 sage: f1(x) = 2 1352 sage: f2(x) = 1 1353 sage: f3(x) = 1 1354 sage: f4(x) = 2 1355 sage: f = Piecewise([[(pi,pi/2),f1],[(pi/2,0),f2],[(0,pi/2),f3],[(pi/2,pi),f4]]) 1356 sage: P = f.plot_fourier_series_partial_sum_hann(3,pi,5,5) # long time 1357 sage: f1(x) = 1 1358 sage: f2(x) = 2 1359 sage: f = Piecewise([[(pi,pi/2),f1],[(pi/2,pi),f2]]) 1360 sage: P = f.plot_fourier_series_partial_sum_hann(15,pi,5,5) # long time 1361 1362 Remember, to view this type show(P) or P.save("path/myplot.png") 1363 and then open it in a graphics viewer such as GIMP. 1364 """ 1365 from sage.plot.all import plot 1366 return plot(self.fourier_series_partial_sum_hann(N,L), xmin, xmax, **kwds) 1367 1368 def plot_fourier_series_partial_sum_filtered(self,N,L,F,xmin,xmax, **kwds): 1369 r""" 1370 Plots the partial sum 1371 1372 .. math:: 1373 1374 f(x) \sim \frac{a_0}{2} + \sum_{n=1}^N F_n*[a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})], 1375 1376 1377 over xmin x xmin, where `F = [F_1,F_2, ..., F_{N}]` is a 1378 list of length `N` consisting of real numbers. This can be 1379 used to plot FS solutions to the heat and wave PDEs. 1380 1381 EXAMPLE:: 1382 1383 sage: f1(x) = 2 1384 sage: f2(x) = 1 1385 sage: f3(x) = 1 1386 sage: f4(x) = 2 1387 sage: f = Piecewise([[(pi,pi/2),f1],[(pi/2,0),f2],[(0,pi/2),f3],[(pi/2,pi),f4]]) 1388 sage: P = f.plot_fourier_series_partial_sum_filtered(3,pi,[1]*3,5,5) # long time 1389 sage: f1(x) = 1 1390 sage: f2(x) = 2 1391 sage: f = Piecewise([[(pi,pi/2),f1],[(pi/2,0),f2],[(0,pi/2),f1],[(pi/2,pi),f2]]) 1392 sage: P = f.plot_fourier_series_partial_sum_filtered(15,pi,[1]*15,5,5) # long time 1393 1394 Remember, to view this type show(P) or P.save("path/myplot.png") 1395 and then open it in a graphics viewer such as GIMP. 1396 """ 1397 from sage.plot.all import plot 1398 return plot(self.fourier_series_partial_sum_filtered(N,L,F), xmin, xmax, **kwds) 1399 1400 def fourier_series_value(self,x,L): 1401 r""" 1402 Returns the value of the Fourier series coefficient of self at 1403 `x`, 1404 1405 1406 .. math:: 1407 1408 f(x) \sim \frac{a_0}{2} + \sum_{n=1}^\infty [a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})], \ \ \ L<x<L. 1409 1410 1411 This method applies to piecewise nonpolynomial functions as well. 1412 1413 INPUT: 1414 1415 1416  ``self``  the function f(x), defined over L x L 1417 1418  ``x``  a real number 1419 1420  ``L``  (the period)/2 1421 1422 1423 OUTPUT: `(f^*(x+)+f^*(x)/2`, where `f^*` denotes 1424 the function `f` extended to `\RR` with period 1425 `2L` (Dirichlet's Theorem for Fourier series). 1426 1427 EXAMPLES:: 1428 1429 sage: f1(x) = 1 1430 sage: f2(x) = 1x 1431 sage: f3(x) = exp(x) 1432 sage: f4(x) = sin(2*x) 1433 sage: f = Piecewise([[(10,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]]) 1434 sage: f.fourier_series_value(101,10) 1435 1/2 1436 sage: f.fourier_series_value(100,10) 1437 1 1438 sage: f.fourier_series_value(10,10) 1439 1/2*sin(20) + 1/2 1440 sage: f.fourier_series_value(20,10) 1441 1 1442 sage: f.fourier_series_value(30,10) 1443 1/2*sin(20) + 1/2 1444 sage: f1(x) = 1 1445 sage: f2(x) = 2 1446 sage: f = Piecewise([[(pi,0),lambda x:0],[(0,pi/2),f1],[(pi/2,pi),f2]]) 1447 sage: f.fourier_series_value(1,pi) 1448 0 1449 sage: f.fourier_series_value(20,pi) 1450 1 1451 sage: f.fourier_series_value(pi/2,pi) 1452 1/2 1453 """ 1454 xnew = x  int(RR(x/(2*L)))*2*L 1455 endpts = self.end_points() 1456 if xnew == endpts[0] or xnew == endpts[1]: 1457 return (self.functions()[0](endpts[0]) + self.functions()[1](endpts[1]))/2 1458 else: 1459 return self(xnew) 1460 1461 def cosine_series_coefficient(self,n,L): 1462 r""" 1463 Returns the nth cosine series coefficient of 1464 `\cos(n\pi x/L)`, `a_n`. 1465 1466 INPUT: 1467 1468 1469  ``self``  the function f(x), defined over 0 x L (no 1470 checking is done to insure this) 1471 1472  ``n``  an integer n=0 1473 1474  ``L``  (the period)/2 1475 1476 1477 OUTPUT: 1478 `a_n = \frac{2}{L}\int_{L}^L f(x)\cos(n\pi x/L)dx` such 1479 that 1480 1481 .. math:: 1482 1483 f(x) \sim \frac{a_0}{2} + \sum_{n=1}^\infty a_n\cos(\frac{n\pi x}{L}),\ \ 0<x<L. 1484 1485 1486 1487 EXAMPLES:: 1488 1489 sage: f(x) = x 1490 sage: f = Piecewise([[(0,1),f]]) 1491 sage: f.cosine_series_coefficient(2,1) 1492 0 1493 sage: f.cosine_series_coefficient(3,1) 1494 4/9/pi^2 1495 sage: f1(x) = 1 1496 sage: f2(x) = 2 1497 sage: f = Piecewise([[(0,pi/2),f1],[(pi/2,pi),f2]]) 1498 sage: f.cosine_series_coefficient(2,pi) 1499 0 1500 sage: f.cosine_series_coefficient(3,pi) 1501 2/pi 1502 sage: f.cosine_series_coefficient(111,pi) 1503 2/37/pi 1504 sage: f1 = lambda x: x*(pix) 1505 sage: f = Piecewise([[(0,pi),f1]]) 1506 sage: f.cosine_series_coefficient(0,pi) 1507 1/3*pi^2 1508 1509 """ 1510 from sage.all import cos, pi 1511 x = var('x') 1512 result = sum([(2*f(x)*cos(pi*x*n/L)/L).integrate(x, a, b) 1513 for (a,b), f in self.list()]) 1514 if is_Expression(result): 1515 return result.simplify_trig() 1516 return result 1517 1518 1519 def sine_series_coefficient(self,n,L): 1520 r""" 1521 Returns the nth sine series coefficient of 1522 `\sin(n\pi x/L)`, `b_n`. 1523 1524 INPUT: 1525 1526  ``self``  the function f(x), defined over 0 x L (no 1527 checking is done to insure this) 1528 1529  ``n``  an integer n0 1530 1531  ``L``  (the period)/2 1532 1533 OUTPUT: 1534 1535 `b_n = \frac{2}{L}\int_{L}^L f(x)\sin(n\pi x/L)dx` such 1536 that 1537 1538 .. math:: 1539 1540 f(x) \sim \sum_{n=1}^\infty b_n\sin(\frac{n\pi x}{L}),\ \ 0<x<L. 1541 1542 EXAMPLES:: 1543 1544 sage: f(x) = 1 1545 sage: f = Piecewise([[(0,1),f]]) 1546 sage: f.sine_series_coefficient(2,1) 1547 0 1548 sage: f.sine_series_coefficient(3,1) 1549 4/3/pi 1550 """ 1551 from sage.all import sin, pi 1552 x = var('x') 1553 result = sum([(2*f(x)*sin(pi*x*n/L)/L).integrate(x, a, b) 1554 for (a,b), f in self.list()]) 1555 if is_Expression(result): 1556 return result.simplify_trig() 1557 return result 1558 1559 def laplace(self, x='x', s='t'): 1560 r""" 1561 Returns the Laplace transform of self with respect to the variable 1562 var. 1563 1564 INPUT: 1565 1566 1567  ``x``  variable of self 1568 1569  ``s``  variable of Laplace transform. 1570 1571 1572 We assume that a piecewise function is 0 outside of its domain and 1573 that the leftmost endpoint of the domain is 0. 1574 1575 EXAMPLES:: 1576 1577 sage: x, s, w = var('x, s, w') 1578 sage: f = Piecewise([[(0,1),1],[(1,2), 1x]]) 1579 sage: f.laplace(x, s) 1580 e^(s)/s + (s + 1)*e^(2*s)/s^2 + 1/s  e^(s)/s^2 1581 sage: f.laplace(x, w) 1582 e^(w)/w + (w + 1)*e^(2*w)/w^2 + 1/w  e^(w)/w^2 1583 1584 :: 1585 1586 sage: y, t = var('y, t') 1587 sage: f = Piecewise([[(1,2), 1y]]) 1588 sage: f.laplace(y, t) 1589 (t + 1)*e^(2*t)/t^2  e^(t)/t^2 1590 1591 :: 1592 1593 sage: s = var('s') 1594 sage: t = var('t') 1595 sage: f1(t) = t 1596 sage: f2(t) = 2 1597 sage: f = Piecewise([[(0,1),f1],[(1,infinity),f2]]) 1598 sage: f.laplace(t,s) 1599 (s + 1)*e^(s)/s^2 + 2*e^(s)/s  1/s^2 1600 """ 1601 from sage.all import assume, exp, forget 1602 x = var(x) 1603 s = var(s) 1604 assume(s>0) 1605 result = sum([(SR(f)*exp(s*x)).integral(x,a,b) 1606 for (a,b),f in self.list()]) 1607 forget(s>0) 1608 return result 1609 1610 def _make_compatible(self, other): 1611 """ 1612 Returns self and other extended to be defined on the same domain as 1613 well as a refinement of their intervals. This is used for adding 1614 and multiplying piecewise functions. 1615 1616 EXAMPLES:: 1617 1618 sage: R.<x> = QQ[] 1619 sage: f1 = Piecewise([[(0, 2), x]]) 1620 sage: f2 = Piecewise([[(1, 3), x^2]]) 1621 sage: f1._make_compatible(f2) 1622 (Piecewise defined function with 2 parts, [[(0, 2), x], [(2, 3), 0]], 1623 Piecewise defined function with 2 parts, [[(0, 1), 0], [(1, 3), x^2]], 1624 [(0, 1), (1, 2), (2, 3)]) 1625 """ 1626 a1, b1 = self.domain() 1627 a2, b2 = other.domain() 1628 a = min(a1, a2) 1629 b = max(b1, b2) 1630 F = self.extend_by_zero_to(a,b) 1631 G = other.extend_by_zero_to(a,b) 1632 endpts = list(set(F.end_points()).union(set(G.end_points()))) 1633 endpts.sort() 1634 return F, G, zip(endpts, endpts[1:]) 1635 1636 def __add__(self,other): 1637 """ 1638 Returns the piecewise defined function which is the sum of self and 1639 other. Does not require both domains be the same. 1640 1641 EXAMPLES:: 1642 1643 sage: x = PolynomialRing(QQ,'x').gen() 1644 sage: f1 = x^0 1645 sage: f2 = 1x 1646 sage: f3 = 2*x 1647 sage: f4 = 10x 1648 sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]]) 1649 sage: g1 = x2 1650 sage: g2 = x5 1651 sage: g = Piecewise([[(0,5),g1],[(5,10),g2]]) 1652 sage: h = f+g 1653 sage: h 1654 Piecewise defined function with 5 parts, [[(0, 1), x  1], [(1, 2), 1], [(2, 3), 3*x  2], [(3, 5), 8], [(5, 10), 5]] 1655 1656 Note that in this case the functions must be defined using 1657 polynomial expressions *not* using the lambda notation. 1658 """ 1659 F, G, intervals = self._make_compatible(other) 1660 fcn = [] 1661 for a,b in intervals: 1662 fcn.append([(a,b), F.which_function(b)+G.which_function(b)]) 1663 return Piecewise(fcn) 1664 1665 def __mul__(self,other): 1666 r""" 1667 Returns the piecewise defined function which is the product of one 1668 piecewise function (self) with another one (other). 1669 1670 EXAMPLES:: 1671 1672 sage: x = PolynomialRing(QQ,'x').gen() 1673 sage: f1 = x^0 1674 sage: f2 = 1x 1675 sage: f3 = 2*x 1676 sage: f4 = 10x 1677 sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]]) 1678 sage: g1 = x2 1679 sage: g2 = x5 1680 sage: g = Piecewise([[(0,5),g1],[(5,10),g2]]) 1681 sage: h = f*g 1682 sage: h 1683 Piecewise defined function with 5 parts, [[(0, 1), x  2], [(1, 2), x^2 + 3*x  2], [(2, 3), 2*x^2  4*x], [(3, 5), x^2 + 12*x  20], [(5, 10), x^2 + 15*x  50]] 1684 sage: g*(11/2) 1685 Piecewise defined function with 2 parts, [[(0, 5), 11/2*x  11], [(5, 10), 11/2*x  55/2]] 1686 1687 Note that in this method the functions must be defined using 1688 polynomial expressions *not* using the lambda notation. 1689 """ 1690 ## needed for scalar multiplication 1691 if isinstance(other,Rational) or isinstance(other,Integer): 1692 return Piecewise([[(a,b), other*f] for (a,b),f in self.list()]) 1693 else: 1694 F, G, intervals = self._make_compatible(other) 1695 fcn = [] 1696 for a,b in intervals: 1697 fcn.append([(a,b),F.which_function(b)*G.which_function(b)]) 1698 return Piecewise(fcn) 1699 1700 __rmul__ = __mul__ 1701 1702 def __eq__(self,other): 1703 """ 1704 Implements Boolean == operator. 1705 1706 EXAMPLES:: 1707 1708 sage: f1 = x^0 1709 sage: f2 = 1x 1710 sage: f3 = 2*x 1711 sage: f4 = 10x 1712 sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]]) 1713 sage: g = Piecewise([[(0,1),1],[(1,2),f2],[(2,3),f3],[(3,10),f4]]) 1714 sage: f==g 1715 True 1716 """ 1717 return self.list()==other.list() 
sage/functions/piecewise.py
diff git a/sage/functions/piecewise.py b/sage/functions/piecewise.py
a b 1 1 r""" 2 2 Piecewisedefined Functions 3 3 4 Sage implements a very simple class of piecewisedefined functions. 5 Functions may be any type of symbolic expression. Infinite 6 intervals are not supported. The endpoints of each interval must 7 line up. 4 This module implement piecewise functions in a single variable. See 5 :mod:`sage.sets.real_set` for more information about how to construct 6 subsets of the real line for the domains. 7 8 EXAMPLES:: 9 10 sage: f = piecewise([((0,1), x^3), ([1,0], x^2)]); f 11 piecewise(x>x^3 on (0, 1), x>x^2 on [1, 0]; x) 12 sage: 2*f 13 2*piecewise(x>x^3 on (0, 1), x>x^2 on [1, 0]; x) 14 sage: f(x=1/2) 15 1/8 16 sage: plot(f) # not tested 8 17 9 18 TODO: 10 19 11 20  Implement max/min location and values, 12 21 13  Need: parent object  ring of piecewise functions14 15  This class should derive from an elementtype class, and should16 define ``_add_``, ``_mul_``, etc. That will automatically take care17 of left multiplication and proper coercion. The coercion mentioned18 below for scalar mult on right is bad, since it only allows ints and19 rationals. The right way is to use an element class and only define20 ``_mul_``, and have a parent, so anything gets coerced properly.21 22 22 AUTHORS: 23 23 24 24  David Joyner (200604): initial version … … 44 44 45 45  Paul Butler (200901): added indefinite integration and default_variable 46 46 47  Volker Braun (2013): Complete rewrite 48 47 49 TESTS:: 48 50 49 sage: R.<x> = QQ[] 50 sage: f = Piecewise([[(0,1),1*x^0]]) 51 sage: 2*f 52 Piecewise defined function with 1 parts, [[(0, 1), 2]] 51 sage: fast_callable(f, vars=[x])(0.5) 52 0.125000000000... 53 53 """ 54 54 55 55 #***************************************************************************** 56 56 # Copyright (C) 2006 William Stein <wstein@gmail.com> 57 57 # 2006 David Joyner <wdjoyner@gmail.com> 58 # 2013 Volker Braun <vbraun.name@gmail.com> 58 59 # 59 # Distributed under the terms of the GNU General Public License (GPL) 60 # 61 # This code is distributed in the hope that it will be useful, 62 # but WITHOUT ANY WARRANTY; without even the implied warranty of 63 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 64 # General Public License for more details. 65 # 66 # The full text of the GPL is available at: 67 # 60 # Distributed under the terms of the GNU General Public License (GPL), 61 # version 2 or any later version. The full text of the GPL is available at: 68 62 # http://www.gnu.org/licenses/ 69 63 #***************************************************************************** 70 64 71 from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing 72 from sage.misc.sage_eval import sage_eval 73 from sage.rings.all import QQ, RR, Integer, Rational, infinity 74 from sage.calculus.functional import derivative 75 from sage.symbolic.expression import is_Expression 76 from sage.symbolic.assumptions import assume, forget 65 from sage.symbolic.function import BuiltinFunction 66 from sage.sets.real_set import RealSet 67 from sage.symbolic.ring import SR 77 68 78 from sage.calculus.calculus import SR, maxima79 from sage.calculus.all import var80 69 81 def piecewise(list_of_pairs, var=None):82 """83 Returns a piecewise function from a list of (interval, function)84 pairs.85 86 ``list_of_pairs`` is a list of pairs (I, fcn), where87 fcn is a Sage function (such as a polynomial over RR, or functions88 using the lambda notation), and I is an interval such as I = (1,3).89 Two consecutive intervals must share a common endpoint.90 91 If the optional ``var`` is specified, then any symbolic expressions92 in the list will be converted to symbolic functions using93 ``fcn.function(var)``. (This says which variable is considered to94 be "piecewise".)95 70 96 We assume that these definitions are consistent (ie, no checking is97 done).98 99 EXAMPLES::100 101 sage: f1(x) = 1102 sage: f2(x) = 2103 sage: f = Piecewise([[(0,pi/2),f1],[(pi/2,pi),f2]])104 sage: f(1)105 1106 sage: f(3)107 2108 sage: f = Piecewise([[(0,1),x], [(1,2),x^2]], x); f109 Piecewise defined function with 2 parts, [[(0, 1), x > x], [(1, 2), x > x^2]]110 sage: f(0.9)111 0.900000000000000112 sage: f(1.1)113 1.21000000000000114 """115 return PiecewisePolynomial(list_of_pairs, var=var)116 71 117 Piecewise = piecewise 72 class PiecewiseFunction(BuiltinFunction): 118 73 119 class PiecewisePolynomial: 120 """ 121 Returns a piecewise function from a list of (interval, function) 122 pairs. 123 124 EXAMPLES:: 125 126 sage: f1(x) = 1 127 sage: f2(x) = 2 128 sage: f = Piecewise([[(0,pi/2),f1],[(pi/2,pi),f2]]) 129 sage: f(1) 130 1 131 sage: f(3) 132 2 133 """ 134 def __init__(self, list_of_pairs, var=None): 135 r""" 136 ``list_of_pairs`` is a list of pairs (I, fcn), where 137 fcn is a Sage function (such as a polynomial over RR, or functions 138 using the lambda notation), and I is an interval such as I = (1,3). 139 Two consecutive intervals must share a common endpoint. 140 141 If the optional ``var`` is specified, then any symbolic 142 expressions in the list will be converted to symbolic 143 functions using ``fcn.function(var)``. (This says which 144 variable is considered to be "piecewise".) 145 146 We assume that these definitions are consistent (ie, no checking is 147 done). 74 def __init__(self): 75 """ 76 Piecewise function 148 77 149 78 EXAMPLES:: 150 79 151 sage: f1(x) = 1 152 sage: f2(x) = 1  x 153 sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]) 154 sage: f.list() 155 [[(0, 1), x > 1], [(1, 2), x > x + 1]] 156 sage: f.length() 157 2 80 sage: var('x, y') 81 (x, y) 82 sage: f = piecewise([((0,1), x^2*y), ([1,0], x*y^2)], var=x); f 83 piecewise(x>x^2*y on (0, 1), x>x*y^2 on [1, 0]; x) 84 sage: f(1/2) 85 1/4*y 86 sage: f(1/2) 87 1/2*y^2 158 88 """ 159 self._length = len(list_of_pairs)160 self._intervals = [x[0] for x in list_of_pairs]161 functions = [x[1] for x in list_of_pairs]162 if var is not None: 163 for i in range(len(functions)):164 if is_Expression(functions[i]):165 functions[i] = functions[i].function(var)166 self._functions = functions 167 # We regenerate self._list in case self._functions was modified168 # above. This also protects us in case somebody mutates a list169 # after they use it as an argument to piecewise().170 self._list = [[self._intervals[i], self._functions[i]] for i in range(self._length)]89 BuiltinFunction.__init__(self, "piecewise", 90 latex_name="piecewise", 91 conversions=dict(), nargs=2) 92 93 def __call__(self, function_pieces, **kwds): 94 r""" 95 Piecewise functions 96 97 INPUT: 98 99  ``function_pieces``  a list of pairs consisting of a 100 domain and a symbolic function. 171 101 172 def list(self): 102  ``var=x``  a symbolic variable or ``None`` (default). The 103 real variable in which the function is piecewise in. 104 105 OUTPUT: 106 107 A piecewisedefined function. A ``ValueError`` will be raised 108 if the domains of the pieces are not pairwise disjoint. 109 110 EXAMPLES:: 111 112 sage: my_abs = piecewise([((1, 0), x), ([0, 1], x)], var=x); my_abs 113 piecewise(x>x on (1, 0), x>x on [0, 1]; x) 114 sage: [ my_abs(i/5) for i in range(4, 5)] 115 [4/5, 3/5, 2/5, 1/5, 0, 1/5, 2/5, 3/5, 4/5] 116 117 TESTS:: 118 119 sage: piecewise([([1, 0], x), ([0, 1], x)], var=x) 120 Traceback (most recent call last): 121 ... 122 ValueError: domains must be pairwise disjoint 123 124 sage: step = piecewise([((1, 0), 1), ([0, 0], 0), ((0, 1), 1)], var=x); step 125 piecewise(x>1 on (1, 0), x>0 on {0}, x>1 on (0, 1); x) 126 sage: step(1/2), step(0), step(1/2) 127 (1, 0, 1) 173 128 """ 174 Returns the pieces of this function as a list of functions. 129 #print 'pf_call', function_pieces, kwds 130 var = kwds.pop('var', None) 131 parameters = [] 132 domain_list = [] 133 for piece in function_pieces: 134 domain, function = piece 135 if not isinstance(domain, RealSet): 136 domain = RealSet(domain) 137 if domain.is_empty(): 138 continue 139 function = SR(function) 140 if var is None and len(function.variables()) > 0: 141 var = function.variables()[0] 142 parameters.append((domain, function)) 143 domain_list.append(domain) 144 if not RealSet.are_pairwise_disjoint(*domain_list): 145 raise ValueError('domains must be pairwise disjoint') 146 if var is None: 147 var = self.default_variable() 148 parameters = SR._force_pyobject(tuple(parameters), recursive=False) 149 return BuiltinFunction.__call__(self, parameters, var, **kwds) 150 151 def _print_(self, parameters, variable): 152 """ 153 Return a string representation 154 155 OUTPUT: 156 157 String. 175 158 176 159 EXAMPLES:: 177 160 178 sage: f1(x) = 1 179 sage: f2(x) = 1  x 180 sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]) 181 sage: f.list() 182 [[(0, 1), x > 1], [(1, 2), x > x + 1]] 161 183 162 """ 184 return self._list 185 186 def length(self): 163 s = 'piecewise(' 164 args = [] 165 for domain, func in parameters: 166 args.append('{0}>{1} on {2}'.format(str(variable), str(func), str(domain))) 167 s += ', '.join(args) + '; {0})'.format(str(variable)) 168 return s 169 170 def _subs_(self, subs_map, options, parameters, x): 187 171 """ 188 Returns the number of pieces of this function.172 Callback from Pynac `subs()` 189 173 190 EXAMPLES:: 174 EXAMPLES: 175 176 If the substitution changes the piecewise variable, it must 177 evaluate to a number so that we know which component we are 178 on:: 179 180 sage: p = piecewise([((2, 0), x), ([0, 2], x)], var=x) 181 sage: p.subs(x=1) 182 1 183 sage: (10+p).subs(x=1) 184 11 185 186 Auxiliary variables can be substituted arbitrarily:: 187 188 sage: var('x,y') 189 (x, y) 190 sage: p = piecewise([((2, 0), x^y), ([0, 2], xy)], var=x); p 191 piecewise(x>x^y on (2, 0), x>x  y on [0, 2]; x) 192 sage: p.subs(y=sin(y)) 193 piecewise(x>x^sin(y) on (2, 0), x>x  sin(y) on [0, 2]; x) 194 """ 195 point = subs_map.apply_to(x, 0) 196 # print 'point =', point 197 if point == x: 198 # substitution only in auxiliary variables 199 new_params = [] 200 for domain, func in parameters: 201 new_params.append((domain, subs_map.apply_to(func, 0))) 202 return piecewise(new_params, var=x) 203 if not (point.is_numeric() and point.is_real()): 204 raise ValueError('substituting the piecewise variable must result in real number') 191 205 192 sage: f1(x) = 1 193 sage: f2(x) = 1  x 194 sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]) 195 sage: f.length() 196 2 206 for domain, func in parameters: 207 if domain.contains(point): 208 return subs_map.apply_to(func, 0) 209 raise ValueError('point is not in the domain') 210 211 @staticmethod 212 def in_operands(ex): 197 213 """ 198 return self._length 199 200 def __repr__(self): 201 """ 202 EXAMPLES:: 203 204 sage: f1(x) = 1 205 sage: f2(x) = 1  x 206 sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]); f 207 Piecewise defined function with 2 parts, [[(0, 1), x > 1], [(1, 2), x > x + 1]] 208 """ 209 return 'Piecewise defined function with %s parts, %s'%( 210 self.length(),self.list()) 211 212 def _latex_(self): 213 r""" 214 EXAMPLES:: 215 216 sage: f1(x) = 1 217 sage: f2(x) = 1  x 218 sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]) 219 sage: latex(f) 220 \begin{cases} 221 x \ {\mapsto}\ 1 &\text{on $(0, 1)$}\cr 222 x \ {\mapsto}\ x + 1 &\text{on $(1, 2)$}\cr 223 \end{cases} 224 225 :: 226 227 sage: f(x) = sin(x*pi/2) 228 sage: g(x) = 1(x1)^2 229 sage: h(x) = x 230 sage: P = Piecewise([[(0,1), f], [(1,3),g], [(3,5), h]]) 231 sage: latex(P) 232 \begin{cases} 233 x \ {\mapsto}\ \sin\left(\frac{1}{2} \, \pi x\right) &\text{on $(0, 1)$}\cr 234 x \ {\mapsto}\ {\left(x  1\right)}^{2} + 1 &\text{on $(1, 3)$}\cr 235 x \ {\mapsto}\ x &\text{on $(3, 5)$}\cr 236 \end{cases} 237 """ 238 from sage.misc.latex import latex 239 tex = ['\\begin{cases}\n'] 240 for (left, right), f in self.list(): 241 tex.append('%s &\\text{on $(%s, %s)$}\\cr\n' % (latex(f), left, right)) 242 tex.append(r'\end{cases}') 243 return ''.join(tex) 214 Return whether a symbolic expression contains a piecewise 215 function as operand 244 216 245 def intervals(self):246 """247 A piecewise nonpolynomial example.248 249 EXAMPLES::250 251 sage: f1(x) = 1252 sage: f2(x) = 1x253 sage: f3(x) = exp(x)254 sage: f4(x) = sin(2*x)255 sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]])256 sage: f.intervals()257 [(0, 1), (1, 2), (2, 3), (3, 10)]258 """259 return self._intervals260 261 def domain(self):262 """263 Returns the domain of the function.264 265 EXAMPLES::266 267 sage: f1(x) = 1268 sage: f2(x) = 1x269 sage: f3(x) = exp(x)270 sage: f4(x) = sin(2*x)271 sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]])272 sage: f.domain()273 (0, 10)274 """275 endpoints = sum(self.intervals(), ())276 return (min(endpoints), max(endpoints))277 278 def functions(self):279 """280 Returns the list of functions (the "pieces").281 282 EXAMPLES::283 284 sage: f1(x) = 1285 sage: f2(x) = 1x286 sage: f3(x) = exp(x)287 sage: f4(x) = sin(2*x)288 sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]])289 sage: f.functions()290 [x > 1, x > x + 1, x > e^x, x > sin(2*x)]291 """292 return self._functions293 294 def extend_by_zero_to(self,xmin=1000,xmax=1000):295 """296 This function simply returns the piecewise defined function which297 is extended by 0 so it is defined on all of (xmin,xmax). This is298 needed to add two piecewise functions in a reasonable way.299 300 EXAMPLES::301 302 sage: f1(x) = 1303 sage: f2(x) = 1  x304 sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])305 sage: f.extend_by_zero_to(1, 3)306 Piecewise defined function with 4 parts, [[(1, 0), 0], [(0, 1), x > 1], [(1, 2), x > x + 1], [(2, 3), 0]]307 """308 zero = QQ['x'](0)309 list_of_pairs = self.list()310 a, b = self.domain()311 if xmin < a:312 list_of_pairs = [[(xmin, a), zero]] + list_of_pairs313 if xmax > b:314 list_of_pairs = list_of_pairs + [[(b, xmax), zero]]315 return Piecewise(list_of_pairs)316 317 def unextend(self):318 """319 This removes any parts in the front or back of the function which320 is zero (the inverse to extend_by_zero_to).321 322 EXAMPLES::323 324 sage: R.<x> = QQ[]325 sage: f = Piecewise([[(3,1),1+2+x],[(1,1),1x^2]])326 sage: e = f.extend_by_zero_to(10,10); e327 Piecewise defined function with 4 parts, [[(10, 3), 0], [(3, 1), x + 3], [(1, 1), x^2 + 1], [(1, 10), 0]]328 sage: d = e.unextend(); d329 Piecewise defined function with 2 parts, [[(3, 1), x + 3], [(1, 1), x^2 + 1]]330 sage: d==f331 True332 """333 list_of_pairs = self.list()334 funcs = self.functions()335 if funcs[0] == 0:336 list_of_pairs = list_of_pairs[1:]337 if funcs[1] == 0:338 list_of_pairs = list_of_pairs[:1]339 return Piecewise(list_of_pairs)340 341 def _riemann_sum_helper(self, N, func, initial=0):342 """343 A helper function for computing Riemann sums.344 345 217 INPUT: 346 218 347 348  ``N``  the number of subdivisions 349 350  ``func``  a function to apply to the endpoints of 351 each subdivision 352 353  ``initial``  the starting value 354 355 356 EXAMPLES:: 357 358 sage: f1(x) = x^2 ## example 1 359 sage: f2(x) = 5x^2 360 sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]) 361 sage: f._riemann_sum_helper(6, lambda x0, x1: (x1x0)*f(x1)) 362 19/6 363 """ 364 a,b = self.domain() 365 rsum = initial 366 h = (ba)/N 367 for i in range(N): 368 x0 = a+i*h 369 x1 = a+(i+1)*h 370 rsum += func(x0, x1) 371 return rsum 219  ``ex``  a symbolic expression. 372 220 373 def riemann_sum_integral_approximation(self,N,mode=None): 374 """ 375 Returns the piecewise line function defined by the Riemann sums in 376 numerical integration based on a subdivision into N subintervals. 377 378 Set mode="midpoint" for the height of the rectangles to be 379 determined by the midpoint of the subinterval; set mode="right" for 380 the height of the rectangles to be determined by the righthand 381 endpoint of the subinterval; the default is mode="left" (the height 382 of the rectangles to be determined by the lefthand endpoint of 383 the subinterval). 384 385 EXAMPLES:: 386 387 sage: f1(x) = x^2 ## example 1 388 sage: f2(x) = 5x^2 389 sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]) 390 sage: f.riemann_sum_integral_approximation(6) 391 17/6 392 sage: f.riemann_sum_integral_approximation(6,mode="right") 393 19/6 394 sage: f.riemann_sum_integral_approximation(6,mode="midpoint") 395 3 396 sage: f.integral(definite=True) 397 3 398 """ 399 if mode is None: 400 return self._riemann_sum_helper(N, lambda x0, x1: (x1x0)*self(x0)) 401 elif mode == "right": 402 return self._riemann_sum_helper(N, lambda x0, x1: (x1x0)*self(x1)) 403 elif mode == "midpoint": 404 return self._riemann_sum_helper(N, lambda x0, x1: (x1x0)*self((x0+x1)/2)) 405 else: 406 raise ValueError, "invalid mode" 221 OUTPUT: 407 222 408 def riemann_sum(self,N,mode=None): 409 """ 410 Returns the piecewise line function defined by the Riemann sums in 411 numerical integration based on a subdivision into N subintervals. 412 Set mode="midpoint" for the height of the rectangles to be 413 determined by the midpoint of the subinterval; set mode="right" for 414 the height of the rectangles to be determined by the righthand 415 endpoint of the subinterval; the default is mode="left" (the height 416 of the rectangles to be determined by the lefthand endpoint of 417 the subinterval). 418 419 EXAMPLES:: 420 421 sage: f1(x) = x^2 422 sage: f2(x) = 5x^2 423 sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]) 424 sage: f.riemann_sum(6,mode="midpoint") 425 Piecewise defined function with 6 parts, [[(0, 1/3), 1/36], [(1/3, 2/3), 1/4], [(2/3, 1), 25/36], [(1, 4/3), 131/36], [(4/3, 5/3), 11/4], [(5/3, 2), 59/36]] 426 427 :: 428 429 sage: f = Piecewise([[(1,1),(1x^2).function(x)]]) 430 sage: rsf = f.riemann_sum(7) 431 sage: P = f.plot(rgbcolor=(0.7,0.1,0.5), plot_points=40) 432 sage: Q = rsf.plot(rgbcolor=(0.7,0.6,0.6), plot_points=40) 433 sage: L = add([line([[a,0],[a,f(x=a)]],rgbcolor=(0.7,0.6,0.6)) for (a,b),f in rsf.list()]) 434 sage: P + Q + L 435 436 :: 437 438 sage: f = Piecewise([[(1,1),(1/2+xx^3)]], x) ## example 3 439 sage: rsf = f.riemann_sum(8) 440 sage: P = f.plot(rgbcolor=(0.7,0.1,0.5), plot_points=40) 441 sage: Q = rsf.plot(rgbcolor=(0.7,0.6,0.6), plot_points=40) 442 sage: L = add([line([[a,0],[a,f(x=a)]],rgbcolor=(0.7,0.6,0.6)) for (a,b),f in rsf.list()]) 443 sage: P + Q + L 444 """ 445 if mode is None: 446 rsum = self._riemann_sum_helper(N, lambda x0,x1: [[(x0,x1),SR(self(x0))]], 447 initial=[]) 448 elif mode == "right": 449 rsum = self._riemann_sum_helper(N, lambda x0,x1: [[(x0,x1),SR(self(x1))]], 450 initial=[]) 451 elif mode == "midpoint": 452 rsum = self._riemann_sum_helper(N, lambda x0,x1: [[(x0,x1),SR(self((x0+x1)/2))]], 453 initial=[]) 454 else: 455 raise ValueError, "invalid mode" 456 return Piecewise(rsum) 457 458 def trapezoid(self,N): 459 """ 460 Returns the piecewise line function defined by the trapezoid rule 461 for numerical integration based on a subdivision into N 462 subintervals. 463 464 EXAMPLES:: 465 466 sage: R.<x> = QQ[] 467 sage: f1 = x^2 468 sage: f2 = 5x^2 469 sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]) 470 sage: f.trapezoid(4) 471 Piecewise defined function with 4 parts, [[(0, 1/2), 1/2*x], [(1/2, 1), 9/2*x  2], [(1, 3/2), 1/2*x + 2], [(3/2, 2), 7/2*x + 8]] 472 473 :: 474 475 sage: R.<x> = QQ[] 476 sage: f = Piecewise([[(1,1),1x^2]]) 477 sage: tf = f.trapezoid(4) 478 sage: P = f.plot(rgbcolor=(0.7,0.1,0.5), plot_points=40) 479 sage: Q = tf.plot(rgbcolor=(0.7,0.6,0.6), plot_points=40) 480 sage: L = add([line([[a,0],[a,f(a)]],rgbcolor=(0.7,0.6,0.6)) for (a,b),f in tf.list()]) 481 sage: P+Q+L 482 483 :: 484 485 sage: R.<x> = QQ[] 486 sage: f = Piecewise([[(1,1),1/2+xx^3]]) ## example 3 487 sage: tf = f.trapezoid(6) 488 sage: P = f.plot(rgbcolor=(0.7,0.1,0.5), plot_points=40) 489 sage: Q = tf.plot(rgbcolor=(0.7,0.6,0.6), plot_points=40) 490 sage: L = add([line([[a,0],[a,f(a)]],rgbcolor=(0.7,0.6,0.6)) for (a,b),f in tf.list()]) 491 sage: P+Q+L 492 493 TESTS: 494 495 Use variables other than x (:trac:`13836`):: 496 497 sage: R.<y> = QQ[] 498 sage: f1 = y^2 499 sage: f2 = 5y^2 500 sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]) 501 sage: f.trapezoid(4) 502 Piecewise defined function with 4 parts, [[(0, 1/2), 1/2*y], [(1/2, 1), 9/2*y  2], [(1, 3/2), 1/2*y + 2], [(3/2, 2), 7/2*y + 8]] 503 504 """ 505 x = QQ[self.default_variable()].gen() 506 def f(x0, x1): 507 f0, f1 = self(x0), self(x1) 508 return [[(x0,x1),f0+(f1f0)*(x1x0)**(1)*(xx0)]] 509 rsum = self._riemann_sum_helper(N, f, initial=[]) 510 return Piecewise(rsum) 511 512 def trapezoid_integral_approximation(self,N): 513 """ 514 Returns the approximation given by the trapezoid rule for numerical 515 integration based on a subdivision into N subintervals. 516 517 EXAMPLES:: 518 519 sage: f1(x) = x^2 ## example 1 520 sage: f2(x) = 1(1x)^2 521 sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]) 522 sage: P = f.plot(rgbcolor=(0.7,0.1,0.5), plot_points=40) 523 sage: tf = f.trapezoid(6) 524 sage: Q = tf.plot(rgbcolor=(0.7,0.6,0.6), plot_points=40) 525 sage: ta = f.trapezoid_integral_approximation(6) 526 sage: t = text('trapezoid approximation = %s'%ta, (1.5, 0.25)) 527 sage: a = f.integral(definite=True) 528 sage: tt = text('area under curve = %s'%a, (1.5, 0.5)) 529 sage: P + Q + t + tt 530 531 :: 532 533 sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]) ## example 2 534 sage: tf = f.trapezoid(4) 535 sage: ta = f.trapezoid_integral_approximation(4) 536 sage: Q = tf.plot(rgbcolor=(0.7,0.6,0.6), plot_points=40) 537 sage: t = text('trapezoid approximation = %s'%ta, (1.5, 0.25)) 538 sage: a = f.integral(definite=True) 539 sage: tt = text('area under curve = %s'%a, (1.5, 0.5)) 540 sage: P+Q+t+tt 541 """ 542 def f(x0, x1): 543 f0, f1 = self(x0), self(x1) 544 return ((f1+f0)/2)*(x1x0) 545 return self._riemann_sum_helper(N, f) 546 547 def critical_points(self): 548 """ 549 Return the critical points of this piecewise function. 550 551 .. warning:: 552 553 Uses maxima, which prints the warning to use results with 554 caution. Only works for piecewise functions whose parts are 555 polynomials with real critical not occurring on the 556 interval endpoints. 557 558 EXAMPLES:: 559 560 sage: R.<x> = QQ[] 561 sage: f1 = x^0 562 sage: f2 = 10*x  x^2 563 sage: f3 = 3*x^4  156*x^3 + 3036*x^2  26208*x 564 sage: f = Piecewise([[(0,3),f1],[(3,10),f2],[(10,20),f3]]) 565 sage: expected = [5, 12, 13, 14] 566 sage: all(abs(ea) < 0.001 for e,a in zip(expected, f.critical_points())) 567 True 568 569 TESTS: 570 571 Use variables other than x (:trac:`13836`):: 572 573 sage: R.<y> = QQ[] 574 sage: f1 = y^0 575 sage: f2 = 10*y  y^2 576 sage: f3 = 3*y^4  156*y^3 + 3036*y^2  26208*y 577 sage: f = Piecewise([[(0,3),f1],[(3,10),f2],[(10,20),f3]]) 578 sage: expected = [5, 12, 13, 14] 579 sage: all(abs(ea) < 0.001 for e,a in zip(expected, f.critical_points())) 580 True 581 """ 582 from sage.calculus.calculus import maxima 583 x = QQ[self.default_variable()].gen() 584 crit_pts = [] 585 for (a,b), f in self.list(): 586 for root in maxima.allroots(SR(f).diff(x)==0): 587 root = float(root.rhs()) 588 if a < root < b: 589 crit_pts.append(root) 590 return crit_pts 591 592 def base_ring(self): 593 """ 594 Returns the base ring of the function pieces. This 595 is useful when this class is extended. 223 Boolean 596 224 597 225 EXAMPLES:: 598 226 599 sage: f1(x) = 1 600 sage: f2(x) = 1x 601 sage: f3(x) = x^25 602 sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3]]) 603 sage: base_ring(f) 604 Symbolic Ring 227 sage: f = piecewise([([0,0], sin(x)), ((0,2), cos(x))]); f 228 piecewise(x>sin(x) on {0}, x>cos(x) on (0, 2); x) 229 sage: piecewise.in_operands(f) 230 True 231 sage: piecewise.in_operands(1+sin(f)) 232 True 233 sage: piecewise.in_operands(1+sin(0*f)) 234 False 235 """ 236 def is_piecewise(ex): 237 result = ex.operator() is piecewise 238 for op in ex.operands(): 239 result = result or is_piecewise(op) 240 return result 241 return is_piecewise(ex) 605 242 606 ::607 243 608 sage: R.<x> = QQ[] 609 sage: f1 = x^0 610 sage: f2 = 10*x  x^2 611 sage: f3 = 3*x^4  156*x^3 + 3036*x^2  26208*x 612 sage: f = Piecewise([[(0,3),f1],[(3,10),f2],[(10,20),f3]]) 613 sage: f.base_ring() 614 Rational Field 244 @staticmethod 245 def simplify(): 615 246 """ 616 return (self.functions()[0]).base_ring() 617 618 def end_points(self): 619 """ 620 Returns a list of all interval endpoints for this function. 621 622 EXAMPLES:: 623 624 sage: f1(x) = 1 625 sage: f2(x) = 1x 626 sage: f3(x) = x^25 627 sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3]]) 628 sage: f.end_points() 629 [0, 1, 2, 3] 630 """ 631 intervals = self.intervals() 632 return [ intervals[0][0] ] + [b for a,b in intervals] 633 634 def __call__(self,x0): 635 """ 636 Evaluates self at x0. Returns the average value of the jump if x0 637 is an interior endpoint of one of the intervals of self and the 638 usual value otherwise. 639 640 EXAMPLES:: 641 642 sage: f1(x) = 1 643 sage: f2(x) = 1x 644 sage: f3(x) = exp(x) 645 sage: f4(x) = sin(2*x) 646 sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]]) 647 sage: f(0.5) 648 1 649 sage: f(5/2) 650 e^(5/2) 651 sage: f(5/2).n() 652 12.1824939607035 653 sage: f(1) 654 1/2 655 """ 656 #x0 = QQ(x0) ## does not allow for evaluation at pi 657 n = self.length() 658 endpts = self.end_points() 659 for i in range(1,n): 660 if x0 == endpts[i]: 661 return (self.functions()[i1](x0) + self.functions()[i](x0))/2 662 if x0 == endpts[0]: 663 return self.functions()[0](x0) 664 if x0 == endpts[n]: 665 return self.functions()[n1](x0) 666 for i in range(n): 667 if endpts[i] < x0 < endpts[i+1]: 668 return self.functions()[i](x0) 669 raise ValueError,"Value not defined outside of domain." 670 671 def which_function(self,x0): 672 """ 673 Returns the function piece used to evaluate self at x0. 674 675 EXAMPLES:: 676 677 sage: f1(z) = z 678 sage: f2(x) = 1x 679 sage: f3(y) = exp(y) 680 sage: f4(t) = sin(2*t) 681 sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]]) 682 sage: f.which_function(3/2) 683 x > x + 1 684 """ 685 for (a,b), f in self.list(): 686 if a <= x0 <= b: 687 return f 688 raise ValueError,"Function not defined outside of domain." 689 690 def default_variable(self): 691 r""" 692 Return the default variable. The default variable is defined as the 693 first variable in the first piece that has a variable. If no pieces have 694 a variable (each piece is a constant value), `x` is returned. 695 696 The result is cached. 697 698 AUTHOR: Paul Butler 699 700 EXAMPLES:: 701 702 sage: f1(x) = 1 703 sage: f2(x) = 5*x 704 sage: p = Piecewise([[(0,1),f1],[(1,4),f2]]) 705 sage: p.default_variable() 706 x 707 708 sage: f1 = 3*var('y') 709 sage: p = Piecewise([[(0,1),4],[(1,4),f1]]) 710 sage: p.default_variable() 711 y 712 713 """ 714 try: 715 return self.__default_variable 716 except AttributeError: 717 pass 718 for _, fun in self._list: 719 try: 720 fun = SR(fun) 721 if fun.variables(): 722 v = fun.variables()[0] 723 self.__default_variable = v 724 return v 725 except TypeError: 726 # pass if fun is lambda function 727 pass 728 # default to x 729 v = var('x') 730 self.__default_value = v 731 return v 732 733 def integral(self, x=None, a=None, b=None, definite=False): 734 r""" 735 By default, returns the indefinite integral of the function. 736 If definite=True is given, returns the definite integral. 737 738 AUTHOR: 739 740  Paul Butler 741 742 EXAMPLES:: 743 744 sage: f1(x) = 1x 745 sage: f = Piecewise([[(0,1),1],[(1,2),f1]]) 746 sage: f.integral(definite=True) 747 1/2 748 749 :: 750 751 sage: f1(x) = 1 752 sage: f2(x) = 2 753 sage: f = Piecewise([[(0,pi/2),f1],[(pi/2,pi),f2]]) 754 sage: f.integral(definite=True) 755 1/2*pi 756 757 sage: f1(x) = 2 758 sage: f2(x) = 3  x 759 sage: f = Piecewise([[(2, 0), f1], [(0, 3), f2]]) 760 sage: f.integral() 761 Piecewise defined function with 2 parts, [[(2, 0), x > 2*x + 4], [(0, 3), x > 1/2*x^2 + 3*x + 4]] 762 763 sage: f1(y) = 1 764 sage: f2(y) = y + 3 765 sage: f3(y) = y  1 766 sage: f4(y) = y^2  1 767 sage: f5(y) = 3 768 sage: f = Piecewise([[(4,3),f1],[(3,2),f2],[(2,0),f3],[(0,2),f4],[(2,3),f5]]) 769 sage: F = f.integral(y) 770 sage: F 771 Piecewise defined function with 5 parts, [[(4, 3), y > y  4], [(3, 2), y > 1/2*y^2 + 3*y + 7/2], [(2, 0), y > 1/2*y^2  y  1/2], [(0, 2), y > 1/3*y^3  y  1/2], [(2, 3), y > 3*y  35/6]] 772 773 Ensure results are consistent with FTC:: 774 775 sage: F(3)  F(4) 776 1 777 sage: F(1)  F(3) 778 1 779 sage: F(2)  F(0) 780 2/3 781 sage: f.integral(y, 0, 2) 782 2/3 783 sage: F(3)  F(4) 784 19/6 785 sage: f.integral(y, 4, 3) 786 19/6 787 sage: f.integral(definite=True) 788 19/6 789 790 :: 791 792 sage: f1(y) = (y+3)^2 793 sage: f2(y) = y+3 794 sage: f3(y) = 3 795 sage: f = Piecewise([[(infinity, 3), f1], [(3, 0), f2], [(0, infinity), f3]]) 796 sage: f.integral() 797 Piecewise defined function with 3 parts, [[(Infinity, 3), y > 1/3*y^3 + 3*y^2 + 9*y + 9], [(3, 0), y > 1/2*y^2 + 3*y + 9/2], [(0, +Infinity), y > 3*y + 9/2]] 798 799 :: 800 801 sage: f1(x) = e^(abs(x)) 802 sage: f = Piecewise([[(infinity, infinity), f1]]) 803 sage: f.integral(definite=True) 804 2 805 sage: f.integral() 806 Piecewise defined function with 1 parts, [[(Infinity, +Infinity), x > 1/2*((sgn(x)  1)*e^(2*x)  2*e^x*sgn(x) + sgn(x) + 1)*e^(x)  1]] 807 808 :: 809 810 sage: f = Piecewise([((0, 5), cos(x))]) 811 sage: f.integral() 812 Piecewise defined function with 1 parts, [[(0, 5), x > sin(x)]] 813 814 815 TESTS: 816 817 Verify that piecewise integrals of zero work (trac #10841):: 818 819 sage: f0(x) = 0 820 sage: f = Piecewise([[(0,1),f0]]) 821 sage: f.integral(x,0,1) 822 0 823 sage: f = Piecewise([[(0,1), 0]]) 824 sage: f.integral(x,0,1) 825 0 826 sage: f = Piecewise([[(0,1), SR(0)]]) 827 sage: f.integral(x,0,1) 828 0 829 830 """ 831 if a != None and b != None: 832 F = self.integral(x) 833 return F(b)  F(a) 834 835 if a != None or b != None: 836 raise TypeError, 'only one endpoint given' 837 838 area = 0 # cumulative definite integral of parts to the left of the current interval 839 integrand_pieces = self.list() 840 integrand_pieces.sort() 841 new_pieces = [] 842 843 if x == None: 844 x = self.default_variable() 845 846 # The integral is computed by iterating over the pieces in order. 847 # The definite integral for each piece is calculated and accumulated in `area`. 848 # Thus at any time, `area` represents the definite integral of all the pieces 849 # encountered so far. The indefinite integral of each piece is also calculated, 850 # and the `area` before each piece is added to the piece. 851 # 852 # If a definite integral is requested, `area` is returned. 853 # Otherwise, a piecewise function is constructed from the indefinite integrals 854 # and returned. 855 # 856 # An exception is made if integral is called on a piecewise function 857 # that starts at infinity. In this case, we do not try to calculate the 858 # definite integral of the first piece, and the value of `area` remains 0 859 # after the first piece. 860 861 for (start, end), fun in integrand_pieces: 862 fun = SR(fun) 863 if start == infinity and not definite: 864 fun_integrated = fun.integral(x, end, x) 865 else: 866 try: 867 assume(start < x) 868 except ValueError: # Assumption is redundant 869 pass 870 fun_integrated = fun.integral(x, start, x) + area 871 forget(start < x) 872 if definite or end != infinity: 873 area += fun.integral(x, start, end) 874 new_pieces.append([(start, end), SR(fun_integrated).function(x)]) 875 876 if definite: 877 return SR(area) 878 else: 879 return Piecewise(new_pieces) 880 881 def convolution(self, other): 882 """ 883 Returns the convolution function, 884 `f*g(t)=\int_{\infty}^\infty f(u)g(tu)du`, for compactly 885 supported `f,g`. 886 887 EXAMPLES:: 888 889 sage: x = PolynomialRing(QQ,'x').gen() 890 sage: f = Piecewise([[(0,1),1*x^0]]) ## example 0 891 sage: g = f.convolution(f) 892 sage: h = f.convolution(g) 893 sage: P = f.plot(); Q = g.plot(rgbcolor=(1,1,0)); R = h.plot(rgbcolor=(0,1,1)); 894 sage: # Type show(P+Q+R) to view 895 sage: f = Piecewise([[(0,1),1*x^0],[(1,2),2*x^0],[(2,3),1*x^0]]) ## example 1 896 sage: g = f.convolution(f) 897 sage: h = f.convolution(g) 898 sage: P = f.plot(); Q = g.plot(rgbcolor=(1,1,0)); R = h.plot(rgbcolor=(0,1,1)); 899 sage: # Type show(P+Q+R) to view 900 sage: f = Piecewise([[(1,1),1]]) ## example 2 901 sage: g = Piecewise([[(0,3),x]]) 902 sage: f.convolution(g) 903 Piecewise defined function with 3 parts, [[(1, 1), 0], [(1, 2), 3/2*x], [(2, 4), 3/2*x]] 904 sage: g = Piecewise([[(0,3),1*x^0],[(3,4),2*x^0]]) 905 sage: f.convolution(g) 906 Piecewise defined function with 5 parts, [[(1, 1), x + 1], [(1, 2), 3], [(2, 3), x], [(3, 4), x + 8], [(4, 5), 2*x + 10]] 907 """ 908 f = self 909 g = other 910 M = min(min(f.end_points()),min(g.end_points())) 911 N = max(max(f.end_points()),max(g.end_points())) 912 R2 = PolynomialRing(QQ,2,names=["tt","uu"]) 913 tt,uu = R2.gens() 914 conv = 0 915 f0 = f.functions()[0] 916 g0 = g.functions()[0] 917 R1 = f0.parent() 918 xx = R1.gen() 919 var = repr(xx) 920 if len(f.intervals())==1 and len(g.intervals())==1: 921 f = f.unextend() 922 g = g.unextend() 923 a1 = f.intervals()[0][0] 924 a2 = f.intervals()[0][1] 925 b1 = g.intervals()[0][0] 926 b2 = g.intervals()[0][1] 927 i1 = repr(f0).replace(var,repr(uu)) 928 i2 = repr(g0).replace(var,"("+repr(ttuu)+")") 929 cmd1 = "integrate((%s)*(%s),%s,%s,%s)"%(i1,i2, uu, a1, ttb1) ## if a1+b1 < tt < a2+b1 930 cmd2 = "integrate((%s)*(%s),%s,%s,%s)"%(i1,i2, uu, ttb2, ttb1) ## if a1+b2 < tt < a2+b1 931 cmd3 = "integrate((%s)*(%s),%s,%s,%s)"%(i1,i2, uu, ttb2, a2) ## if a1+b2 < tt < a2+b2 932 cmd4 = "integrate((%s)*(%s),%s,%s,%s)"%(i1,i2, uu, a1, a2) ## if a2+b1 < tt < a1+b2 933 conv1 = maxima.eval(cmd1) 934 conv2 = maxima.eval(cmd2) 935 conv3 = maxima.eval(cmd3) 936 conv4 = maxima.eval(cmd4) 937 # this is a very, very, very ugly hack 938 x = PolynomialRing(QQ,'x').gen() 939 fg1 = sage_eval(conv1.replace("tt",var), {'x':x}) ## should be = R2(conv1) 940 fg2 = sage_eval(conv2.replace("tt",var), {'x':x}) ## should be = R2(conv2) 941 fg3 = sage_eval(conv3.replace("tt",var), {'x':x}) ## should be = R2(conv3) 942 fg4 = sage_eval(conv4.replace("tt",var), {'x':x}) ## should be = R2(conv4) 943 if a1b1<a2b2: 944 if a2+b1!=a1+b2: 945 h = Piecewise([[(a1+b1,a1+b2),fg1],[(a1+b2,a2+b1),fg4],[(a2+b1,a2+b2),fg3]]) 946 else: 947 h = Piecewise([[(a1+b1,a1+b2),fg1],[(a1+b2,a2+b2),fg3]]) 948 else: 949 if a1+b2!=a2+b1: 950 h = Piecewise([[(a1+b1,a2+b1),fg1],[(a2+b1,a1+b2),fg2],[(a1+b2,a2+b2),fg3]]) 951 else: 952 h = Piecewise([[(a1+b1,a2+b1),fg1],[(a2+b1,a2+b2),fg3]]) 953 return h 954 955 if len(f.intervals())>1 or len(g.intervals())>1: 956 z = Piecewise([[(3*abs(NM),3*abs(NM)),0*xx**0]]) 957 ff = f.functions() 958 gg = g.functions() 959 intvlsf = f.intervals() 960 intvlsg = g.intervals() 961 for i in range(len(ff)): 962 for j in range(len(gg)): 963 f0 = Piecewise([[intvlsf[i],ff[i]]]) 964 g0 = Piecewise([[intvlsg[j],gg[j]]]) 965 h = g0.convolution(f0) 966 z = z + h 967 return z.unextend() 968 969 def derivative(self): 970 r""" 971 Returns the derivative (as computed by maxima) 972 Piecewise(I,`(d/dx)(self_I)`), as I runs over the 973 intervals belonging to self. self must be piecewise polynomial. 974 975 EXAMPLES:: 976 977 sage: f1(x) = 1 978 sage: f2(x) = 1x 979 sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]) 980 sage: f.derivative() 981 Piecewise defined function with 2 parts, [[(0, 1), x > 0], [(1, 2), x > 1]] 982 sage: f1(x) = 1 983 sage: f2(x) = 2 984 sage: f = Piecewise([[(0,pi/2),f1],[(pi/2,pi),f2]]) 985 sage: f.derivative() 986 Piecewise defined function with 2 parts, [[(0, 1/2*pi), x > 0], [(1/2*pi, pi), x > 0]] 987 988 :: 989 990 sage: f = Piecewise([[(0,1), (x * 2)]], x) 991 sage: f.derivative() 992 Piecewise defined function with 1 parts, [[(0, 1), x > 2]] 993 """ 994 x = self.default_variable() 995 dlist = [[(a, b), derivative(f(x), x).function(x)] for (a,b),f in self.list()] 996 return Piecewise(dlist) 997 998 def tangent_line(self, pt): 999 """ 1000 Computes the linear function defining the tangent line of the 1001 piecewise function self. 1002 1003 EXAMPLES:: 1004 1005 sage: f1(x) = x^2 1006 sage: f2(x) = 5x^3+x 1007 sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]) 1008 sage: tf = f.tangent_line(0.9) ## tangent line at x=0.9 1009 sage: P = f.plot(rgbcolor=(0.7,0.1,0.5), plot_points=40) 1010 sage: Q = tf.plot(rgbcolor=(0.7,0.2,0.2), plot_points=40) 1011 sage: P + Q 1012 """ 1013 pt = QQ(pt) 1014 R = QQ[self.default_variable()] 1015 x = R.gen() 1016 der = self.derivative() 1017 tanline = (xpt)*der(pt)+self(pt) 1018 dlist = [[(a, b), tanline] for (a,b),f in self.list()] 1019 return Piecewise(dlist) 1020 1021 def plot(self, *args, **kwds): 1022 """ 1023 Returns the plot of self. 1024 1025 Keyword arguments are passed onto the plot command for each piece 1026 of the function. E.g., the plot_points keyword affects each 1027 segment of the plot. 1028 1029 EXAMPLES:: 1030 1031 sage: f1(x) = 1 1032 sage: f2(x) = 1x 1033 sage: f3(x) = exp(x) 1034 sage: f4(x) = sin(2*x) 1035 sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]]) 1036 sage: P = f.plot(rgbcolor=(0.7,0.1,0), plot_points=40) 1037 sage: P 1038 1039 Remember: to view this, type show(P) or P.save("path/myplot.png") 1040 and then open it in a graphics viewer such as GIMP. 1041 1042 TESTS: 1043 1044 We should not add each piece to the legend individually, since 1045 this creates duplicates (:trac:`12651`). This tests that only 1046 one of the graphics objects in the plot has a non``None`` 1047 ``legend_label``:: 1048 1049 sage: f1 = sin(x) 1050 sage: f2 = cos(x) 1051 sage: f = piecewise([[(1,0), f1],[(0,1), f2]]) 1052 sage: p = f.plot(legend_label='$f(x)$') 1053 sage: lines = [ 1054 ... line 1055 ... for line in p._objects 1056 ... if line.options()['legend_label'] is not None ] 1057 sage: len(lines) 1058 1 1059 """ 1060 from sage.plot.all import plot, Graphics 1061 1062 g = Graphics() 1063 1064 for i, ((a,b), f) in enumerate(self.list()): 1065 # If it's the first piece, pass all arguments. Otherwise, 1066 # filter out 'legend_label' so that we don't add each 1067 # piece to the legend separately (trac #12651). 1068 if i != 0 and 'legend_label' in kwds: 1069 del kwds['legend_label'] 1070 1071 g += plot(f, a, b, *args, **kwds) 1072 1073 return g 1074 1075 def fourier_series_cosine_coefficient(self,n,L): 1076 r""" 1077 Returns the nth Fourier series coefficient of 1078 `\cos(n\pi x/L)`, `a_n`. 1079 1080 INPUT: 1081 1082 1083  ``self``  the function f(x), defined over L x L 1084 1085  ``n``  an integer n=0 1086 1087  ``L``  (the period)/2 1088 1089 1090 OUTPUT: 1091 `a_n = \frac{1}{L}\int_{L}^L f(x)\cos(n\pi x/L)dx` 1092 1093 EXAMPLES:: 1094 1095 sage: f(x) = x^2 1096 sage: f = Piecewise([[(1,1),f]]) 1097 sage: f.fourier_series_cosine_coefficient(2,1) 1098 pi^(2) 1099 sage: f(x) = x^2 1100 sage: f = Piecewise([[(pi,pi),f]]) 1101 sage: f.fourier_series_cosine_coefficient(2,pi) 1102 1 1103 sage: f1(x) = 1 1104 sage: f2(x) = 2 1105 sage: f = Piecewise([[(pi,pi/2),f1],[(pi/2,pi),f2]]) 1106 sage: f.fourier_series_cosine_coefficient(5,pi) 1107 3/5/pi 1108 """ 1109 from sage.all import cos, pi 1110 x = var('x') 1111 result = sum([(f(x)*cos(pi*x*n/L)/L).integrate(x, a, b) 1112 for (a,b), f in self.list()]) 1113 if is_Expression(result): 1114 return result.simplify_trig() 1115 return result 1116 1117 def fourier_series_sine_coefficient(self,n,L): 1118 r""" 1119 Returns the nth Fourier series coefficient of 1120 `\sin(n\pi x/L)`, `b_n`. 1121 1122 INPUT: 1123 1124 1125  ``self``  the function f(x), defined over L x L 1126 1127  ``n``  an integer n0 1128 1129  ``L``  (the period)/2 1130 1131 1132 OUTPUT: 1133 `b_n = \frac{1}{L}\int_{L}^L f(x)\sin(n\pi x/L)dx` 1134 1135 EXAMPLES:: 1136 1137 sage: f(x) = x^2 1138 sage: f = Piecewise([[(1,1),f]]) 1139 sage: f.fourier_series_sine_coefficient(2,1) # L=1, n=2 1140 0 1141 """ 1142 from sage.all import sin, pi 1143 x = var('x') 1144 result = sum([(f(x)*sin(pi*x*n/L)/L).integrate(x, a, b) 1145 for (a,b), f in self.list()]) 1146 if is_Expression(result): 1147 return result.simplify_trig() 1148 return result 1149 1150 def _fourier_series_helper(self, N, L, scale_function): 1151 r""" 1152 A helper function for the construction of Fourier series. The 1153 argument scale_function is a function which takes in n, 1154 representing the `n^{th}` coefficient, and return an 1155 expression to scale the sine and cosine coefficients by. 1156 1157 EXAMPLES:: 1158 1159 sage: f(x) = x^2 1160 sage: f = Piecewise([[(1,1),f]]) 1161 sage: f._fourier_series_helper(3, 1, lambda n: 1) 1162 cos(2*pi*x)/pi^2  4*cos(pi*x)/pi^2 + 1/3 1163 """ 1164 from sage.all import pi, sin, cos, srange 1165 x = self.default_variable() 1166 a0 = self.fourier_series_cosine_coefficient(0,L) 1167 result = a0/2 + sum([(self.fourier_series_cosine_coefficient(n,L)*cos(n*pi*x/L) + 1168 self.fourier_series_sine_coefficient(n,L)*sin(n*pi*x/L))* 1169 scale_function(n) 1170 for n in srange(1,N)]) 1171 return result.expand() 1172 1173 1174 def fourier_series_partial_sum(self,N,L): 1175 r""" 1176 Returns the partial sum 1177 1178 .. math:: 1179 1180 f(x) \sim \frac{a_0}{2} + \sum_{n=1}^N [a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})], 1181 1182 as a string. 1183 1184 EXAMPLE:: 1185 1186 sage: f(x) = x^2 1187 sage: f = Piecewise([[(1,1),f]]) 1188 sage: f.fourier_series_partial_sum(3,1) 1189 cos(2*pi*x)/pi^2  4*cos(pi*x)/pi^2 + 1/3 1190 sage: f1(x) = 1 1191 sage: f2(x) = 2 1192 sage: f = Piecewise([[(pi,pi/2),f1],[(pi/2,pi),f2]]) 1193 sage: f.fourier_series_partial_sum(3,pi) 1194 3*cos(x)/pi  3*sin(2*x)/pi + 3*sin(x)/pi  1/4 1195 """ 1196 return self._fourier_series_helper(N, L, lambda n: 1) 1197 1198 def fourier_series_partial_sum_cesaro(self,N,L): 1199 r""" 1200 Returns the Cesaro partial sum 1201 1202 .. math:: 1203 1204 f(x) \sim \frac{a_0}{2} + \sum_{n=1}^N (1n/N)*[a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})], 1205 1206 1207 as a string. This is a "smoother" partial sum  the Gibbs 1208 phenomenon is mollified. 1209 1210 EXAMPLE:: 1211 1212 sage: f(x) = x^2 1213 sage: f = Piecewise([[(1,1),f]]) 1214 sage: f.fourier_series_partial_sum_cesaro(3,1) 1215 1/3*cos(2*pi*x)/pi^2  8/3*cos(pi*x)/pi^2 + 1/3 1216 sage: f1(x) = 1 1217 sage: f2(x) = 2 1218 sage: f = Piecewise([[(pi,pi/2),f1],[(pi/2,pi),f2]]) 1219 sage: f.fourier_series_partial_sum_cesaro(3,pi) 1220 2*cos(x)/pi  sin(2*x)/pi + 2*sin(x)/pi  1/4 1221 """ 1222 return self._fourier_series_helper(N, L, lambda n: 1n/N) 1223 1224 def fourier_series_partial_sum_hann(self,N,L): 1225 r""" 1226 Returns the Hannfiltered partial sum (named after von Hann, not 1227 Hamming) 1228 1229 .. math:: 1230 1231 f(x) \sim \frac{a_0}{2} + \sum_{n=1}^N H_N(n)*[a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})], 1232 1233 as a string, where `H_N(x) = (1+\cos(\pi x/N))/2`. This is 1234 a "smoother" partial sum  the Gibbs phenomenon is mollified. 1235 1236 EXAMPLE:: 1237 1238 sage: f(x) = x^2 1239 sage: f = Piecewise([[(1,1),f]]) 1240 sage: f.fourier_series_partial_sum_hann(3,1) 1241 1/4*cos(2*pi*x)/pi^2  3*cos(pi*x)/pi^2 + 1/3 1242 sage: f1(x) = 1 1243 sage: f2(x) = 2 1244 sage: f = Piecewise([[(pi,pi/2),f1],[(pi/2,pi),f2]]) 1245 sage: f.fourier_series_partial_sum_hann(3,pi) 1246 9/4*cos(x)/pi  3/4*sin(2*x)/pi + 9/4*sin(x)/pi  1/4 1247 """ 1248 from sage.all import cos, pi 1249 return self._fourier_series_helper(N, L, lambda n: (1+cos(pi*n/N))/2) 1250 1251 def fourier_series_partial_sum_filtered(self,N,L,F): 1252 r""" 1253 Returns the "filtered" partial sum 1254 1255 .. math:: 1256 1257 f(x) \sim \frac{a_0}{2} + \sum_{n=1}^N F_n*[a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})], 1258 1259 as a string, where `F = [F_1,F_2, ..., F_{N}]` is a list 1260 of length `N` consisting of real numbers. This can be used 1261 to plot FS solutions to the heat and wave PDEs. 1262 1263 EXAMPLE:: 1264 1265 sage: f(x) = x^2 1266 sage: f = Piecewise([[(1,1),f]]) 1267 sage: f.fourier_series_partial_sum_filtered(3,1,[1,1,1]) 1268 cos(2*pi*x)/pi^2  4*cos(pi*x)/pi^2 + 1/3 1269 sage: f1(x) = 1 1270 sage: f2(x) = 2 1271 sage: f = Piecewise([[(pi,pi/2),f1],[(pi/2,pi),f2]]) 1272 sage: f.fourier_series_partial_sum_filtered(3,pi,[1,1,1]) 1273 3*cos(x)/pi  3*sin(2*x)/pi + 3*sin(x)/pi  1/4 1274 """ 1275 return self._fourier_series_helper(N, L, lambda n: F[n]) 1276 1277 def plot_fourier_series_partial_sum(self,N,L,xmin,xmax, **kwds): 1278 r""" 1279 Plots the partial sum 1280 1281 .. math:: 1282 1283 f(x) \sim \frac{a_0}{2} + sum_{n=1}^N [a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})], 1284 1285 over xmin x xmin. 1286 1287 EXAMPLE:: 1288 1289 sage: f1(x) = 2 1290 sage: f2(x) = 1 1291 sage: f3(x) = 1 1292 sage: f4(x) = 2 1293 sage: f = Piecewise([[(pi,pi/2),f1],[(pi/2,0),f2],[(0,pi/2),f3],[(pi/2,pi),f4]]) 1294 sage: P = f.plot_fourier_series_partial_sum(3,pi,5,5) # long time 1295 sage: f1(x) = 1 1296 sage: f2(x) = 2 1297 sage: f = Piecewise([[(pi,pi/2),f1],[(pi/2,pi),f2]]) 1298 sage: P = f.plot_fourier_series_partial_sum(15,pi,5,5) # long time 1299 1300 Remember, to view this type show(P) or P.save("path/myplot.png") 1301 and then open it in a graphics viewer such as GIMP. 1302 """ 1303 from sage.plot.all import plot 1304 return plot(self.fourier_series_partial_sum(N,L), xmin, xmax, **kwds) 1305 1306 def plot_fourier_series_partial_sum_cesaro(self,N,L,xmin,xmax, **kwds): 1307 r""" 1308 Plots the partial sum 1309 1310 .. math:: 1311 1312 f(x) \sim \frac{a_0}{2} + \sum_{n=1}^N (1n/N)*[a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})], 1313 1314 1315 over xmin x xmin. This is a "smoother" partial sum  the Gibbs 1316 phenomenon is mollified. 1317 1318 EXAMPLE:: 1319 1320 sage: f1(x) = 2 1321 sage: f2(x) = 1 1322 sage: f3(x) = 1 1323 sage: f4(x) = 2 1324 sage: f = Piecewise([[(pi,pi/2),f1],[(pi/2,0),f2],[(0,pi/2),f3],[(pi/2,pi),f4]]) 1325 sage: P = f.plot_fourier_series_partial_sum_cesaro(3,pi,5,5) # long time 1326 sage: f1(x) = 1 1327 sage: f2(x) = 2 1328 sage: f = Piecewise([[(pi,pi/2),f1],[(pi/2,pi),f2]]) 1329 sage: P = f.plot_fourier_series_partial_sum_cesaro(15,pi,5,5) # long time 1330 1331 Remember, to view this type show(P) or P.save("path/myplot.png") 1332 and then open it in a graphics viewer such as GIMP. 1333 """ 1334 from sage.plot.all import plot 1335 return plot(self.fourier_series_partial_sum_cesaro(N,L), xmin, xmax, **kwds) 1336 1337 def plot_fourier_series_partial_sum_hann(self,N,L,xmin,xmax, **kwds): 1338 r""" 1339 Plots the partial sum 1340 1341 .. math:: 1342 1343 f(x) \sim \frac{a_0}{2} + \sum_{n=1}^N H_N(n)*[a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})], 1344 1345 1346 over xmin x xmin, where H_N(x) = (0.5)+(0.5)\*cos(x\*pi/N) is the 1347 Nth Hann filter. 1348 1349 EXAMPLE:: 1350 1351 sage: f1(x) = 2 1352 sage: f2(x) = 1 1353 sage: f3(x) = 1 1354 sage: f4(x) = 2 1355 sage: f = Piecewise([[(pi,pi/2),f1],[(pi/2,0),f2],[(0,pi/2),f3],[(pi/2,pi),f4]]) 1356 sage: P = f.plot_fourier_series_partial_sum_hann(3,pi,5,5) # long time 1357 sage: f1(x) = 1 1358 sage: f2(x) = 2 1359 sage: f = Piecewise([[(pi,pi/2),f1],[(pi/2,pi),f2]]) 1360 sage: P = f.plot_fourier_series_partial_sum_hann(15,pi,5,5) # long time 1361 1362 Remember, to view this type show(P) or P.save("path/myplot.png") 1363 and then open it in a graphics viewer such as GIMP. 1364 """ 1365 from sage.plot.all import plot 1366 return plot(self.fourier_series_partial_sum_hann(N,L), xmin, xmax, **kwds) 1367 1368 def plot_fourier_series_partial_sum_filtered(self,N,L,F,xmin,xmax, **kwds): 1369 r""" 1370 Plots the partial sum 1371 1372 .. math:: 1373 1374 f(x) \sim \frac{a_0}{2} + \sum_{n=1}^N F_n*[a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})], 1375 1376 1377 over xmin x xmin, where `F = [F_1,F_2, ..., F_{N}]` is a 1378 list of length `N` consisting of real numbers. This can be 1379 used to plot FS solutions to the heat and wave PDEs. 1380 1381 EXAMPLE:: 1382 1383 sage: f1(x) = 2 1384 sage: f2(x) = 1 1385 sage: f3(x) = 1 1386 sage: f4(x) = 2 1387 sage: f = Piecewise([[(pi,pi/2),f1],[(pi/2,0),f2],[(0,pi/2),f3],[(pi/2,pi),f4]]) 1388 sage: P = f.plot_fourier_series_partial_sum_filtered(3,pi,[1]*3,5,5) # long time 1389 sage: f1(x) = 1 1390 sage: f2(x) = 2 1391 sage: f = Piecewise([[(pi,pi/2),f1],[(pi/2,0),f2],[(0,pi/2),f1],[(pi/2,pi),f2]]) 1392 sage: P = f.plot_fourier_series_partial_sum_filtered(15,pi,[1]*15,5,5) # long time 1393 1394 Remember, to view this type show(P) or P.save("path/myplot.png") 1395 and then open it in a graphics viewer such as GIMP. 1396 """ 1397 from sage.plot.all import plot 1398 return plot(self.fourier_series_partial_sum_filtered(N,L,F), xmin, xmax, **kwds) 1399 1400 def fourier_series_value(self,x,L): 1401 r""" 1402 Returns the value of the Fourier series coefficient of self at 1403 `x`, 1404 1405 1406 .. math:: 1407 1408 f(x) \sim \frac{a_0}{2} + \sum_{n=1}^\infty [a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})], \ \ \ L<x<L. 1409 1410 1411 This method applies to piecewise nonpolynomial functions as well. 1412 1413 INPUT: 1414 1415 1416  ``self``  the function f(x), defined over L x L 1417 1418  ``x``  a real number 1419 1420  ``L``  (the period)/2 1421 1422 1423 OUTPUT: `(f^*(x+)+f^*(x)/2`, where `f^*` denotes 1424 the function `f` extended to `\RR` with period 1425 `2L` (Dirichlet's Theorem for Fourier series). 1426 1427 EXAMPLES:: 1428 1429 sage: f1(x) = 1 1430 sage: f2(x) = 1x 1431 sage: f3(x) = exp(x) 1432 sage: f4(x) = sin(2*x) 1433 sage: f = Piecewise([[(10,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]]) 1434 sage: f.fourier_series_value(101,10) 1435 1/2 1436 sage: f.fourier_series_value(100,10) 1437 1 1438 sage: f.fourier_series_value(10,10) 1439 1/2*sin(20) + 1/2 1440 sage: f.fourier_series_value(20,10) 1441 1 1442 sage: f.fourier_series_value(30,10) 1443 1/2*sin(20) + 1/2 1444 sage: f1(x) = 1 1445 sage: f2(x) = 2 1446 sage: f = Piecewise([[(pi,0),lambda x:0],[(0,pi/2),f1],[(pi/2,pi),f2]]) 1447 sage: f.fourier_series_value(1,pi) 1448 0 1449 sage: f.fourier_series_value(20,pi) 1450 1 1451 sage: f.fourier_series_value(pi/2,pi) 1452 1/2 1453 """ 1454 xnew = x  int(RR(x/(2*L)))*2*L 1455 endpts = self.end_points() 1456 if xnew == endpts[0] or xnew == endpts[1]: 1457 return (self.functions()[0](endpts[0]) + self.functions()[1](endpts[1]))/2 1458 else: 1459 return self(xnew) 1460 1461 def cosine_series_coefficient(self,n,L): 1462 r""" 1463 Returns the nth cosine series coefficient of 1464 `\cos(n\pi x/L)`, `a_n`. 1465 1466 INPUT: 1467 1468 1469  ``self``  the function f(x), defined over 0 x L (no 1470 checking is done to insure this) 1471 1472  ``n``  an integer n=0 1473 1474  ``L``  (the period)/2 1475 1476 1477 OUTPUT: 1478 `a_n = \frac{2}{L}\int_{L}^L f(x)\cos(n\pi x/L)dx` such 1479 that 1480 1481 .. math:: 1482 1483 f(x) \sim \frac{a_0}{2} + \sum_{n=1}^\infty a_n\cos(\frac{n\pi x}{L}),\ \ 0<x<L. 1484 1485 1486 1487 EXAMPLES:: 1488 1489 sage: f(x) = x 1490 sage: f = Piecewise([[(0,1),f]]) 1491 sage: f.cosine_series_coefficient(2,1) 1492 0 1493 sage: f.cosine_series_coefficient(3,1) 1494 4/9/pi^2 1495 sage: f1(x) = 1 1496 sage: f2(x) = 2 1497 sage: f = Piecewise([[(0,pi/2),f1],[(pi/2,pi),f2]]) 1498 sage: f.cosine_series_coefficient(2,pi) 1499 0 1500 sage: f.cosine_series_coefficient(3,pi) 1501 2/pi 1502 sage: f.cosine_series_coefficient(111,pi) 1503 2/37/pi 1504 sage: f1 = lambda x: x*(pix) 1505 sage: f = Piecewise([[(0,pi),f1]]) 1506 sage: f.cosine_series_coefficient(0,pi) 1507 1/3*pi^2 1508 1509 """ 1510 from sage.all import cos, pi 1511 x = var('x') 1512 result = sum([(2*f(x)*cos(pi*x*n/L)/L).integrate(x, a, b) 1513 for (a,b), f in self.list()]) 1514 if is_Expression(result): 1515 return result.simplify_trig() 1516 return result 1517 1518 1519 def sine_series_coefficient(self,n,L): 1520 r""" 1521 Returns the nth sine series coefficient of 1522 `\sin(n\pi x/L)`, `b_n`. 1523 1524 INPUT: 1525 1526  ``self``  the function f(x), defined over 0 x L (no 1527 checking is done to insure this) 1528 1529  ``n``  an integer n0 1530 1531  ``L``  (the period)/2 247 Combine piecewise operands into single piecewise function 1532 248 1533 249 OUTPUT: 1534 250 1535 `b_n = \frac{2}{L}\int_{L}^L f(x)\sin(n\pi x/L)dx` such 1536 that 251 A piecewise function whose operands are not piecewiese if 252 possible, that is, as long as the piecewise variable is the same. 253 254 """ 255 raise NotImplementedError 256 257 258 class EvaluationMethods: 259 260 def expression_at(cls, self, parameters, variable, point): 261 """ 262 Return the expression defining the piecewise function at 263 ``value`` 264 265 INPUT: 266 267  ``point``  a real number. 268 269 OUTPUT: 270 271 The symbolic expression defining the function value at the 272 given ``point``. 273 274 EXAMPLES:: 275 276 sage: f = piecewise([([0,0], sin(x)), ((0,2), cos(x))]); f 277 piecewise(x>sin(x) on {0}, x>cos(x) on (0, 2); x) 278 sage: f.expression_at(0) 279 sin(x) 280 sage: f.expression_at(1) 281 cos(x) 282 sage: f.expression_at(2) 283 Traceback (most recent call last): 284 ... 285 ValueError: point is not in the domain 286 """ 287 for domain, func in parameters: 288 if domain.contains(point): 289 return func 290 raise ValueError('point is not in the domain') 291 292 which_function = expression_at 293 294 def domains(cls, self, parameters, variable): 295 """ 296 Return the individual domains 297 298 See also :meth:`~expressions`. 299 300 OUTPUT: 301 302 The collection of domains of the component functions as a 303 tuple of :class:`~sage.sets.real_set.RealSet`. 304 305 EXAMPLES:: 306 307 sage: f = piecewise([([0,0], sin(x)), ((0,2), cos(x))]); f 308 piecewise(x>sin(x) on {0}, x>cos(x) on (0, 2); x) 309 sage: f.domains() 310 ({0}, (0, 2)) 311 """ 312 return tuple(dom for dom, fun in parameters) 313 314 def domain(cls, self, parameters, variable): 315 """ 316 Return the domain 317 318 OUTPUT: 319 320 The union of the domains of the individual pieces as a 321 :class:`~sage.sets.real_set.RealSet`. 322 323 EXAMPLES:: 324 325 sage: f = piecewise([([0,0], sin(x)), ((0,2), cos(x))]); f 326 piecewise(x>sin(x) on {0}, x>cos(x) on (0, 2); x) 327 sage: f.domain() 328 [0, 2) 329 """ 330 intervals = [] 331 for domain, func in parameters: 332 intervals += list(domain) 333 return RealSet(*intervals) 334 335 def __len__(cls, self, parameters, variable): 336 """ 337 Return the number of "pieces" 338 339 OUTPUT: 340 341 Integer. 342 343 EXAMPLES:: 344 345 sage: f = piecewise([([0,0], sin(x)), ((0,2), cos(x))]); f 346 piecewise(x>sin(x) on {0}, x>cos(x) on (0, 2); x) 347 sage: len(f) 348 2 349 """ 350 return len(parameters) 351 352 def expressions(cls, self, parameters, variable): 353 """ 354 Return the individual domains 355 356 See also :meth:`~domains`. 357 358 OUTPUT: 359 360 The collection of expressions of the component functions. 361 362 EXAMPLES:: 363 364 sage: f = piecewise([([0,0], sin(x)), ((0,2), cos(x))]); f 365 piecewise(x>sin(x) on {0}, x>cos(x) on (0, 2); x) 366 sage: f.expressions() 367 (sin(x), cos(x)) 368 """ 369 return tuple(fun for dom, fun in parameters) 370 371 def iteritems(cls, self, parameters, variable): 372 for pair in parameters: 373 yield pair 374 375 def __call__(cls, self, parameters, variable, value=None, **kwds): 376 """ 377 Call the piecewise function 378 379 EXAMPLES:: 380 381 sage: f = piecewise([([0,0], sin(x)), ((0,2), cos(x))]); f 382 piecewise(x>sin(x) on {0}, x>cos(x) on (0, 2); x) 383 sage: f(0) 384 0 385 sage: f(1) 386 cos(1) 387 sage: f(2) 388 Traceback (most recent call last): 389 ... 390 ValueError: point is not in the domain 391 """ 392 self = piecewise(parameters, var=variable) 393 substitution = dict() 394 for k, v in kwds.iteritems(): 395 substitution[SR.var(k)] = v 396 if value is not None: 397 substitution[variable] = value 398 return self.subs(substitution) 399 400 def _fast_float_(cls, self, *args): 401 """ 402 Do not support the old ``fast_float`` 403 404 OUTPUT: 405 406 This method raises ``NotImplementedError`` so that 407 plotting uses the newer `fast_callable` implementation. 408 409 EXAMPLES:: 410 411 sage: f = piecewise([([0,0], sin(x)), ((0,2), cos(x))]) 412 sage: f._fast_float_() 413 Traceback (most recent call last): 414 ... 415 NotImplementedError 416 """ 417 raise NotImplementedError 418 419 def _fast_callable_(cls, self, parameters, variable, etb): 420 """ 421 Override the ``fast_callable`` 422 423 OUTPUT: 424 425 A :class:`~sage.ext.fast_callable.ExpressionCall` 426 representing the piecewise function in the expression 427 tree. 428 429 EXAMPLES:: 430 431 sage: p = piecewise([((1, 0), x), ([0, 1], x)], var=x) 432 sage: from sage.ext.fast_callable import ExpressionTreeBuilder 433 sage: etb = ExpressionTreeBuilder(vars=['x']) 434 sage: p._fast_callable_(etb) 435 {CommutativeRings.element_class}(v_0) 436 """ 437 # print 'ev_fast_cal', parameters, variable, etb 438 self = piecewise(parameters, var=variable) 439 return etb.call(self, variable) 440 441 def restriction(cls, self, parameters, variable, restricted_domain): 442 """ 443 Restrict the domain 444 445 INPUT: 446 447  ``restricted_domain``  a 448 :class:`~sage.sets.real_set.RealSet` or something that 449 defines one. 450 451 OUTPUT: 452 453 A new piecewise function obtained by restricting the domain. 454 455 EXAMPLES:: 456 457 sage: f = piecewise([((oo, oo), x)]); f 458 piecewise(x>x on (oo, +oo); x) 459 sage: f.restriction([[1,1], [3,3]]) 460 piecewise(x>x on [1, 1] + {3}; x) 461 """ 462 restricted_domain = RealSet(*restricted_domain) 463 new_param = [] 464 for domain, func in parameters: 465 domain = domain.intersection(restricted_domain) 466 new_param.append((domain, func)) 467 return piecewise(new_param, var=variable) 468 469 def extension(cls, self, parameters, variable, extension, extension_domain=None): 470 """ 471 Extend the function 472 473 INPUT: 474 475  ``extension``  a symbolic expression 476 477  ``extension_domain``  a 478 :class:`~sage.sets.real_set.RealSet` or ``None`` 479 (default). The domain of the extension. By default, the 480 entire complement of the current domain. 481 482 EXAMPLES:: 483 484 sage: f = piecewise([((1,1), x)]); f 485 piecewise(x>x on (1, 1); x) 486 sage: f(3) 487 Traceback (most recent call last): 488 ... 489 ValueError: point is not in the domain 490 491 sage: g = f.extension(0); g 492 piecewise(x>x on (1, 1), x>0 on (oo, 1] + [1, +oo); x) 493 sage: g(3) 494 0 495 496 sage: h = f.extension(1, RealSet.unbounded_above_closed(1)); h 497 piecewise(x>x on (1, 1), x>1 on [1, +oo); x) 498 sage: h(3) 499 1 500 """ 501 self = piecewise(parameters, var=variable) 502 if extension_domain is None: 503 extension_domain = self.domain().complement() 504 ext = ((extension_domain, SR(extension)),) 505 return piecewise(parameters + ext, var=variable) 1537 506 1538 .. math:: 1539 1540 f(x) \sim \sum_{n=1}^\infty b_n\sin(\frac{n\pi x}{L}),\ \ 0<x<L. 1541 1542 EXAMPLES:: 1543 1544 sage: f(x) = 1 1545 sage: f = Piecewise([[(0,1),f]]) 1546 sage: f.sine_series_coefficient(2,1) 1547 0 1548 sage: f.sine_series_coefficient(3,1) 1549 4/3/pi 1550 """ 1551 from sage.all import sin, pi 1552 x = var('x') 1553 result = sum([(2*f(x)*sin(pi*x*n/L)/L).integrate(x, a, b) 1554 for (a,b), f in self.list()]) 1555 if is_Expression(result): 1556 return result.simplify_trig() 1557 return result 507 def pieces(cls, self, parameters, variable): 508 """ 509 Return the "pieces" 1558 510 1559 def laplace(self, x='x', s='t'): 1560 r""" 1561 Returns the Laplace transform of self with respect to the variable 1562 var. 1563 1564 INPUT: 1565 1566 1567  ``x``  variable of self 1568 1569  ``s``  variable of Laplace transform. 1570 1571 1572 We assume that a piecewise function is 0 outside of its domain and 1573 that the leftmost endpoint of the domain is 0. 1574 1575 EXAMPLES:: 1576 1577 sage: x, s, w = var('x, s, w') 1578 sage: f = Piecewise([[(0,1),1],[(1,2), 1x]]) 1579 sage: f.laplace(x, s) 1580 e^(s)/s + (s + 1)*e^(2*s)/s^2 + 1/s  e^(s)/s^2 1581 sage: f.laplace(x, w) 1582 e^(w)/w + (w + 1)*e^(2*w)/w^2 + 1/w  e^(w)/w^2 1583 1584 :: 1585 1586 sage: y, t = var('y, t') 1587 sage: f = Piecewise([[(1,2), 1y]]) 1588 sage: f.laplace(y, t) 1589 (t + 1)*e^(2*t)/t^2  e^(t)/t^2 1590 1591 :: 1592 1593 sage: s = var('s') 1594 sage: t = var('t') 1595 sage: f1(t) = t 1596 sage: f2(t) = 2 1597 sage: f = Piecewise([[(0,1),f1],[(1,infinity),f2]]) 1598 sage: f.laplace(t,s) 1599 (s + 1)*e^(s)/s^2 + 2*e^(s)/s  1/s^2 1600 """ 1601 from sage.all import assume, exp, forget 1602 x = var(x) 1603 s = var(s) 1604 assume(s>0) 1605 result = sum([(SR(f)*exp(s*x)).integral(x,a,b) 1606 for (a,b),f in self.list()]) 1607 forget(s>0) 1608 return result 511 OUTPUT: 1609 512 1610 def _make_compatible(self, other): 1611 """ 1612 Returns self and other extended to be defined on the same domain as 1613 well as a refinement of their intervals. This is used for adding 1614 and multiplying piecewise functions. 1615 1616 EXAMPLES:: 1617 1618 sage: R.<x> = QQ[] 1619 sage: f1 = Piecewise([[(0, 2), x]]) 1620 sage: f2 = Piecewise([[(1, 3), x^2]]) 1621 sage: f1._make_compatible(f2) 1622 (Piecewise defined function with 2 parts, [[(0, 2), x], [(2, 3), 0]], 1623 Piecewise defined function with 2 parts, [[(0, 1), 0], [(1, 3), x^2]], 1624 [(0, 1), (1, 2), (2, 3)]) 1625 """ 1626 a1, b1 = self.domain() 1627 a2, b2 = other.domain() 1628 a = min(a1, a2) 1629 b = max(b1, b2) 1630 F = self.extend_by_zero_to(a,b) 1631 G = other.extend_by_zero_to(a,b) 1632 endpts = list(set(F.end_points()).union(set(G.end_points()))) 1633 endpts.sort() 1634 return F, G, zip(endpts, endpts[1:]) 1635 1636 def __add__(self,other): 1637 """ 1638 Returns the piecewise defined function which is the sum of self and 1639 other. Does not require both domains be the same. 1640 1641 EXAMPLES:: 1642 1643 sage: x = PolynomialRing(QQ,'x').gen() 1644 sage: f1 = x^0 1645 sage: f2 = 1x 1646 sage: f3 = 2*x 1647 sage: f4 = 10x 1648 sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]]) 1649 sage: g1 = x2 1650 sage: g2 = x5 1651 sage: g = Piecewise([[(0,5),g1],[(5,10),g2]]) 1652 sage: h = f+g 1653 sage: h 1654 Piecewise defined function with 5 parts, [[(0, 1), x  1], [(1, 2), 1], [(2, 3), 3*x  2], [(3, 5), 8], [(5, 10), 5]] 1655 1656 Note that in this case the functions must be defined using 1657 polynomial expressions *not* using the lambda notation. 1658 """ 1659 F, G, intervals = self._make_compatible(other) 1660 fcn = [] 1661 for a,b in intervals: 1662 fcn.append([(a,b), F.which_function(b)+G.which_function(b)]) 1663 return Piecewise(fcn) 1664 1665 def __mul__(self,other): 1666 r""" 1667 Returns the piecewise defined function which is the product of one 1668 piecewise function (self) with another one (other). 1669 1670 EXAMPLES:: 1671 1672 sage: x = PolynomialRing(QQ,'x').gen() 1673 sage: f1 = x^0 1674 sage: f2 = 1x 1675 sage: f3 = 2*x 1676 sage: f4 = 10x 1677 sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]]) 1678 sage: g1 = x2 1679 sage: g2 = x5 1680 sage: g = Piecewise([[(0,5),g1],[(5,10),g2]]) 1681 sage: h = f*g 1682 sage: h 1683 Piecewise defined function with 5 parts, [[(0, 1), x  2], [(1, 2), x^2 + 3*x  2], [(2, 3), 2*x^2  4*x], [(3, 5), x^2 + 12*x  20], [(5, 10), x^2 + 15*x  50]] 1684 sage: g*(11/2) 1685 Piecewise defined function with 2 parts, [[(0, 5), 11/2*x  11], [(5, 10), 11/2*x  55/2]] 1686 1687 Note that in this method the functions must be defined using 1688 polynomial expressions *not* using the lambda notation. 1689 """ 1690 ## needed for scalar multiplication 1691 if isinstance(other,Rational) or isinstance(other,Integer): 1692 return Piecewise([[(a,b), other*f] for (a,b),f in self.list()]) 1693 else: 1694 F, G, intervals = self._make_compatible(other) 1695 fcn = [] 1696 for a,b in intervals: 1697 fcn.append([(a,b),F.which_function(b)*G.which_function(b)]) 1698 return Piecewise(fcn) 513 A tuple of piecewise functions, each having only a single 514 expression. 1699 515 1700 __rmul__ = __mul__516 EXAMPLES:: 1701 517 1702 def __eq__(self,other): 1703 """ 1704 Implements Boolean == operator. 518 sage: p = piecewise([((1, 0), x), ([0, 1], x)], var=x) 519 sage: p.pieces() 520 (piecewise(x>x on (1, 0); x), 521 piecewise(x>x on [0, 1]; x)) 522 """ 523 result = [] 524 for domain, func in parameters: 525 result.append(piecewise([(domain, func)], var=variable)) 526 return tuple(result) 1705 527 1706 EXAMPLES:: 1707 1708 sage: f1 = x^0 1709 sage: f2 = 1x 1710 sage: f3 = 2*x 1711 sage: f4 = 10x 1712 sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]]) 1713 sage: g = Piecewise([[(0,1),1],[(1,2),f2],[(2,3),f3],[(3,10),f4]]) 1714 sage: f==g 1715 True 1716 """ 1717 return self.list()==other.list() 528 529 530 531 piecewise = PiecewiseFunction() 532 533 def Piecewise(*args, **kwds): 534 """ 535 Deprecated spelling of ``piecewise`` 536 537 EXAMPLES:: 538 539 sage: Piecewise([((1, 0), x), ([0, 1], x)], var=x) 540 doctest:...: DeprecationWarning: use lowercase piecewise instead 541 See http://trac.sagemath.org/14801 for details. 542 piecewise(x>x on (1, 0), x>x on [0, 1]; x) 543 """ 544 from sage.misc.superseded import deprecation 545 deprecation(14801, 'use lowercase piecewise instead') 546 return piecewise(*args, **kwds) 547 548 549 
sage/symbolic/expression_conversions.py
diff git a/sage/symbolic/expression_conversions.py b/sage/symbolic/expression_conversions.py
a b 20 20 from sage.symbolic.pynac import I 21 21 from sage.functions.all import exp 22 22 from sage.symbolic.operators import arithmetic_operators, relation_operators, FDerivativeOperator 23 from sage.functions.piecewise import piecewise 23 24 from sage.rings.number_field.number_field_element_quadratic import NumberFieldElement_quadratic 24 25 GaussianField = I.pyobject().parent() 25 26 … … 1181 1182 try: 1182 1183 return self.ff.fast_float_constant(float(ex)) 1183 1184 except TypeError: 1184 raise ValueError, "free variable: %s" % repr(ex)1185 raise NotImplementedError, "free variable: %s" % repr(ex) 1185 1186 1186 1187 def arithmetic(self, ex, operator): 1187 1188 """