Calculus

Here are some examples of calculus symbolic computations using Sage.

Assumptions

By default a symbolic variable is as general as possible: it may stand for a positive, negative, or complex value. Many calculus results hold only under extra hypotheses, and Sage will not apply them unless you say so. Use assume() to attach a hypothesis to a variable, assumptions() to list what is currently assumed, and forget() to remove assumptions again. For example, \(\sqrt{x^2} = x\) holds only for positive \(x\):

sage: assume(x > 0)
sage: bool(sqrt(x^2) == x)
True
sage: assumptions()
[x > 0]
sage: forget()
sage: bool(sqrt(x^2) == x)
False

>>> from sage.all import *
>>> assume(x > Integer(0))
>>> bool(sqrt(x**Integer(2)) == x)
True
>>> assumptions()
[x > 0]
>>> forget()
>>> bool(sqrt(x**Integer(2)) == x)
False

assume(x > 0)
bool(sqrt(x^2) == x)
assumptions()
forget()
bool(sqrt(x^2) == x)

Assumptions decide integrals and sums whose answers depend on a sign or range condition:

sage: n = var('n')
sage: assume(n + 1 > 0)
sage: integral(x^n, x)
x^(n + 1)/(n + 1)
sage: forget()

>>> from sage.all import *
>>> n = var('n')
>>> assume(n + Integer(1) > Integer(0))
>>> integral(x**n, x)
x^(n + 1)/(n + 1)
>>> forget()

n = var('n')
assume(n + 1 > 0)
integral(x^n, x)
forget()

sage: var('q, a, k')
(q, a, k)
sage: assume(abs(q) < 1)
sage: sum(a*q^k, k, 0, oo)
-a/(q - 1)
sage: forget()

>>> from sage.all import *
>>> var('q, a, k')
(q, a, k)
>>> assume(abs(q) < Integer(1))
>>> sum(a*q**k, k, Integer(0), oo)
-a/(q - 1)
>>> forget()

var('q, a, k')
assume(abs(q) < 1)
sum(a*q^k, k, 0, oo)
forget()

A variable can also be declared to have a property such as 'integer', 'real', 'positive', or 'even':

sage: n = var('n')
sage: assume(n, 'integer')
sage: sin(n*pi)
0
sage: forget()

>>> from sage.all import *
>>> n = var('n')
>>> assume(n, 'integer')
>>> sin(n*pi)
0
>>> forget()

n = var('n')
assume(n, 'integer')
sin(n*pi)
forget()

Pass each relationship as a separate assumption. A chained comparison is not an error, but it is silently truncated to its first relation (Python evaluates 0 < x < 1 as (0 < x) and (x < 1), and the undecided first comparison short-circuits the and), so only 0 < x gets recorded:

sage: assume(0 < x < 1)
sage: assumptions()
[0 < x]
sage: forget()
sage: assume(0 < x, x < 1)
sage: assumptions()
[0 < x, x < 1]
sage: forget()

>>> from sage.all import *
>>> assume(Integer(0) < x < Integer(1))
>>> assumptions()
[0 < x]
>>> forget()
>>> assume(Integer(0) < x, x < Integer(1))
>>> assumptions()
[0 < x, x < 1]
>>> forget()

assume(0 < x < 1)
assumptions()
forget()
assume(0 < x, x < 1)
assumptions()
forget()

For a temporary hypothesis, use assuming in a with block: the assumption is in force inside the block and automatically forgotten afterwards:

sage: solve(x^2 == 4, x)
[x == -2, x == 2]
sage: with assuming(x > 0):
....:     solve(x^2 == 4, x)
[x == 2]
sage: assumptions()
[]

>>> from sage.all import *
>>> solve(x**Integer(2) == Integer(4), x)
[x == -2, x == 2]
>>> with assuming(x > Integer(0)):
...     solve(x**Integer(2) == Integer(4), x)
[x == 2]
>>> assumptions()
[]

solve(x^2 == 4, x)
with assuming(x > 0):
    solve(x^2 == 4, x)
assumptions()

A restriction can also be built into a variable when it is created, with the domain keyword (this too registers an assumption):

sage: t = var('t', domain='positive')
sage: abs(t)
t
sage: forget()

>>> from sage.all import *
>>> t = var('t', domain='positive')
>>> abs(t)
t
>>> forget()

t = var('t', domain='positive')
abs(t)
forget()

If you have used Maple’s assume/about or Mathematica’s Assuming, this is the corresponding mechanism in Sage; see sage.symbolic.assumptions for the full reference, including forget() of individual assumptions.

Differentiation

Differentiation:

sage: var('x k w')
(x, k, w)
sage: f = x^3 * e^(k*x) * sin(w*x); f
x^3*e^(k*x)*sin(w*x)
sage: f.diff(x)
w*x^3*cos(w*x)*e^(k*x) + k*x^3*e^(k*x)*sin(w*x) + 3*x^2*e^(k*x)*sin(w*x)
sage: latex(f.diff(x))
w x^{3} \cos\left(w x\right) e^{\left(k x\right)} + k x^{3} e^{\left(k x\right)} \sin\left(w x\right) + 3 \, x^{2} e^{\left(k x\right)} \sin\left(w x\right)

>>> from sage.all import *
>>> var('x k w')
(x, k, w)
>>> f = x**Integer(3) * e**(k*x) * sin(w*x); f
x^3*e^(k*x)*sin(w*x)
>>> f.diff(x)
w*x^3*cos(w*x)*e^(k*x) + k*x^3*e^(k*x)*sin(w*x) + 3*x^2*e^(k*x)*sin(w*x)
>>> latex(f.diff(x))
w x^{3} \cos\left(w x\right) e^{\left(k x\right)} + k x^{3} e^{\left(k x\right)} \sin\left(w x\right) + 3 \, x^{2} e^{\left(k x\right)} \sin\left(w x\right)

var('x k w')
f = x^3 * e^(k*x) * sin(w*x); f
f.diff(x)
latex(f.diff(x))

If you type view(f.diff(x)) another window will open up displaying the compiled output. In the notebook, you can enter

var('x k w')
f = 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 differentiate and integrate 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: p.derivative()
2*x
sage: p.integral()
1/3*x^3 + x

>>> from sage.all import *
>>> R = PolynomialRing(QQ,"x")
>>> x = R.gen()
>>> p = x**Integer(2) + Integer(1)
>>> p.derivative()
2*x
>>> p.integral()
1/3*x^3 + x

R = PolynomialRing(QQ,"x")
x = R.gen()
p = x^2 + 1
p.derivative()
p.integral()

on the command line. At this point you can also type view(p.derivative()) or view(p.integral()) to open a new window with output typeset by LaTeX.

Critical points

You can find critical points of a piecewise defined function:

sage: x = PolynomialRing(RationalField(), 'x').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]

>>> from sage.all import *
>>> x = PolynomialRing(RationalField(), 'x').gen()
>>> f1 = x**Integer(0)
>>> f2 = Integer(1)-x
>>> f3 = Integer(2)*x
>>> f4 = Integer(10)*x-x**Integer(2)
>>> f = piecewise([((Integer(0),Integer(1)),f1), ((Integer(1),Integer(2)),f2), ((Integer(2),Integer(3)),f3), ((Integer(3),Integer(10)),f4)])
>>> f.critical_points()
[5.0]

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

Power series

Sage offers several ways to construct and work with power series.

To get Taylor series from function expressions use the method .taylor() on the expression:

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

>>> from sage.all import *
>>> var('f0 k x')
(f0, k, x)
>>> g = f0/sinh(k*x)**Integer(4)
>>> g.taylor(x, Integer(0), Integer(3))
-62/945*f0*k^2*x^2 + 11/45*f0 - 2/3*f0/(k^2*x^2) + f0/(k^4*x^4)

var('f0 k x')
g = f0/sinh(k*x)^4
g.taylor(x, 0, 3)

Formal power series expansions of functions can be had with the .series() method:

sage: (1/(2-cos(x))).series(x,7)
1 + (-1/2)*x^2 + 7/24*x^4 + (-121/720)*x^6 + Order(x^7)

>>> from sage.all import *
>>> (Integer(1)/(Integer(2)-cos(x))).series(x,Integer(7))
1 + (-1/2)*x^2 + 7/24*x^4 + (-121/720)*x^6 + Order(x^7)

(1/(2-cos(x))).series(x,7)

Certain manipulations on such series are hard to perform at the moment, however. There are two alternatives: either use the Maxima subsystem of Sage for full symbolic functionality:

sage: f = log(sin(x)/x)
sage: f.taylor(x, 0, 10)
-1/467775*x^10 - 1/37800*x^8 - 1/2835*x^6 - 1/180*x^4 - 1/6*x^2
sage: maxima(f).powerseries(x,0)._sage_()
sum(2^(2*i... - 1)*(-1)^i...*x^(2*i...)*bern(2*i...)/(i...*factorial(2*i...)), i..., 1, +Infinity)

>>> from sage.all import *
>>> f = log(sin(x)/x)
>>> f.taylor(x, Integer(0), Integer(10))
-1/467775*x^10 - 1/37800*x^8 - 1/2835*x^6 - 1/180*x^4 - 1/6*x^2
>>> maxima(f).powerseries(x,Integer(0))._sage_()
sum(2^(2*i... - 1)*(-1)^i...*x^(2*i...)*bern(2*i...)/(i...*factorial(2*i...)), i..., 1, +Infinity)

f = log(sin(x)/x)
f.taylor(x, 0, 10)
maxima(f).powerseries(x,0)._sage_()

Or you can use the formal power series rings for fast computation. These are missing symbolic functions, on the other hand:

sage: R.<w> = QQ[[]]
sage: ps = w + 17/2*w^2 + 15/4*w^4 + O(w^6); ps
w + 17/2*w^2 + 15/4*w^4 + O(w^6)
sage: ps.exp()
1 + w + 9*w^2 + 26/3*w^3 + 265/6*w^4 + 413/10*w^5 + O(w^6)
sage: (1+ps).log()
w + 8*w^2 - 49/6*w^3 - 193/8*w^4 + 301/5*w^5 + O(w^6)
sage: (ps^1000).coefficients()
[1, 8500, 36088875, 102047312625, 1729600092867375/8]

>>> from sage.all import *
>>> R = QQ[['w']]; (w,) = R._first_ngens(1)
>>> ps = w + Integer(17)/Integer(2)*w**Integer(2) + Integer(15)/Integer(4)*w**Integer(4) + O(w**Integer(6)); ps
w + 17/2*w^2 + 15/4*w^4 + O(w^6)
>>> ps.exp()
1 + w + 9*w^2 + 26/3*w^3 + 265/6*w^4 + 413/10*w^5 + O(w^6)
>>> (Integer(1)+ps).log()
w + 8*w^2 - 49/6*w^3 - 193/8*w^4 + 301/5*w^5 + O(w^6)
>>> (ps**Integer(1000)).coefficients()
[1, 8500, 36088875, 102047312625, 1729600092867375/8]

R.<w> = QQ[[]]
ps = w + 17/2*w^2 + 15/4*w^4 + O(w^6); ps
ps.exp()
(1+ps).log()
(ps^1000).coefficients()

Integration

Numerical integration is discussed in Riemann and trapezoid sums for integrals below.

Sage can integrate some simple functions on its own:

sage: f = x^3
sage: f.integral(x)
1/4*x^4
sage: integral(x^3,x)
1/4*x^4
sage: f = x*sin(x^2)
sage: integral(f,x)
-1/2*cos(x^2)

>>> from sage.all import *
>>> f = x**Integer(3)
>>> f.integral(x)
1/4*x^4
>>> integral(x**Integer(3),x)
1/4*x^4
>>> f = x*sin(x**Integer(2))
>>> integral(f,x)
-1/2*cos(x^2)

f = x^3
f.integral(x)
integral(x^3,x)
f = x*sin(x^2)
integral(f,x)

Sage can also compute symbolic definite integrals involving limits.

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

>>> from sage.all import *
>>> var('x, k, w')
(x, k, w)
>>> f = x**Integer(3) * e**(k*x) * sin(w*x)
>>> f.integrate(x)
((24*k^3*w - 24*k*w^3 - (k^6*w + 3*k^4*w^3 + 3*k^2*w^5 + w^7)*x^3 + 6*(k^5*w + 2*k^3*w^3 + k*w^5)*x^2 - 6*(3*k^4*w + 2*k^2*w^3 - w^5)*x)*cos(w*x)*e^(k*x) - (6*k^4 - 36*k^2*w^2 + 6*w^4 - (k^7 + 3*k^5*w^2 + 3*k^3*w^4 + k*w^6)*x^3 + 3*(k^6 + k^4*w^2 - k^2*w^4 - w^6)*x^2 - 6*(k^5 - 2*k^3*w^2 - 3*k*w^4)*x)*e^(k*x)*sin(w*x))/(k^8 + 4*k^6*w^2 + 6*k^4*w^4 + 4*k^2*w^6 + w^8)
>>> integrate(Integer(1)/x**Integer(2), x, Integer(1), infinity)
1

var('x, k, w')
f = x^3 * e^(k*x) * sin(w*x)
f.integrate(x)
integrate(1/x^2, x, 1, infinity)

Convolution

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, 'x').gen()
sage: f = piecewise([((0,1),1*x^0)])
sage: g = f.convolution(f)
sage: h = f.convolution(g)
sage: set_verbose(-1)
sage: P = f.plot(); Q = g.plot(rgbcolor=(1,1,0)); R = h.plot(rgbcolor=(0,1,1))

>>> from sage.all import *
>>> x = PolynomialRing(QQ, 'x').gen()
>>> f = piecewise([((Integer(0),Integer(1)),Integer(1)*x**Integer(0))])
>>> g = f.convolution(f)
>>> h = f.convolution(g)
>>> set_verbose(-Integer(1))
>>> P = f.plot(); Q = g.plot(rgbcolor=(Integer(1),Integer(1),Integer(0))); R = h.plot(rgbcolor=(Integer(0),Integer(1),Integer(1)))

x = PolynomialRing(QQ, 'x').gen()
f = piecewise([((0,1),1*x^0)])
g = f.convolution(f)
h = f.convolution(g)
set_verbose(-1)
P = f.plot(); Q = g.plot(rgbcolor=(1,1,0)); R = h.plot(rgbcolor=(0,1,1))

To view this, type show(P+Q+R).

Riemann and trapezoid sums for integrals

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

  • compute (for plotting purposes) the piecewise linear function defined by the trapezoid rule for numerical integration based on a subdivision into \(N\) subintervals

  • the approximation given by the trapezoid rule,

  • compute (for plotting purposes) the piecewise constant function defined by the Riemann sums (left-hand, right-hand, or midpoint) in numerical integration based on a subdivision into \(N\) subintervals,

  • the approximation given by the Riemann sum approximation.

sage: f1(x) = x^2
sage: f2(x) = 5-x^2
sage: f = piecewise([[[0,1], f1], [RealSet.open_closed(1,2), f2]])
sage: t = f.trapezoid(2); t
piecewise(x|-->1/2*x on (0, 1/2), x|-->3/2*x - 1/2 on (1/2, 1), x|-->7/2*x - 5/2 on (1, 3/2), x|-->-7/2*x + 8 on (3/2, 2); x)
sage: t.integral()
piecewise(x|-->1/4*x^2 on (0, 1/2), x|-->3/4*x^2 - 1/2*x + 1/8 on (1/2, 1), x|-->7/4*x^2 - 5/2*x + 9/8 on (1, 3/2), x|-->-7/4*x^2 + 8*x - 27/4 on (3/2, 2); x)
sage: t.integral(definite=True)
9/4

>>> from sage.all import *
>>> __tmp__=var("x"); f1 = symbolic_expression(x**Integer(2)).function(x)
>>> __tmp__=var("x"); f2 = symbolic_expression(Integer(5)-x**Integer(2)).function(x)
>>> f = piecewise([[[Integer(0),Integer(1)], f1], [RealSet.open_closed(Integer(1),Integer(2)), f2]])
>>> t = f.trapezoid(Integer(2)); t
piecewise(x|-->1/2*x on (0, 1/2), x|-->3/2*x - 1/2 on (1/2, 1), x|-->7/2*x - 5/2 on (1, 3/2), x|-->-7/2*x + 8 on (3/2, 2); x)
>>> t.integral()
piecewise(x|-->1/4*x^2 on (0, 1/2), x|-->3/4*x^2 - 1/2*x + 1/8 on (1/2, 1), x|-->7/4*x^2 - 5/2*x + 9/8 on (1, 3/2), x|-->-7/4*x^2 + 8*x - 27/4 on (3/2, 2); x)
>>> t.integral(definite=True)
9/4

f1(x) = x^2
f2(x) = 5-x^2
f = piecewise([[[0,1], f1], [RealSet.open_closed(1,2), f2]])
t = f.trapezoid(2); t
t.integral()
t.integral(definite=True)

Laplace transforms

If you have a piecewise-defined polynomial function then there is a “native” command for computing Laplace transforms. This calls Maxima but it’s worth noting that Maxima cannot handle (using the direct interface illustrated in the last few examples) this type of computation.

sage: var('x s')
(x, s)
sage: f1(x) = 1
sage: f2(x) = 1-x
sage: f = piecewise([((0,1),f1), ((1,2),f2)])
sage: f.laplace(x, s)
-e^(-s)/s + (s + 1)*e^(-2*s)/s^2 + 1/s - e^(-s)/s^2

>>> from sage.all import *
>>> var('x s')
(x, s)
>>> __tmp__=var("x"); f1 = symbolic_expression(Integer(1)).function(x)
>>> __tmp__=var("x"); f2 = symbolic_expression(Integer(1)-x).function(x)
>>> f = piecewise([((Integer(0),Integer(1)),f1), ((Integer(1),Integer(2)),f2)])
>>> f.laplace(x, s)
-e^(-s)/s + (s + 1)*e^(-2*s)/s^2 + 1/s - e^(-s)/s^2

var('x s')
f1(x) = 1
f2(x) = 1-x
f = piecewise([((0,1),f1), ((1,2),f2)])
f.laplace(x, s)

For other “reasonable” functions, Laplace transforms can be computed using the Maxima interface:

sage: var('k, s, t')
(k, s, t)
sage: f = 1/exp(k*t)
sage: f.laplace(t,s)
1/(k + s)

>>> from sage.all import *
>>> var('k, s, t')
(k, s, t)
>>> f = Integer(1)/exp(k*t)
>>> f.laplace(t,s)
1/(k + s)

var('k, s, t')
f = 1/exp(k*t)
f.laplace(t,s)

is one way to compute LT’s and

sage: var('s, t')
(s, t)
sage: f = t^5*exp(t)*sin(t)
sage: L = laplace(f, t, s); L
3840*(s - 1)^5/(s^2 - 2*s + 2)^6 - 3840*(s - 1)^3/(s^2 - 2*s + 2)^5 +
720*(s - 1)/(s^2 - 2*s + 2)^4

>>> from sage.all import *
>>> var('s, t')
(s, t)
>>> f = t**Integer(5)*exp(t)*sin(t)
>>> L = laplace(f, t, s); L
3840*(s - 1)^5/(s^2 - 2*s + 2)^6 - 3840*(s - 1)^3/(s^2 - 2*s + 2)^5 +
720*(s - 1)/(s^2 - 2*s + 2)^4

var('s, t')
f = t^5*exp(t)*sin(t)
L = laplace(f, t, s); L

is another way.

Ordinary differential equations

Symbolically solving ODEs can be done using Sage interface with Maxima. See

sage:desolvers?

>>> from sage.all import *
sage:desolvers?

for available commands. Numerical solution of ODEs can be done using Sage interface with Octave (an experimental package), or routines in the GSL (Gnu Scientific Library).

An example, how to solve ODE’s symbolically in Sage using the Maxima interface (do not type the ....:):

sage: y=function('y')(x); desolve(diff(y,x,2) + 3*x == y, dvar = y, ics = [1,1,1])
3*x - 2*e^(x - 1)
sage: desolve(diff(y,x,2) + 3*x == y, dvar = y)
_K2*e^(-x) + _K1*e^x + 3*x
sage: desolve(diff(y,x) + 3*x == y, dvar = y)
(3*(x + 1)*e^(-x) + _C)*e^x
sage: desolve(diff(y,x) + 3*x == y, dvar = y, ics = [1,1]).expand()
3*x - 5*e^(x - 1) + 3

sage: f=function('f')(x); desolve_laplace(diff(f,x,2) == 2*diff(f,x)-f, dvar = f, ics = [0,1,2])
x*e^x + e^x

sage: desolve_laplace(diff(f,x,2) == 2*diff(f,x)-f, dvar = f)
-x*e^x*f(0) + x*e^x*D[0](f)(0) + e^x*f(0)

>>> from sage.all import *
>>> y=function('y')(x); desolve(diff(y,x,Integer(2)) + Integer(3)*x == y, dvar = y, ics = [Integer(1),Integer(1),Integer(1)])
3*x - 2*e^(x - 1)
>>> desolve(diff(y,x,Integer(2)) + Integer(3)*x == y, dvar = y)
_K2*e^(-x) + _K1*e^x + 3*x
>>> desolve(diff(y,x) + Integer(3)*x == y, dvar = y)
(3*(x + 1)*e^(-x) + _C)*e^x
>>> desolve(diff(y,x) + Integer(3)*x == y, dvar = y, ics = [Integer(1),Integer(1)]).expand()
3*x - 5*e^(x - 1) + 3

>>> f=function('f')(x); desolve_laplace(diff(f,x,Integer(2)) == Integer(2)*diff(f,x)-f, dvar = f, ics = [Integer(0),Integer(1),Integer(2)])
x*e^x + e^x

>>> desolve_laplace(diff(f,x,Integer(2)) == Integer(2)*diff(f,x)-f, dvar = f)
-x*e^x*f(0) + x*e^x*D[0](f)(0) + e^x*f(0)

y=function('y')(x); desolve(diff(y,x,2) + 3*x == y, dvar = y, ics = [1,1,1])
desolve(diff(y,x,2) + 3*x == y, dvar = y)
desolve(diff(y,x) + 3*x == y, dvar = y)
desolve(diff(y,x) + 3*x == y, dvar = y, ics = [1,1]).expand()
f=function('f')(x); desolve_laplace(diff(f,x,2) == 2*diff(f,x)-f, dvar = f, ics = [0,1,2])
desolve_laplace(diff(f,x,2) == 2*diff(f,x)-f, dvar = f)

If you have Octave and gnuplot installed,

sage: octave.de_system_plot(['x+y','x-y'], [1,-1], [0,2]) # optional - octave

>>> from sage.all import *
>>> octave.de_system_plot(['x+y','x-y'], [Integer(1),-Integer(1)], [Integer(0),Integer(2)]) # optional - octave

octave.de_system_plot(['x+y','x-y'], [1,-1], [0,2]) # optional - octave

yields the two plots \((t,x(t)), (t,y(t))\) on the same graph (the \(t\)-axis is the horizontal axis) of the system of ODEs

\[x' = x+y, x(0) = 1; y' = x-y, y(0) = -1,\]

for \(0 <= t <= 2\). The same result can be obtained by using desolve_system_rk4:

sage: x, y, t = var('x y t')
sage: P=desolve_system_rk4([x+y, x-y], [x,y], ics=[0,1,-1], ivar=t, end_points=2)
sage: p1 = list_plot([[i,j] for i,j,k in P], plotjoined=True)
sage: p2 = list_plot([[i,k] for i,j,k in P], plotjoined=True, color='red')
sage: p1+p2
Graphics object consisting of 2 graphics primitives

>>> from sage.all import *
>>> x, y, t = var('x y t')
>>> P=desolve_system_rk4([x+y, x-y], [x,y], ics=[Integer(0),Integer(1),-Integer(1)], ivar=t, end_points=Integer(2))
>>> p1 = list_plot([[i,j] for i,j,k in P], plotjoined=True)
>>> p2 = list_plot([[i,k] for i,j,k in P], plotjoined=True, color='red')
>>> p1+p2
Graphics object consisting of 2 graphics primitives

x, y, t = var('x y t')
P=desolve_system_rk4([x+y, x-y], [x,y], ics=[0,1,-1], ivar=t, end_points=2)
p1 = list_plot([[i,j] for i,j,k in P], plotjoined=True)
p2 = list_plot([[i,k] for i,j,k in P], plotjoined=True, color='red')
p1+p2

Another way this system can be solved is to use the command desolve_system.

sage: t=var('t'); x=function('x',t); y=function('y',t)
sage: des = [diff(x,t) == x+y, diff(y,t) == x-y]
sage: desolve_system(des, [x,y], ics = [0, 1, -1])
[x(t) == cosh(sqrt(2)*t), y(t) == sqrt(2)*sinh(sqrt(2)*t) - cosh(sqrt(2)*t)]

>>> from sage.all import *
>>> t=var('t'); x=function('x',t); y=function('y',t)
>>> des = [diff(x,t) == x+y, diff(y,t) == x-y]
>>> desolve_system(des, [x,y], ics = [Integer(0), Integer(1), -Integer(1)])
[x(t) == cosh(sqrt(2)*t), y(t) == sqrt(2)*sinh(sqrt(2)*t) - cosh(sqrt(2)*t)]

t=var('t'); x=function('x',t); y=function('y',t)
des = [diff(x,t) == x+y, diff(y,t) == x-y]
desolve_system(des, [x,y], ics = [0, 1, -1])

The output of this command is not a pair of functions.

Finally, 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)

>>> from sage.all import *
>>> R = PowerSeriesRing(QQ, default_prec=Integer(10), names=('t',)); (t,) = R._first_ngens(1)
>>> a = Integer(2) - Integer(3)*t + Integer(4)*t**Integer(2) + O(t**Integer(10))
>>> b = Integer(3) - Integer(4)*t**Integer(2) + O(t**Integer(7))
>>> f = a.solve_linear_de(prec=Integer(5), b=b, f0=Integer(3)/Integer(5))
>>> f
3/5 + 21/5*t + 33/10*t^2 - 38/15*t^3 + 11/24*t^4 + O(t^5)
>>> f.derivative() - a*f - b
O(t^4)

R.<t> = PowerSeriesRing(QQ, default_prec=10)
a = 2 - 3*t + 4*t^2 + O(t^10)
b = 3 - 4*t^2 + O(t^7)
f = a.solve_linear_de(prec=5, b=b, f0=3/5)
f
f.derivative() - a*f - b

Fourier series of periodic functions

Let \(f\) be a real-valued periodic function of period \(2L\). The Fourier series of \(f\) is

\[S(x) = \frac{a_0}{2} + \sum_{n=1}^\infty \left[a_n\cos\left(\frac{n\pi x}{L}\right) + b_n\sin\left(\frac{n\pi x}{L}\right)\right]\]

where

\[a_n = \frac{1}{L}\int_{-L}^L f(x)\cos\left(\frac{n\pi x}{L}\right) dx,\]

and

\[b_n = \frac{1}{L}\int_{-L}^L f(x)\sin\left(\frac{n\pi x}{L}\right) dx,\]

The Fourier coefficients \(a_n\) and \(b_n\) are computed by declaring \(f\) as a piecewise-defined function over one period and invoking the methods fourier_series_cosine_coefficient and fourier_series_sine_coefficient, while the partial sums are obtained via fourier_series_partial_sum:

sage: f = piecewise([((0,pi/2), -1), ((pi/2,pi), 2)])
sage: f.fourier_series_cosine_coefficient(0)
1
sage: f.fourier_series_sine_coefficient(5)
-6/5/pi
sage: s5 = f.fourier_series_partial_sum(5); s5
-6/5*sin(10*x)/pi - 2*sin(6*x)/pi - 6*sin(2*x)/pi + 1/2
sage: plot(f, (0,pi)) + plot(s5, (x,0,pi), color='red')
Graphics object consisting of 2 graphics primitives

>>> from sage.all import *
>>> f = piecewise([((Integer(0),pi/Integer(2)), -Integer(1)), ((pi/Integer(2),pi), Integer(2))])
>>> f.fourier_series_cosine_coefficient(Integer(0))
1
>>> f.fourier_series_sine_coefficient(Integer(5))
-6/5/pi
>>> s5 = f.fourier_series_partial_sum(Integer(5)); s5
-6/5*sin(10*x)/pi - 2*sin(6*x)/pi - 6*sin(2*x)/pi + 1/2
>>> plot(f, (Integer(0),pi)) + plot(s5, (x,Integer(0),pi), color='red')
Graphics object consisting of 2 graphics primitives

f = piecewise([((0,pi/2), -1), ((pi/2,pi), 2)])
f.fourier_series_cosine_coefficient(0)
f.fourier_series_sine_coefficient(5)
s5 = f.fourier_series_partial_sum(5); s5
plot(f, (0,pi)) + plot(s5, (x,0,pi), color='red')
_images/calculus-1.svg