Up One Level Constructions in SAGE (How to) Contents Index
Previous: Up: Constructions in SAGE (How Next:
Subsections

2. Calculus

Here are some examples of calculus symbolic computations using SAGE. They use the Maxima interface.

2.1 Using maxima

Differentiation: 

sage: f = maxima('x^3 * %e^(k*x) * sin(w*x)'); f
x^3*%e^(k*x)*sin(w*x)
sage: f.diff('x')
k*x^3*%e^(k*x)*sin(w*x) + 3*x^2*%e^(k*x)*sin(w*x) + w*x^3*%e^(k*x)*cos(w*x)
sage: latex(f.diff('x'))
 k x^3 e^{k x} \sin \left(w x\right)+3 x^2 e^{k x} \sin \left(w x  \right)+w x^3 e^{k x} \cos \left(w x\right)
If you type view(f.diff('x')) another window will open up displaying the compiled latex output. In the SAGE notebook, you can enter 
f = maxima('x^3 * %e^(k*x) * sin(w*x)')
show(f)
show(f.diff('x'))
into a cell and press shift-enter for a similar result. You can also call Maxima indirectly using the commands 
R = PolynomialRing(QQ,"x")
x = R.gen()
p = x^2 + 1
show(p.derivative())
show(p.integral())
in a notebook cell, or 
sage: R = PolynomialRing(QQ,"x")
sage: x = R.gen()
sage: p = x^2 + 1
sage: view(p.derivative())
sage: view(p.integral())
on the command line.

Definite integrals involving limits is also possible: 

sage: f = maxima('x^3 * %e^(k*x) * sin(w*x)')
sage: f.integrate('x')
(((k*w^6 + 3*k^3*w^4 + 3*k^5*w^2 + k^7)*x^3 + (3*w^6 + 3*k^2*w^4 -
3*k^4*w^2 - 3*k^6)*x^2 + ( - 18*k*w^4 - 12*k^3*w^2 + 6*k^5)*x - 6*w^4 +
36*k^2*w^2 - 6*k^4)*%e^(k*x)*sin(w*x) + (( - w^7 - 3*k^2*w^5 - 3*k^4*w^3 -
k^6*w)*x^3 + (6*k*w^5 + 12*k^3*w^3 + 6*k^5*w)*x^2 + (6*w^5 - 12*k^2*w^3 -
18*k^4*w)*x - 24*k*w^3 + 24*k^3*w)*%e^(k*x)*cos(w*x))/(w^8 + 4*k^2*w^6 +
6*k^4*w^4 + 4*k^6*w^2 + k^8)
sage: f = maxima('1/x^2')
sage: f.integrate('x', 1, 'inf')
1

You can find the convolution of any piecewise defined function with another (off the domain of definition, they are assumed to be zero). Here is $ f$ , $ f*f$ , and $ f*f*f$ , where $ f(x)=1$ , $ 0<x<1$ : 

sage: x = PolynomialRing(QQ).gen()
sage: f = Piecewise([[(0,1),1*x^0]])
sage: g = f.convolution(f)
sage: h = f.convolution(g)
sage: P = f.plot(); Q = g.plot(rgbcolor=(1,1,0)); R = h.plot(rgbcolor=(0,1,1))
To view this, type show(PQ+R)+.

You can find critical points of a piecewise defined function: 

sage: x = PolynomialRing(RationalField()).gen()
sage: f1 = x^0
sage: f2 = 1-x
sage: f3 = 2*x
sage: f4 = 10*x-x^2
sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]])
sage: f.critical_points()
[5.0]

Taylor series: 

sage: g = maxima('f0/sinh(k*x)^4')
sage: g.taylor('x', 0, 3)
f0/(k^4*x^4) - 2*f0/(3*k^2*x^2) + 11*f0/45 - 62*k^2*f0*x^2/945
sage: g.powerseries('x',0)
 16*f0*('sum((2^(2*i1 - 1) - 1)*bern(2*i1)*k^(2*i1 - 1)*x^(2*i1 - 1)/(2*i1)!,i1,0,inf))^4
Of course, you can view the latexed version of this using view(g.powerseries('x',0)).

The Maclaurin and power series of $ \log({\frac{\sin(x)}{x}})$ : 

sage: f = maxima("log(sin(x)/x)")
sage: f.powerseries("x",0)
 ('sum(( - 1)^i2*2^(2*i2)*bern(2*i2)*x^(2*i2)/(i2*(2*i2)!),i2,1,inf))/2
sage: maxima.eval("taylor (log(sin(x)/x), x, 0, 10);")
 '-x^2/6 - x^4/180 - x^6/2835 - x^8/37800 - x^10/467775'
sage: [bernoulli(2*i) for i in range(1,7)]
 [1/6, -1/30, 1/42, -1/30, 5/66, -691/2730]

Laplace transforms: 

sage: f = maxima('1/exp(k*t)')
sage: f.laplace("t","s")
 1/(s + k)
is one way to compute LT's and 
sage: _ = maxima.eval("f(t) := t^5*exp(t)*sin(t)")
sage: maxima("laplace(f(t),t,s)")
360*(2*s - 2)/(s^2 - 2*s + 2)^4 - 480*(2*s - 2)^3/(s^2 - 2*s + 2)^5 +
120*(2*s - 2)^5/(s^2 - 2*s + 2)^6
sage: maxima("laplace(f(t),t,s)").display2d()
                                         3                 5
           360 (2 s - 2)    480 (2 s - 2)     120 (2 s - 2)
          --------------- - --------------- + ---------------
            2           4     2           5     2           6
          (s  - 2 s + 2)    (s  - 2 s + 2)    (s  - 2 s + 2)
is another way.

2.2 Ordinary differential equations

If you have Octave and gnuplot installed, 

sage.: octave.de_system_plot(['x+y','x-y'], [1,-1], [0,2])
yields the two plots $ (t,x(t)), (t,y(t))$ on the same graph (the $ t$ -axis is the horizonal axis) of the system of ODEs

$\displaystyle x' = x+y, x(0) = 1;\qquad y' = x-y, y(0) = -1,$   for$\displaystyle \quad 0 < t < 2. $

Another way this system can be solved is to use the command desolve_system in calculus/examples. 

sage: attach os.environ['SAGE_ROOT'] + '/examples/calculus/desolvers.sage'
sage: des = ["'diff(x(t),t)=x(t)+y(t)","'diff(y(t),t)=x(t)-y(t)"]
sage: vars = ["t","x","y"]
sage: ics = [0,1,-1]
sage: desolve_system(des,vars,ics)
 '[x(t) = cosh(sqrt(2)*t),y(t) = 2*sinh(sqrt(2)*t)/sqrt(2) - cosh(sqrt(2)*t)]'
This outputs a string, not a pair of SAGE functions.

You can solve ODE's symbolically in SAGE using the Maxima interface: 

sage.: maxima.de_solve('diff(y,x,2) + 3*x = y', ['x','y'], [1,1,1])
    y = 3*x - 2*%e^(x - 1)
sage.: maxima.de_solve('diff(y,x,2) + 3*x = y', ['x','y'])
    y = %k1*%e^x + %k2*%e^ - x + 3*x
sage.: maxima.de_solve('diff(y,x) + 3*x = y', ['x','y'])
    y = (%c - 3*( - x - 1)*%e^ - x)*%e^x
sage.: maxima.de_solve('diff(y,x) + 3*x = y', ['x','y'],[1,1])
    y =  - %e^ - 1*(5*%e^x - 3*%e*x - 3*%e)

sage.: maxima.clear('x'); maxima.clear('f')
sage.: maxima.de_solve_laplace("diff(f(x),x,2) = 2*diff(f(x),x)-f(x)", ["x","f"], [0,1,2])
   f(x) = x*%e^x + %e^x
            
sage.: maxima.clear('x'); maxima.clear('f')            
sage.: f = maxima.de_solve_laplace("diff(f(x),x,2) = 2*diff(f(x),x)-f(x)", ["x","f"])
sage.: f
   f(x) = x*%e^x*(at('diff(f(x),x,1),x = 0)) - f(0)*x*%e^x + f(0)*%e^x
sage.: f.display2d()
                                               !
                                   x  d        !                  x          x
                        f(x) = x %e  (-- (f(x))!     ) - f(0) x %e  + f(0) %e
                                      dx       !
                                               !x = 0

Finally, SAGE can solve linear DEs using power series: 

sage: R.<t> = PowerSeriesRing(QQ, default_prec=10)
sage: a = 2 - 3*t + 4*t^2 + O(t^10)
sage: b = 3 - 4*t^2 + O(t^7)
sage: f = a.solve_linear_de(prec=5, b=b, f0=3/5)
sage: f
 3/5 + 21/5*t + 33/10*t^2 - 38/15*t^3 + 11/24*t^4 + O(t^5)
sage: f.derivative() - a*f - b
 O(t^4)

2.3 Euler's method for numerically solving an ODE

The goal is to find an approximate solution to the problem

$\displaystyle y'=f(x,y), y(a)=c,$ (2.1)

where $ f(x,y)$ is some given function. We shall try to approximate the value of the solution at $ x=b$ , where $ b>a$ is given.

Tabular idea: Let $ n>0$ be an integer, which we call the step size. This is related to the increment by

$\displaystyle h=\frac{b-a}{n}. $

This can be expressed simplest using a table.

$ x$ $ y$ $ hf(x,y)$
$ a$ $ c$ $ hf(a,c)$
$ a+h$ $ c+hf(a,c)$ &vellip#vdots;
$ a+2h$ &vellip#vdots;  
&vellip#vdots;    
$ b$ ??? xxx

The goal is to fill out all the blanks of the table but the xxx entry and find the ??? entry, which is the Euler's method approximation for $ y(b)$ . This is implemented in eulers_method.sage.

To find an approximation to $ y(1)$ using two steps of Euler's method, where $ y'=5x+y-5$ , $ y(0)=1$ , using SAGE, type: 

 
        sage: attach "eulers_method.sage"
        sage: x,y=PolynomialRing(QQ,2).gens()
        sage: eulers_method(5*x+y-5,0,1,1/2,1)  
         x                    y                  h*f(x,y)
         0                    1                   -2
       1/2                   -1                 -7/4
         1                -11/4                -11/8
        sage: RR = RealField(sci_not=0, prec=4, rnd='RNDU')
        sage: x,y=PolynomialRing(RR,2).gens()
        sage: eulers_method(5*x+y-5,0,1,1/2,1)
         x                    y                  h*f(x,y)
         0                    1                -2.00
       1/2                -1.00                -1.75
         1                -2.75                -1.37

So $ y(1)\approx -11/4 = -2.75$ .

2.4 Fourier series of periodic functions

If $ f(x)$ is a piecewise-defined polynomial function on $ -L<x<L$ then the Fourier series

$\displaystyle 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})], $

converges. In addition to computing the coefficients $ a_n,b_n$ , it will also compute the partial sums (as a string), plot the partial sums (as a function of $ x$ over $ (-L,L)$ , for comparison with the plot of $ f(x)$ itself), compute the value of the FS at a point, and similar computations for the cosine series (if $ f(x)$ is even) and the sine series (if $ f(x)$ is odd). Also, it will plot the partial F.S. Cesaro mean sums (a ``smoother" partial sum illustrating how the Gibbs phenomenon is mollified). 

 
sage: f1 = lambda x:-1
sage: f2 = lambda x:2
sage: f = Piecewise([[(0,pi/2),f1],[(pi/2,pi),f2]])
sage: f.fourier_series_cosine_coefficient(5,pi)
(-3/(5*pi))
sage: f.fourier_series_sine_coefficient(2,pi)
 (-3/pi)
sage: f.fourier_series_partial_sum(3,pi)
 '1/4 + ((-3/pi)*cos(1*pi*x/pi) + (1/pi)*sin(1*pi*x/pi)) + 
  (0*cos(2*pi*x/pi) + (-3/pi)*sin(2*pi*x/pi))'
Type show(f.plot_fourier_series_partial_sum(15,pi,-5,5)) and show(f.plot_fourier_series_partial_sum_cesaro(15,pi,-5,5)) (and be patient) to view the partial sums.

2.5 Riemann and trapezoid sums for integrals

Regarding numerical approximation of $ \int_a^bf(x) dx$ , where $ f$ is a piecewise defined function, SAGE can 

 
sage: f1 = lambda x:x^2      
sage: f2 = lambda x:5-x^2
sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])
sage: f.trapezoid(4)
 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]]
sage: f.riemann_sum_integral_approximation(6,mode="right")
 19/6
sage: f.integral()
 3

Up One Level Constructions in SAGE (How to) Contents Index
Previous: Up: Constructions in SAGE (How Next:
Release 1.4.1, documentation updated on October 18, 2006.
See About this document... for information on suggesting changes.