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)

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]

Integration:

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

Taylor series:

sage: g = maxima('f/sinh(k*x)^4')
sage: g.taylor('x', 0, 3)
f/(k^4*x^4) - 2*f/(3*k^2*x^2) + 11*f/45 - 62*k^2*f*x^2/945

Laplace transforms:

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)

2.2 Ordinary differential equations

If you have Octave 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.
$

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

You can also use some of the programs included in the subdirectory examples/calculus but you must ``load'' or ``attach'' then yourself:

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$ .

See About this document... for information on suggesting changes.