Fourier Series

A Fourier series is a decomposition of any periodic function or periodic signal into the sum of a (possibly infinite) set of simple oscillating functions, namely, sines and cosines (or, equivalently, complex exponentials). Fourier series are also central to the original proof of the Nyquist-Shannon sampling theorem. The study of Fourier series is a branch of Fourier and harmonic analysis.
The Fourier series is named in honour of a French mathematician and polician Jean-Baptise Joseph Fourier (1768 - 1830). He was twice a governor (prefect) of Genoble, and imprisoned briefly during the Terror but in 1795 was appointed to the École Normale Supérieure, and subsequently succeeded Joseph-Louis Lagrange at the École Polytechnique. Fourier accompanied Napoleon Bonaparte on his Egyptian expedition in 1798, as scientific adviser, and was appointed secretary of the Institut d'Égypte. Cut off from France by the English fleet, he organized the workshops on which the French army had to rely for their munitions of war. He also contributed several mathematical papers to the Egyptian Institute (also called the Cairo Institute) which Napoleon founded.

 

 

Sturm-Liouville Problems

Computing the solution of some given equation is one of the fundamental problems of numerical analysis. If a real-valued function is differentiable and the derivative is known, then Newton's method may be used to find an approximation to its roots.

 

 

Singular Sturm-Liouville Problems

 

 

 

Green Functions

Computing the solution of some given equation is one of the fundamental problems of numerical analysis. If a real-valued function is differentiable and the derivative is known, then Newton's method may be used to find an approximation to its roots.

 

Orthogonal Expansions

Numerical integration methods can prove useful if the integrand is known only at certain points or the antiderivate is very difficult or even mpossible to find. In education, calculating numerical methods by hand may be useful in some cases, but computer programs are usually better suited in finding patterns and comparing different methods. In the next example, three numerical integration methods are implemented in Sage: the midpoint rule, the trapezoidal rule and Simpson's rule. The differences between the exact value of integration and the approximation are tabulated by the number of subintervals n

 

 

Bessel Expansion

The Bessel functions \( J_{\nu} (x) \) are used to express eigenfunctions \( \phi_n (x) = J_{\nu} (x\,\alpha_{\nu , n}/\ell ) \) of the singular Sturm-Liouville problem:

\[ \frac{{\text d}}{{\text d}x} \left( x\,\frac{{\text d}y}{{\text d}x} \right) + \left( \lambda \,x - \frac{\nu^2}{x} \right) y= 0 , \qquad y(0) <\infty , \quad y(\ell )=0. \]

Let \( \alpha_{\nu ,x} \) be the n-th positive root of the equation \( J_{\nu} (x) =0, \) then \( \lambda_{\nu , n} = \alpha^2_{\nu ,n} /\ell^2 \) is the corresponding eigenvalue.

 

 

Chebyshev Expansion

Computing the solution of some given equation is one of the fundamental problems of numerical analysis. If a real-valued function is differentiable and the derivative is known, then Newton's method may be used to find an approximation to its roots.

 

 

Legendre Expansion

Computing the solution of some given equation is one of the fundamental problems of numerical analysis. If a real-valued function is differentiable and the derivative is known, then Newton's method may be used to find an approximation to its roots.

 

 

Hermite Expansion

Computing the solution of some given equation is one of the fundamental problems of numerical analysis. If a real-valued function is differentiable and the derivative is known, then Newton's method may be used to find an approximation to its roots.

 

 

Laguerre Expansion

Computing the solution of some given equation is one of the fundamental problems of numerical analysis. If a real-valued function is differentiable and the derivative is known, then Newton's method may be used to find an approximation to its roots.

sage: f(x) = x^2-3
sage: a = 0.0
sage: b = 2.0
sage: table = []
sage: exact = integrate(f(x), x, a, b)
sage: for n in [4, 10, 20, 50, 100]:
....:     h = (b-a)/n
....:     midpoint = sum([f(a+(i+1/2)*h)*h for i in range(n)])
....:     trapezoid = h/2*(f(a) + 2*sum([f(a+i*h) for i in range(1,n)])
....:     + f(b))
....:     simpson = h/3*(f(a) + sum([4*f(a+i*h) for i in range(1,n,2)])
....:     + sum([2*f(a+i*h) for i in range (2,n,2)]) + f(b))
....:     table.append([n, h.n(digits=2), (midpoint-exact).n(digits=6),
....:     (trapezoid-exact).n(digits=6), (simpson-exact).n(digits=6)])
....:     html.table(table, header=["n", "h", "Midpoint rule",
....:     "Trapezoidal rule", "Simpson's rule"])

There are also built-in methods for numerical integration in Sage. For instance, it is possible to automatically produce piecewise-defined
line functions defined by the trapezoidal rule or the midpoint rule. These functions can be used to visualize different geometric interpretations of the numerical integration methods. In the next example, midpoint rule is used to calculate an approximation for the definite integral of the function
f(x) =x^2 -5x + 10
over the interval [0;10] using six subintervals

sage: f(x) = x^2-5*x+10
sage: f = Piecewise([[(0,10), f]])
sage: g = f.riemann_sum(6, mode="midpoint")
sage: F = f.plot(color="blue")
sage: R = add([line([[a,0],[a,f(x=a)],[b,f(x=b)],[b,0]], color="red") for (a,b), f in g.list()]
sage: show(F+R)

The Trapezoid Rule approximates the area under a given curve by finding the area under a linear approximation to the curve. The linear segments are given by the secant lines through the endpoints of the subintervals: for f(x)=x^2+1 on [0,2] with 4 subintervals, this looks like:

sage: x = var('x') 
sage: f1(x) = x^2 + 1 
sage: f = Piecewise([[(0,2),f1]]) 
sage: trapezoid_sum = f.trapezoid(4) 
sage: P = f.plot(rgbcolor=(0,0,1), plot_points=40)
sage: Q = trapezoid_sum.plot(rgbcolor=(1,0,0), plot_points=40) 
sage: L = add([line([[a,0],[a,f(x=a)]], rgbcolor=(1,0,0)) for (a,b), f in trapezoid_sum.list()])
sage: M = line([[2,0],[2,f1(2)]], rgbcolor=(1,0,0)) 
sage: show(P + Q + L + M)

To integrate the function cos(2*x)^4 +sin(x) we do

sage: numerical_integral(lambda x: cos(2*x)^5 + sin(x),  0, 2*pi)
(2.3561944901923457, 5.322071926740605e-14)

The input can be any callable:

sage: numerical_integral(lambda x: cos(2*x)^5 + sin(x),  0, 2*pi)
(2.3561944901923457, 5.322071926740605e-14)

We check this with a symbolic integration:

sage: (cos(2*x)^4+sin(x)).integral(x,0,2*pi)
3/4*pi

It is possible to integrate on infinite intervals as well by using +Infinity or -Infinity in the interval argument. For example:

sage: f = exp(-2*x)
sage: numerical_integral(f, 0, +Infinity)  
(0.4999999999993456, 2.5048509541974486e-07)

We can also numerically integrate symbolic expressions using either this function (which uses GSL) or the native integration (which uses Maxima):

sage: (x^3*sin(1/x)).nintegral(x,0,pi)
(9.874626194990983, 6.585622218041451e-08, 315, 0)

Fourier series of periodic functions

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

\[f(x) \sim \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]\]

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)
-3*cos(x)/pi - 3*sin(2*x)/pi + sin(x)/pi + 1/4

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.

Fourier Series

There are two equivalent forms of Fourier series---one is complex form

\[ f(x) \,\sim\, \sum_{k=-\infty}^{\infty} \alpha_k e^{k{\bf j} \pi x/\ell} \]

where

\[ \alpha_k = \frac{1}{2\ell} \int_{-\ell}^{\ell} f(x)\, e^{-k{\bf j} \pi x/\ell} \,{\text d} x , \qquad k=0, \pm 1, \pm 2, \ldots ; \]

and another is trigonometric form

\[ f(x) \,\sim\, \frac{a_0}{2} + \sum_{k=1}^{\infty} \left[ a_k \cos \frac{k \pi x}{\ell} + b_k \sin \frac{k \pi x}{\ell} \right] , \]

where

\begin{align*} a_0 &= \frac{1}{\ell} \int_{-\ell}^{\ell} f(x)\,{\text d} x , \\ a_k &= \frac{1}{\ell} \int_{-\ell}^{\ell} f(x)\, \cos \frac{k \pi x}{\ell} \,{\text d} x , \qquad k= 1, 2, 3, \ldots ; \\ b_k &= \frac{1}{\ell} \int_{-\ell}^{\ell} f(x)\, \sin \frac{k \pi x}{\ell} \,{\text d} x , \qquad k= 1, 2, 3, \ldots . \end{align*}

which are based on Euler formulas:

\[ \sin \theta = \frac{1}{2{\bf j}} \,e^{{\bf j}\theta} - \frac{1}{2{\bf j}} \,e^{-{\bf j}\theta} = \Im \,e^{{\bf j}\theta} , \qquad \cos \theta = \frac{1}{2} \,e^{{\bf j}\theta} - \frac{1}{2} \,e^{-{\bf j}\theta} = \Re \,e^{{\bf j}\theta} . \]

If \(f(x)\) is a piecewise-defined polynomial function on \(-L<x<L,\) then the Fourier series converges pointwise except points of discontionuities. More general, if \(f(x)\) is a piecewise-defined polynomial function on \(-L<x<L,\) then the Fourier series

\[f(x) \sim \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]\]

converges.

Example 2.5.1. Consider a function \( f(x) = x^2 \) on the interval \( [0,2]. \) Expanding this function periodically with the period \( T= 2 = 2{\ell} \) , we get

sage: var('x');
            sage: f = x^2
            sage: fp=Piecewise([[(-1,0),(x+2)^2], [(0,1), f]])
sage: L=1
            sage: fp.plot()

Now we calculate Fourier coefficients of the given function \( f(x) = x^2 .\) To achive this, we first introduce the parameter \( k, \) and then define it as an integer.

sage: var('k') 
sage: from sage.symbolic.assumptions import GenericDeclaration
sage: decl = GenericDeclaration(k, 'integer')
sage: decl.assume()
sage: simplify(sin(k*pi))
0

So Sage knows that k is an integer. Now we calculate Euler--Fourier coefficients:

sage: ak=integrate((x+2)^2*cos(k*pi*x),x,-1,0)+integrate((x)^2*cos(k*pi*x),x,0,1); ak
4/(pi^2*k^2)
sage: ak.simplify_full()
4/(pi^2*k^2)
sage: bk=integrate((x+2)^2*sin(k*pi*x),x,-1,0)+integrate((x)^2*sin(k*pi*x),x,0,1); 
sage: bk.simplify_full()
-2*(2*pi^2*k^2 - 1)/(pi^3*k^3) - 2/(pi^3*k^3)
sage: a0=integrate((x+2)^2,x,-1,0)+integrate((x)^2,x,0,1);a0
                  8/3  
sage: alphak=1/(2)*integrate(f*e^(-i*k*pi*x/L),x,-1,1); 
sage: alphak.simplify_full()
-1/2*(-I*pi^2*k^2 - 2*pi*k + (I*pi^2*k^2 - 2*pi*k - 2*I)*e^(2*I*pi*k) + 2*I)*e^(-I*pi*k)/(pi^3*k^3)

Now we build and plot the partial Fourier series with finite number of terms (10 and 50):

sage: plot(a0/2+sum(ak*cos(k*pi*x/L)+bk*sin(k*pi*x/L),k,1,10), (x, 0, 2))
                    sage: plot(a0/2+sum(ak*cos(k*pi*x/L)+bk*sin(k*pi*x/L),k,1,50), (x, 0, 2))
                  

Also, we can plot partial sums as

sage: s10=a0/2+sum(ak*cos(k*pi*x/L)+bk*sin(k*pi*x/L),k,1,10)
sage: plot(s10,(x,0,2))
sage: s10.plot()
                  

Finally, we plot partial Cesaro sums to surpress Gibbs phenomena:

sage: show(fp.plot_fourier_series_partial_sum_cesaro(10,1,0,2))
                    sage: show(fp.plot_fourier_series_partial_sum_cesaro(50,1,0,2))

Example 2.5.2. Consider a piecewise continuous function

sage: f1 = lambda t: -1
                     sage: f2 = lambda t: 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)
                     -3*cos(t)/pi - 3*sin(2*t)/pi + sin(t)/pi + 1/4
                   

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.

To plot Fourier partial sums, type:

sage: show(f.plot_fourier_series_partial_sum(15,pi,-5,5))
                    

 

 

Complex Fourier Series

The alternative form to the classical Fourier series is its complex counterpart: .

\[ f(x) \,\sim\, \sum_{k=-\infty}^{\infty} \alpha_k e^{k{\bf j} \pi x/\ell} \]

where

\[ \alpha_k = \frac{1}{2\ell} \int_{-\ell}^{\ell} f(x)\, e^{-k{\bf j} \pi x/\ell} \,{\text d} x , \qquad k=0, \pm 1, \pm 2, \ldots ; \]

We repeat computations for our previous example. Since the case k=0 requires a special care, we calculate the corresponding Fourier coefficient for k=0, as well as for positive and negative values of k:

sage: var('x')
sage: L=1;
sage: fp=Piecewise([[(-L,0),(x+2)^2], [(0,L), x^2]])
sage: alphak=1/2*integrate((x+2)^2*e^(-i*k*pi*x/L),x,-1,0)+1/2*integrate(x^2*e^(-i*k*pi*x/L),x,0,1)
sage: alphak.simplify_full()
-1/2*(-I*pi^2*k^2 - 2*pi*k + (I*pi^2*k^2 + 2*pi*k - 2*I)*e^(2*I*pi*k) + (-4*I*pi^2*k^2 - 4*pi*k)*e^(I*pi*k) + 2*I)*e^(-I*pi*k)/(pi^3*k^3)
sage: alpha0=1/2*integrate((x+2)^2,x,-1,0)+1/2*integrate(x^2,x,0,1)
sage: alphakm=1/2*integrate((x+2)^2*e^(i*k*pi*x/L),x,-1,0)+1/2*integrate(x^2*e^(i*k*pi*x/L),x,0,1)
sage: alphakm.simplify_full()
-1/2*(-I*pi^2*k^2 + 2*pi*k + (I*pi^2*k^2 - 2*pi*k - 2*I)*e^(2*I*pi*k) + (4*I*pi^2*k^2 - 4*pi*k)*e^(I*pi*k) + 2*I)*e^(-I*pi*k)/(pi^3*k^3)

Then we build and plot partial sums in the same way as we did for regular Fourier series:

sage: complex10=alpha0 +sum(alphak*e^(i*k*pi*x/L),k,1,10) + sum(alphakm*e^(-i*k*pi*x/L),k,1,10)
sage: complex10.plot()

 

Manual Evaluation

We can build a Fourier series based on our knowledge of Euler-Fourier coefficients.

sage: ak=integrate((x+2)^2*cos(k*pi*x),x,-1,0)+integrate((x)^2*cos(k*pi*x),x,0,1);
sage: a0=integrate((x+2)^2,x,-1,0)+integrate(x^2,x,0,1);
sage: bk=integrate((x+2)^2*sin(k*pi*x),x,-1,0)+integrate((x)^2*sin(k*pi*x),x,0,1);
sage: s10=a0/2+sum(ak*cos(k*pi*x/L)+bk*sin(k*pi*x/L),k,1,10)
     

Sine and Cosine Fourier Approximations

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]
    

Cesaro Summation

The Cesàro sum is defined as the limit of the arithmetic mean of the partial sums of the series.

Cesàro summation is named for the Italian analyst Ernesto Cesàro (1859–1906).

:

sage: show(fp.plot_fourier_series_partial_sum_cesaro(10,1,0,2))
            # or
            sage: c10=a0/2+sum(ak*cos(k*pi*x/L)*(11-k)/11+bk*sin(k*pi*x/L)*(11-k)/11,k,1,10) 
            sage: show(fp.plot_fourier_series_partial_sum_cesaro(50,1,0,2))
sage: plot(c10,(x,0,2))
          

Of course, you can view the LaTeX-ed version of this using view(g.powerseries('x',0)).

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

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: [bernoulli(2*i) for i in range(1,7)]
              [1/6, -1/30, 1/42, -1/30, 5/66, -691/2730]
              sage: maxima(f).powerseries(x,0)    # TODO: write this without maxima
              'sum((-1)^i3*2^(2*i3-1)*bern(2*i3)*_SAGE_VAR_x^(2*i3)/(i3*factorial(2*i3)),i3,1,inf)
            

Gibbs Phenomenon

Let \( F_N (x) \) be the finite Fourier sum with N+1 terms:

\[ F_N (x) = \frac{a_0}{2} + \sum_{k=1}^N \left( a_k \cos \frac{k\pi x}{\ell} + b_k \sin \frac{k\pi x}{\ell} \right) \]

If a function f(x) has a discontinuity at the point \( x_0 \) of amount \( f(x_0 +0) - f(x_0 -0) \) then finite Fourier sums experience overshoot and undershoot in a neighborhood of this point. These undershoots or overshoots cannot be eliminated by increasing the the number of terms in the finite Fourier sum and they approach with \( N\mapsto \infty \) the value

\[ 1.1789797444721675\ldots \left\vert f(x_0 +0) - f(x_0 -0) \right\vert . \]

Such behavior of Fourier series is usually referred to as the Gibbs phenomenon, which was first noticed and analyzed by the English mathematician Henry Wilbraham (1825--1883) in 1848. The term "Gibbs phenomenon" was introduced by the American mathematician Maxime Bocher in 1906. The 'magic' number is the approximation of

\[ \frac{2}{\pi} \,\mbox{Si} (\pi ) \approx 1.1789797444721675\ldots \qquad\mbox{where } \mbox{Si} (x) = \int_0^x \frac{\sin t}{t} \,{\text d}t \]

is sine integral (special function), which is build-in Sage as sin_integral(x).

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)
      

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
    

Applications

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: 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).

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
        

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)
        

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
      

is another way.

Fourier Sine and Cosine Series

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

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)
          

If you have Octave and gnuplot installed,

sage: 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 horizonal 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
            

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)]
            

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)
          

Fourier series of periodic functions

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

\[f(x) \sim \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]\]

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)
            -3*cos(x)/pi - 3*sin(2*x)/pi + sin(x)/pi + 1/4
          

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.