Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
Download

Chapter 4 - Calculus

37 views
unlisted
ubuntu2004-eol
Kernel: SageMath 9.7

Chapter 4 - Calculus (An intro to the zoo of functions via SageMath)


4.1. Introduction to Polynomials (Evaluation and Division)
4.2. Interpolation and applications
4.3. Limits of sequences and functions
4.4. Continuous functions and applications
4.5. Derivatives
4.6. Elementary theorems on derivatives
4.7. Series
4.8. Taylor series
4.9. Integration

This chapter is devoted to introductory notions from calculus and analysis. In calculus the fundamental objects that we deal with are functions.

Below we begin with material on polynomial functions. Recall that a function P(x)P(x) is called polynomial, if it has the form P(x)=anxn+an1xn1++a1x+a0 P(x)=a_{n}x^{n}+a_{n-1}x^{n-1}+\cdots+a_{1}x+a_{0} where nn is a nonnegative integer, and the numbers a0a_0, a1a_1, \ldots, ana_n are constants, called the coefficients of the polynomial. If the leading coefficient ana_{n} is non-zero, an0a_{n}\neq 0, then nn is called the degree of PP. For example, a line f(x)=ax+bf(x)=ax+b with slope a0a\neq 0 is a polynomial of degree 1, hence such a polynomial is a linear function. On the other hand, a quadratic function (parabola) g(x)=ax2+bx+cg(x)=ax^2+bx+c, with a0a\neq 0 is a polynomial of degree 2.

SageMath provides an enjoyable framework to study functions and their applications, and the same applies for polynomials. We begin with an important notion in the theory of polynomials, the so-called polynomial interpolation\mathbb{polynomial \ interpolation}, having a variety of interesting applications in numerical analysis (numerical methods, and hence also to other sciences).

4.1. Introduction to Polynomials (Evaluation and Division)

As we proncounced in this section we will focus on interpolation of polynomials. As we will verify below, for most of the cases that we are interested in, Sage succesfully applies and simplifies the computations in such procedures (as most of the mathematical packages, e.g., Mathematica, Matlab or Maple). Nowadays it is well know that such packages become essentiall for realizing the beauty of numerical methods. Therefore, we should expect that most of them contain well-fixed routines to treat interpolation methods, and our goal below is to describe the situation for SageMath (at least, and this has many similarities with the implementation in Mathematica of polynomial interpolation).

The first series of exercises is about a classical topic, namely the classical Horner scheme, which is about the division of polynomials with applications in their evaluation.


Example

Given the polynomial p(x)=i=0naixi=a0+a1x+a2x2+a3x3++anxn= p(x)=\sum _{i=0}^{n}a_{i}x^{i}=a_{0}+a_{1}x+a_{2}x^{2}+a_{3}x^{3}+\cdots +a_{n}x^{n} = =a0+x(a1+x(a2+x(a3++x(an1+xan)))) \quad = a_{0}+x (a_{1}+x ( a_{2}+x(a_{3}+\cdots +x(a_{n-1}+x a_{n})\cdots )))

define coefficients bib_i as follows:

bn:=anbn1:=an1+bnx0b1:=a1+b2x0b0:=a0+b1x0\begin{aligned}b_{n}&:=a_{n}\\b_{n-1}&:=a_{n-1}+b_{n}x_{0}\\ & \dots \quad \quad \quad & \\ b_{1} &:=a_{1}+b_{2}x_{0}\\b_{0}&:=a_{0}+b_{1}x_{0}\end{aligned}

Then p(x)=(b1+b2x+b3x2+b4x3++bn1xn2+bnxn1)(xx0)+b0, p(x)=(b_{1}+b_{2}x+b_{3}x^{2}+b_{4}x^{3}+\cdots +b_{n-1}x^{n-2}+b_{n}x^{n-1})(x-x_{0})+b_{0}, and it follows that p(x0)=b0. p(x_0) = b_0. Let us now use a bit of programming in Sage, in order to implement the situation.

def horner_division(p, x0): '''p is the list of coeficients starting from degree 0''' n = len(p) b = [0] * (n-1) a = p[n-1] for i in range(n-2, -1, -1): # i goes from n - 2 to 0 b[i] = a a = p[i] + x0*a return b, a # example # p(x) = 1 + 2*x - 3*x**2 + x**3 p = [1, 2, -3, 1] x0 = 2 q, b0 = horner_division(p, x0) print("Quotient:", q) print("Remainder:", b0)
Quotient: [0, -1, 1] Remainder: 1

In other words, x33x2+2x+1=(x2x)(x2)+1 x^3 - 3x^2 + 2x + 1 = (x^2 - x)(x - 2) + 1

Exercise

Use the routine horner_division{\tt{horner\_division}} constructed above, to determine all roots of the polynomial p(x)=x3+103x2+187x1515 p(x) = x^3 + 103x^2 + 187x - 1515.

Solution:

p = [-1515.0, 187.0, 103.0, 1.0] # start by testing small integer values of x0 x0 = -1 q, r = horner_division(p, x0) print(q, r)
[85.0000000000000, 102.000000000000, 1.00000000000000] -1600.00000000000

Remark

We can perform algebraic operations on polynomials in python using the 'numpy' library.

from numpy.polynomial.polynomial import polymul, polydiv, polyadd, polypow p = polymul(polymul([101,1],[-3,1]),[5,1]) print(p) # coeficients of (1 + x)^n are binomial coefficients q = polypow([1,1], 10) print(q)
[-1.515e+03 1.870e+02 1.030e+02 1.000e+00] [ 1. 10. 45. 120. 210. 252. 210. 120. 45. 10. 1.]

4.2. Interpolation and applications

Often, we may have some data but the function f(x)f(x) that generates the data is unknown. In such problems we try to fit a certain class of functions to the data. The process of fitting a function to some given data is called interpolation. The most usual class of functions for such a procedure are the polynomials (this is because polynomials have the nice property of approximating any continuous function -- Weierstrass Approximation Theorem). The process of fitting a polynomial through given data is called polynomial interpolation\mathit{polynomial \ interpolation}. The next series of exercises is about the so-called Lagrange interpolation polynomials\mathit{Lagrange \ interpolation \ polynomials}. The Lagrange interpolation method provides a direct approach for determining interpolated values regardless of the data points spacing, that is, it can be fitted to unequally spaced or equally spaced data.


1) Lagrange interpolation

Consider some distinct points x0,x1,,xnx_0, x_1, \ldots, x_{n} (for simplicity you may think these points as real numbers). Next we want to find a polynomial P(x)P(x) of degree not greater than nn, which takes a prescribed value yiy_i at xix_i, for all i=0,1,,ni=0, 1, \ldots, n. Such a polynomial is given by P(x)=y00(x)++ynn(x)P(x)=y_0\ell_0(x)+\cdots+y_n\ell_{n}(x), with i(x)=(xx0)(xxi1)(xxi+1)(xxn)(xix0)(xixi1)(xixi+1)(xixn) \ell_{i}(x)=\frac{(x-x_0)\cdots(x-x_{i-1})(x-x_{i+1})\cdots(x-x_{n})}{(x_{i}-x_{0})\cdots(x_{i}-x_{i-1})(x_{i}-x_{i+1})\cdots(x_{i}-x_{n})} for any i=0,1,,ni=0, 1, \ldots, n. Since we have i(xi)=1\ell_{i}(x_i)=1 and i(xj)=0\ell_{i}(x_j)=0 for any iji\neq j, it follows that P(xi)=yi i=0,1,,n. P(x_i)=y_i \ \quad \forall i=0, 1, \ldots, n.

The polynomial P(x)P(x) is known as the Lagrange interpolation polynomial, while the polynomials i\ell_{i} are called elementary Lagrange polynomials.

For instance:

\bullet If two points (x0,y0)(x_0, y_0), (x1,y1)(x_1, y_1) are given then one has n=1n=1 and P(x)=y00(x)+y11(x)P(x)=y_{0}\ell_{0}(x)+y_{1}\ell_{1}(x) with 0(x)=xx1x0x1,1(x)=xx0x1x0. \ell_{0}(x)=\frac{x-x_1}{x_0-x_1}\,,\quad \ell_1(x)=\frac{x-x_0}{x_1-x_0}\,. \bullet for n=2n=2 and some given points x0,x1,x2x_0, x_1, x_2 the polynomials 0,1,2\ell_0, \ell_1, \ell_2 are given respectively by 0(x)=(xx1)(xx2)(x0x1)(x0x2),1(x)=(xx0)(xx2)(x1x0)(x1x2),2(x)=(xx0)(xx1)(x2x0)(x2x1).\begin{darray}{rcl} \ell_{0}(x)&=&\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}\,,\\ \ell_{1}(x)&=&\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}\,,\\ \ell_{2}(x)&=&\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}\,. \end{darray}

In an analogous way are treated the cases with n=3,4n=3, 4, etc. If yi=f(xi)y_{i}=f(x_i) for some function ff, then P(x)P(x) is referred to as the Lagrange interpolation polynomial for ff.


Exercise (Lagrange interpolation polynomial)

Using polynomial interpolation, present an approximate formula of the sine function, using the known values of sin\sin at the points 00, π6\frac{\pi}{6}, π4\frac{\pi}{4}, π3\frac{\pi}{3}, π2\frac{\pi}{2}. Make a plot comparing the interpolation polynomial with sin(x)\sin(x).

Solution

We have the table x:0π/6π/4π/3π/2y:01/22/23/21. x: \quad 0 \quad \pi/6 \quad \pi/4 \quad \pi/3 \quad \pi/2 \\ y: \quad 0 \quad 1/2 \quad \sqrt2/2 \quad \sqrt3/2 \quad 1\,. Thus by applying the method based on elementary Lagrange polynomials we can compute the interpolation polunomial P(x)P(x) for the sin function. In SageMath this goes as follows:

nodes=[(0, 0), (pi/6, 1/2), (pi/4, sqrt(2)/2), (pi/3, sqrt(3)/2), (pi/2, 1)] R=PolynomialRing(RR, "x") P=R.lagrange_polynomial(nodes); show(P) A=plot(P, 0, pi, color="blue", legend_label="$P(x)$") B=list_plot({0: 0, pi/6: 1/2, pi/4: sqrt(2)/2, pi/3: sqrt(3)/2, pi/2 : 1}, size=30, figsize=4, color="black") C=plot(sin(x), 0, pi, color="red", legend_label="$sin(x)$") show(A+B+C, figsize=4)

0.0287971124604161x40.204340696021654x3+0.0213730075288600x2+0.995626184275262x+0.000000000000000\displaystyle 0.0287971124604161 x^{4} - 0.204340696021654 x^{3} + 0.0213730075288600 x^{2} + 0.995626184275262 x + 0.000000000000000

Image in a Jupyter notebook

We see from the given plots that for the given values and from 00 to π/2\pi/2, the interpolation polynomial P(x)0.02879x40.20434x3+0.02137x2+0.99562xP(x)\approx 0.02879x^{4}−0.20434x^{3}+0.02137x^{2}+0.99562x is accurate enough.

Exercise (Elementary Lagrange polynomials)

Consider the nodes x0=0x_0=0, x1=1x_1=1, x2=4x_2=4 and the values y0=1y_0=1, y1=2y_1=2, y2=4y_2=4. Write down the corresponding Lagrange interpolation polynomial P(x)P(x), and provide its plot, together with the plots of the elementary Lagrange polynomials 0,1,2\ell_0, \ell_1, \ell_2, along with the given points.

Solution:

Three points x0,x1,x2x_0, x_1, x_2 are given, hence the Lagrange interpolation polynomial P(x)P(x) is at most of degree two. The elementary Lagrange polynomials are given by a(x):=0(x)=(xx1)(xx2)/(x0x1)(x0x2), a(x):=\ell_{0}(x)=(x-x_{1})(x-x_{2})/(x_{0}-x_{1})(x_{0}-x_{2}),

b(x):=1(x)=(xx0)(xx2)/(x1x0)(x1x2),b(x):=\ell_{1}(x)=(x-x_{0})(x-x_{2})/(x_{1}-x_{0})(x_{1}-x_{2}),c(x):=2(x)=(xx0)(xx1)/(x2x0)(x2x1).c(x):=\ell_{2}(x)=(x-x_{0})(x-x_{1})/(x_{2}-x_{0})(x_{2}-x_{1}).

Thus we compute a(x)=(x1)(x4)4a(x)=\frac{(x-1)(x-4)}{4}, b(x)=x(x4)3b(x)=-\frac{x(x-4)}{3} and c(x)=x(x1)12c(x)=\frac{x(x-1)}{12} and now we can write down the interpolation polynomial P(x)P(x), as follows: P(x)=y00(x)+y11(x)+y22(x)=112x2+1312x+1. P(x)=y_0\ell_{0}(x)+y_1\ell_{1}(x)+y_2\ell_{2}(x)=-\frac{1}{12}x^2+\frac{13}{12}x+1\,. To verify this in Sage and plot the same time the given data, we have used the cell:

nodes=[(0, 1), (1, 2), (4, 4)] R=PolynomialRing(QQ, "x") P=R.lagrange_polynomial(nodes); show(P) aa(x)=(x-1)*(x-4)/4 bb(x)=-x*(x-4)/3 cc(x)=x*(x-1)/12 a=plot(P(x), x, 0, 5, color="grey", legend_label="P") b=plot(aa, x, 0, 5, color="black", legend_label="a") c=plot(bb, x, 0, 5, color="green", legend_label="b") d=plot(cc, x, 0, 5, color="blue", legend_label="c") p=point((0, 0), size=50, color="black") pt=text("$x_0$",(0.1, -0.2),color="black", fontsize="12") q=point((1, 0), size=50, color="black") qt=text("$x_1$",(1.01, -0.4),color="black", fontsize="12") r=point((4, 0), size=50, color="black") rt=text("$x_2$",(4.01, -0.4),color="black", fontsize="12") u=point((0, 1), size=50, color="black") ut=text("$a(0)=$",(-0.3, 1.0),color="black", fontsize="12") v=point((1, 1), size=50, color="black") vt=text("$b(1)=1$",(1.3, 1.0),color="black", fontsize="12") w=point((4, 1), size=50, color="black") wt=text("$c(4)=1$",(4.3, 1.0),color="black", fontsize="12") show(a+b+c+d+p+pt+q+qt+r+rt+u+ut+v+vt+w+wt, figsize=8)

112x2+1312x+1\displaystyle -\frac{1}{12} x^{2} + \frac{13}{12} x + 1

Image in a Jupyter notebook

Exercise

Prove that the Lagrange interpolation polynomial for the data given below, is of degree two:

xi2112yi1111\begin{matrix} x_i & -2 & -1 & 1 & 2 \\ \hline y_i & 1 & -1 & -1 & 1 \end{matrix}

Solution:

There are given 4 points x0,x1,x2,x3x_0, x_1, x_2, x_3, so the Lagrange interpolation polynomial should be of degree 3\leq 3. We will use Sage to show that in fact it has degree 22:

R = PolynomialRing(QQ, "x") R.lagrange_polynomial([(-2 ,1) ,(-1 ,-1) ,(1 , -1), (2, 1)])
2/3*x^2 - 5/3

Thus the interpolation polynomial under question is given by P(x)=23x253P(x)=\frac{2}{3}x^2-\frac{5}{3}.

Exercise

Consider the function f(x)=e2x2xf(x)=e^{2x}-2x and the points x0=1x_0=1, x1=1.5x_1=1.5, x2=1.6x_2=1.6. Construct the Lagrange interpolation polynomial P(x)P(x), to approximate f(1.55)f(1.55). Next find the difference between the approximation and the real value.

Solution:

We use SageMath and the command .lagrange_polynomial{\tt{.lagrange\_polynomial}} to obtain the Lagrange interpolation polynomial. We have

nodes=[(1, e^2-2), (1.5, e^(3)-3), (1.6, e^(3.2)-3.2)] R=PolynomialRing(RR, "x") P=R.lagrange_polynomial(nodes); show(P)

31.7949518178380x256.0944178960809x+29.6885221771736\displaystyle 31.7949518178380 x^{2} - 56.0944178960809 x + 29.6885221771736

Thus P(x)31.79495x256.09441x+29.68852P(x)\approx 31.79495x^2-56.09441x+29.68852 and P(x)P(x) is of degree two. Let us now compare the approximation with the real value:

P(1.55)
19.1295461806039
f(x)=e^(2*x)-2*x f(1.55)
19.0979512814416

The difference between the approximation and the real value is given by the absolute value f(1.55)P(1.55)\left|f(1.55)-P(1.55)\right|. We compute

abs(f(1.55)-P(1.55))
0.0315948991622825

** Remark (on the use of scipy)**

We can apply the Lagrange intepolation method in SageMath also by using scipy{\tt{scipy}} and importing the scipy.intrpolate{\tt{scipy.intrpolate}} as follows:

#interpolation using scipy import scipy import scipy.interpolate from scipy.interpolate import lagrange import numpy as np xx=np.array([0, 2, 3, 4]) yy=np.array([1, 3, 2, 1/2]) polL=lagrange(xx, yy) polL
poly1d([ 0.10416667, -1.1875 , 2.95833333, 1. ])

Above we are looking for the Lagrange interpolation polynomials for the data (0,1)(0, 1), (2,3)(2, 3), (3,2)(3, 2), and (4,1/2)(4, 1/2), which we introduced in Sage in terms of xx- and yy-coordinates. The answer returned by Sage should be read as P(x)0.1041x31.1875x2+2.9583x+1P(x)\approx -0.1041x^3-1.1875x^2+2.9583x+1. Finally, note that we may also type

show(polL)

        3         20.1042 x - 1.188 x + 2.958 x + 1\displaystyle \begin{array}{l} \verb| |\verb|3|\verb| |\verb|2|\\ \verb|0.1042|\verb| |\verb|x|\verb| |\verb|-|\verb| |\verb|1.188|\verb| |\verb|x|\verb| |\verb|+|\verb| |\verb|2.958|\verb| |\verb|x|\verb| |\verb|+|\verb| |\verb|1| \end{array}

Exercise for practice

Find via the method of scipy{\tt{scipy}} the Lagrange interpolation polynomial for the data (0,e)(0, e), (0.25,ln(e))(0.25, \sqrt{\ln(e)}), (0.5,e2)(0.5, e^2), (1,ln(e2))(1, \sqrt{\ln(e^2)}).

2) Hermite interpolation

Rougly speaking Hermite interpolation polynomials are polynomials P(x)P(x) that interpolate data that may includes their first derivatives. More particular assume that are given:

  1. n+1n+1 pairwise distinctint nodes, say x0,,xnx_0, \ldots, x_n,

  2. The values y0=P(x0),,yn=P(xn)y_0=P(x_0), \ldots, y_n=P(x_n),

  3. The first derivatives y0,,yny' _ 0, \ldots, y'_ n, where yi:=P(xi)y'_ i:=P'(x_i) and in general by P(x)P'(x) we denote the (first) derivative of a function P(x)P(x).

Then we may express PP by the sum

P(x)=i=0ny0hi(1)(x)+i=0nyihi(2)(x)P(x)=\sum_{i=0}^ny_0\cdot h^{(1)}_ {i}(x)+\sum_{i=0}^ny'_ {i}\cdot h^{(2)}_ {i}(x)

where hi(1)(x)h^{(1)}_ i(x) and hi(2)(x)h^{(2)}_ i(x) (i=0,1,,n)(i=0, 1, \ldots, n) are the fundamental Hermite interpolation polynomials of first and second type, given by

hi(1)(x):=(1(xi)(xi)(xxi))(i(x))2,hi(2)(x):=(xxi)(i(x))2,h^{(1)}_ {i}(x):=\Big(1-\frac{\ell''(x_i)}{\ell'(x_i)}(x-x_i)\Big)(\ell_{i}(x))^2, \quad\quad h^{(2)}_ {i}(x):=(x-x_i)(\ell_{i}(x))^2,

respectively. Here we have (x):=(xx0)(xx1)(xxn)\ell(x):=(x-x_0)(x-x_1)\cdots(x-x_n) and i(x)\ell_i(x) are the fundamental Lagrange polynomials introduced above (for i=0,1,,ni=0, 1, \ldots, n).

Again, as in the case of Lagrange interpolation we can use the above method to interpolate some given function f(x)f(x) at some given points.


Exercise

Consider the function f(x)=sin(x)f(x)=\sin(x). Given 5 distinct points x0,x1,x2,x3,x4x_0, x_1, x_2, x_3, x_4 in the real line, write a code in Sage which will return the Hermite interpolation polynomial P(x)P(x) corresponding to these nodes and the values yi=f(xi)y_i=f(x_i) and yi=f(xi)y'_i=f'(x_i), for i=0,1,2,3,4i=0, 1, 2, 3, 4 (for simplicity you may fix some 5-tuple (x0,x1,x2,x3,x4)(x_0, x_1, x_2, x_3, x_4)).

Solution:

Since the Hermite interpolation method includes derivatives, next we are goint to use commands as derivative(f,x){\tt{derivative(f, x)}}, which gives the derivative of a function ff and derivative(f,x,n){\tt{derivative(f, x, n) }}, returning the nnth derivative of ff (see also below for more details on how one can treat derivatives in SageMath). So, let us fix the nodes

x0=2π,x1=π/2,x2=0,x3=π/2,x4=2π.x_0=-2\pi, \quad x_1=-\pi/2, \quad x_2=0, \quad x_3=\pi/2, \quad x_4=2\pi.

Below we write the code in such a way that changing only the very first line, i.e., giving different values to the 5-tuple (x0,x1,x2,x3,x4)(x_0, x_1, x_2, x_3, x_4), will immediately produce the corresponding Hermite interpolation polynomial PP, together with the plots of the function f(x)f(x) and that of P(x)P(x).

x0=-2*pi; x1=-pi/2; x2=0; x3=pi/2; x4=2*pi l(x)=(x-x0)*(x-x1)*(x-x2)*(x-x3)*(x-x4) d1l(x)=derivative(l, x) d2l(x)=derivative(l(x),x,2) l0(x)=((x-x1)*(x-x2)*(x-x3)*(x-x4))/((x0-x1)*(x0-x2)*(x0-x3)*(x0-x4)) l1(x)=((x-x0)*(x-x2)*(x-x3)*(x-x4))/((x1-x0)*(x1-x2)*(x1-x3)*(x1-x4)) l2(x)=((x-x0)*(x-x1)*(x-x3)*(x-x4))/((x2-x0)*(x2-x1)*(x2-x3)*(x2-x4)) l3(x)=((x-x0)*(x-x1)*(x-x2)*(x-x4))/((x3-x0)*(x3-x1)*(x3-x2)*(x3-x4)) l4(x)=((x-x0)*(x-x1)*(x-x2)*(x-x3))/((x4-x0)*(x4-x1)*(x4-x2)*(x4-x3)) h01(x)=(1-(d2l(x0)/d1l(x0))*(x-x0))*(l0(x))^2 h11(x)=(1-(d2l(x1)/d1l(x1))*(x-x1))*(l1(x))^2 h21(x)=(1-(d2l(x2)/d1l(x2))*(x-x2))*(l2(x))^2 h31(x)=(1-(d2l(x3)/d1l(x3))*(x-x3))*(l3(x))^2 h41(x)=(1-(d2l(x4)/d1l(x4))*(x-x4))*(l4(x))^2 h02(x)=(x-x0)*(l0(x))^2 h12(x)=(x-x1)*(l1(x))^2 h22(x)=(x-x2)*(l2(x))^2 h32(x)=(x-x3)*(l3(x))^2 h42(x)=(x-x4)*(l4(x))^2 f(x)=sin(x) y0=f(x0); y1=f(x1); y2=f(x2); y3=f(x3); y4=f(x4) dy0=derivative(f(x), x)(x=x0) dy1=derivative(f(x), x)(x=x1) dy2=derivative(f(x), x)(x=x2) dy3=derivative(f(x), x)(x=x3) dy4=derivative(f(x), x)(x=x4) P(x)=y0*h01(x)+y1*h11(x)+y2*h21(x)+y3*h31(x)+y4*h41(x)+dy0*h02(x)+dy1*h12(x)+dy2*h22(x)+dy3*h32(x)+dy4*h42(x) show(P(x)) a=plot(f(x), x, x0, x4, color="darkgray", thickness=3, legend_label=" $f(x)$") b=plot(P(x), x, x0, x4, color="darkblue", thickness=1, linestyle="-.", legend_label="$P(x)$") show(a+b)

(2π+x)2(2πx)2(π+2x)2(π2x)2x16π8(2π+x)2(2πx)(π+2x)2(π2x)2x214400π8+(2π+x)(2πx)2(π+2x)2(π2x)2x214400π816(2π+x)2(2πx)2(π2x)2x2(41(π+2x)π+15)3375π8+16(2π+x)2(2πx)2(π+2x)2x2(41(π2x)π+15)3375π8\displaystyle \frac{{\left(2 \, \pi + x\right)}^{2} {\left(2 \, \pi - x\right)}^{2} {\left(\pi + 2 \, x\right)}^{2} {\left(\pi - 2 \, x\right)}^{2} x}{16 \, \pi^{8}} - \frac{{\left(2 \, \pi + x\right)}^{2} {\left(2 \, \pi - x\right)} {\left(\pi + 2 \, x\right)}^{2} {\left(\pi - 2 \, x\right)}^{2} x^{2}}{14400 \, \pi^{8}} + \frac{{\left(2 \, \pi + x\right)} {\left(2 \, \pi - x\right)}^{2} {\left(\pi + 2 \, x\right)}^{2} {\left(\pi - 2 \, x\right)}^{2} x^{2}}{14400 \, \pi^{8}} - \frac{16 \, {\left(2 \, \pi + x\right)}^{2} {\left(2 \, \pi - x\right)}^{2} {\left(\pi - 2 \, x\right)}^{2} x^{2} {\left(\frac{41 \, {\left(\pi + 2 \, x\right)}}{\pi} + 15\right)}}{3375 \, \pi^{8}} + \frac{16 \, {\left(2 \, \pi + x\right)}^{2} {\left(2 \, \pi - x\right)}^{2} {\left(\pi + 2 \, x\right)}^{2} x^{2} {\left(\frac{41 \, {\left(\pi - 2 \, x\right)}}{\pi} + 15\right)}}{3375 \, \pi^{8}}

Image in a Jupyter notebook

In fact one case see the explicit form of Hermite interpolation polynomials of 1st and 2nd type, after typing (for instance, for the 1st type)

show([h01(x), h11(x), h21(x), h31(x)])

[(2πx)2(π+2x)2(π2x)2x2(109(2π+x)π+30)432000π8,16(2π+x)2(2πx)2(π2x)2x2(41(π+2x)π+15)3375π8,(2π+x)2(2πx)2(π+2x)2(π2x)216π8,16(2π+x)2(2πx)2(π+2x)2x2(41(π2x)π+15)3375π8]\displaystyle \left[\frac{{\left(2 \, \pi - x\right)}^{2} {\left(\pi + 2 \, x\right)}^{2} {\left(\pi - 2 \, x\right)}^{2} x^{2} {\left(\frac{109 \, {\left(2 \, \pi + x\right)}}{\pi} + 30\right)}}{432000 \, \pi^{8}}, \frac{16 \, {\left(2 \, \pi + x\right)}^{2} {\left(2 \, \pi - x\right)}^{2} {\left(\pi - 2 \, x\right)}^{2} x^{2} {\left(\frac{41 \, {\left(\pi + 2 \, x\right)}}{\pi} + 15\right)}}{3375 \, \pi^{8}}, \frac{{\left(2 \, \pi + x\right)}^{2} {\left(2 \, \pi - x\right)}^{2} {\left(\pi + 2 \, x\right)}^{2} {\left(\pi - 2 \, x\right)}^{2}}{16 \, \pi^{8}}, \frac{16 \, {\left(2 \, \pi + x\right)}^{2} {\left(2 \, \pi - x\right)}^{2} {\left(\pi + 2 \, x\right)}^{2} x^{2} {\left(\frac{41 \, {\left(\pi - 2 \, x\right)}}{\pi} + 15\right)}}{3375 \, \pi^{8}}\right]

Let us focus on these polynomials for a few, for instance on h0(1)(x)h_{0}^{(1)}(x). Let us find its coefficients:

R.<x> = PolynomialRing(QQ) h01(x).coefficients()
[[31/13500/pi^2, 2], [-139/108000/pi^3, 3], [-677/36000/pi^4, 4], [1519/144000/pi^5, 5], [181/4500/pi^6, 6], [-407/18000/pi^7, 7], [-47/6750/pi^8, 8], [109/27000/pi^9, 9]]

This means that h0(1)(x)h_0^{(1)}(x) is of degree 99 and has no term of first order and no constant coefficient, as we can see also based on the command expand{\tt{expand}}:

expand(h01(x))
31/13500*x^2/pi^2 - 139/108000*x^3/pi^3 - 677/36000*x^4/pi^4 + 1519/144000*x^5/pi^5 + 181/4500*x^6/pi^6 - 407/18000*x^7/pi^7 - 47/6750*x^8/pi^8 + 109/27000*x^9/pi^9

Let us now check h0(2)(x)h_0^{(2)}(x)

h02(x).coefficients()
[[1/1800/pi, 2], [-1/3600/pi^2, 3], [-11/2400/pi^3, 4], [11/4800/pi^4, 5], [1/100/pi^5, 6], [-1/200/pi^6, 7], [-1/450/pi^7, 8], [1/900/pi^8, 9]]
expand(h02(x))
1/1800*x^2/pi - 1/3600*x^3/pi^2 - 11/2400*x^4/pi^3 + 11/4800*x^5/pi^4 + 1/100*x^6/pi^5 - 1/200*x^7/pi^6 - 1/450*x^8/pi^7 + 1/900*x^9/pi^8

Lets us now plot the Hermite polynomials of first and of second type seperately:

p1=plot(h01(x), x, x0, x4, color="green") p2=plot(h11(x), x, x0, x4, color="steelblue") p3=plot(h21(x), x, x0, x4, color="brown") p4=plot(h31(x), x, x0, x4, color="red") show(p1+p2+p3+p4)
Image in a Jupyter notebook
t1=plot(h02(x), x, x0, x4, color="pink", figsize=8) t2=plot(h12(x), x, x0, x4, color="darkblue") t3=plot(h22(x), x, x0, x4, color="black") t4=plot(h32(x), x, x0, x4, color="darkred") show(t1+t2+t3+t4)
Image in a Jupyter notebook

We proceed with another example:

x0=0; x1=0.2; x2=1/e; x3=e; x4=pi**2 l(x)=(x-x0)*(x-x1)*(x-x2)*(x-x3)*(x-x4) d1l(x)=derivative(l, x) d2l(x)=derivative(l(x),x,2) l0(x)=((x-x1)*(x-x2)*(x-x3)*(x-x4))/((x0-x1)*(x0-x2)*(x0-x3)*(x0-x4)) l1(x)=((x-x0)*(x-x2)*(x-x3)*(x-x4))/((x1-x0)*(x1-x2)*(x1-x3)*(x1-x4)) l2(x)=((x-x0)*(x-x1)*(x-x3)*(x-x4))/((x2-x0)*(x2-x1)*(x2-x3)*(x2-x4)) l3(x)=((x-x0)*(x-x1)*(x-x2)*(x-x4))/((x3-x0)*(x3-x1)*(x3-x2)*(x3-x4)) l4(x)=((x-x0)*(x-x1)*(x-x2)*(x-x3))/((x4-x0)*(x4-x1)*(x4-x2)*(x4-x3)) h01(x)=(1-(d2l(x0)/d1l(x0))*(x-x0))*(l0(x))^2 h11(x)=(1-(d2l(x1)/d1l(x1))*(x-x1))*(l1(x))^2 h21(x)=(1-(d2l(x2)/d1l(x2))*(x-x2))*(l2(x))^2 h31(x)=(1-(d2l(x3)/d1l(x3))*(x-x3))*(l3(x))^2 h41(x)=(1-(d2l(x4)/d1l(x4))*(x-x4))*(l4(x))^2 h02(x)=(x-x0)*(l0(x))^2 h12(x)=(x-x1)*(l1(x))^2 h22(x)=(x-x2)*(l2(x))^2 h32(x)=(x-x3)*(l3(x))^2 h42(x)=(x-x4)*(l4(x))^2 f(x)=sin(x) y0=f(x0); y1=f(x1); y2=f(x2); y3=f(x3); y4=f(x4) dy0=derivative(f(x), x)(x=x0) dy1=derivative(f(x), x)(x=x1) dy2=derivative(f(x), x)(x=x2) dy3=derivative(f(x), x)(x=x3) dy4=derivative(f(x), x)(x=x4) P(x)=y0*h01(x)+y1*h11(x)+y2*h21(x)+y3*h31(x)+y4*h41(x)+dy0*h02(x)+dy1*h12(x)+dy2*h22(x)+dy3*h32(x)+dy4*h42(x) a=plot(f(x), x, x0, x4, color="darkgray", thickness=3, legend_label=" $f(x)$") b=plot(P(x), x, x0, x4, color="darkblue", thickness=1, linestyle="-.", legend_label="$P(x)$") show(a+b)
Image in a Jupyter notebook

As we can see P(x)P(x) has much more complicated form than the previous case:

P(x).simplify()
25.0*(pi^2 - x)^2*(x - e)^2*(x - e^(-1))^2*(x - 0.2)^2*x/pi^4 + (pi^2 - x)^2*(x - e)^2*(x - e^(-1))*(x - 0.2)^2*x^2*cos(e^(-1))*e^2/((pi^2 - e^(-1))^2*(e - e^(-1))^2*(e^(-1) - 0.2)^2) + (pi^2 - x)^2*(x - e)*(x - e^(-1))^2*(x - 0.2)^2*x^2*cos(e)*e^(-2)/((pi^2 - e)^2*(e - e^(-1))^2*(e - 0.2)^2) - (pi^2 - x)^2*(x - e^(-1))^2*(x - 0.2)^2*x^2*(2*((pi^2 - e)*(e - e^(-1))*(e - 0.2) + (pi^2 - e)*(e - e^(-1))*e + (pi^2 - e)*(e - 0.2)*e - (e - e^(-1))*(e - 0.2)*e)*(x - e)*e^(-1)/((pi^2 - e)*(e - e^(-1))*(e - 0.2)) - 1)*e^(-2)*sin(e)/((pi^2 - e)^2*(e - e^(-1))^2*(e - 0.2)^2) - (pi^2 - x)^2*(x - e)^2*(x - 0.2)^2*x^2*(2*((pi^2 - e^(-1))*(e - e^(-1))*(e^(-1) - 0.2) + (pi^2 - e^(-1))*(e - e^(-1))*e^(-1) - (pi^2 - e^(-1))*(e^(-1) - 0.2)*e^(-1) - (e - e^(-1))*(e^(-1) - 0.2)*e^(-1))*(x - e^(-1))*e/((pi^2 - e^(-1))*(e - e^(-1))*(e^(-1) - 0.2)) - 1)*e^2*sin(e^(-1))/((pi^2 - e^(-1))^2*(e - e^(-1))^2*(e^(-1) - 0.2)^2) + 0.9800665778412416*(pi^2 - x)^2*(x - e)^2*(x - e^(-1))^2*(x - 0.2)*x^2/((pi^2 - 0.2)^2*(e - 0.2)^2*(-0.2*e^(-1) + 0.04000000000000001)^2) + (pi^2 - x)^2*(x - e)^2*(x - e^(-1))^2*x^2*(-0.1986693307950612*(5.0*(pi^2 - 0.2)*(e - 0.2)*(2*e^(-1) - 0.4) - 5.0*(pi^2 - 0.2)*(0.4*e - 0.08000000000000002) - 5.0*(pi^2 - 0.2)*(0.4*e^(-1) - 0.08000000000000002) - 5.0*(e - 0.2)*(0.4*e^(-1) - 0.08000000000000002))*(x - 0.2)/((pi^2 - 0.2)*(e - 0.2)*(e^(-1) - 0.2)) + 0.1986693307950612)/((pi^2 - 0.2)^2*(e - 0.2)^2*(-0.2*e^(-1) + 0.04000000000000001)^2) - (pi^2 - x)*(x - e)^2*(x - e^(-1))^2*(x - 0.2)^2*x^2*cos(pi^2)/(pi^4*(pi^2 - e)^2*(pi^2 - e^(-1))^2*(pi^2 - 0.2)^2) + (x - e)^2*(x - e^(-1))^2*(x - 0.2)^2*x^2*(2*(pi^2*(pi^2 - e)*(pi^2 - e^(-1)) + pi^2*(pi^2 - e)*(pi^2 - 0.2) + pi^2*(pi^2 - e^(-1))*(pi^2 - 0.2) + (pi^2 - e)*(pi^2 - e^(-1))*(pi^2 - 0.2))*(pi^2 - x)/(pi^2*(pi^2 - e)*(pi^2 - e^(-1))*(pi^2 - 0.2)) + 1)*sin(pi^2)/(pi^4*(pi^2 - e)^2*(pi^2 - e^(-1))^2*(pi^2 - 0.2)^2)
show(_)

25.0(π2x)2(xe)2(xe(1))2(x0.2)2xπ4+(π2x)2(xe)2(xe(1))(x0.2)2x2cos(e(1))e2(π2e(1))2(ee(1))2(e(1)0.2)2+(π2x)2(xe)(xe(1))2(x0.2)2x2cos(e)e(2)(π2e)2(ee(1))2(e0.2)2(π2x)2(xe(1))2(x0.2)2x2(2((π2e)(ee(1))(e0.2)+(π2e)(ee(1))e+(π2e)(e0.2)e(ee(1))(e0.2)e)(xe)e(1)(π2e)(ee(1))(e0.2)1)e(2)sin(e)(π2e)2(ee(1))2(e0.2)2(π2x)2(xe)2(x0.2)2x2(2((π2e(1))(ee(1))(e(1)0.2)+(π2e(1))(ee(1))e(1)(π2e(1))(e(1)0.2)e(1)(ee(1))(e(1)0.2)e(1))(xe(1))e(π2e(1))(ee(1))(e(1)0.2)1)e2sin(e(1))(π2e(1))2(ee(1))2(e(1)0.2)2+0.9800665778412416(π2x)2(xe)2(xe(1))2(x0.2)x2(π20.2)2(e0.2)2(0.2e(1)+0.04000000000000001)2+(π2x)2(xe)2(xe(1))2x2(0.1986693307950612(5.0(π20.2)(e0.2)(2e(1)0.4)5.0(π20.2)(0.4e0.08000000000000002)5.0(π20.2)(0.4e(1)0.08000000000000002)5.0(e0.2)(0.4e(1)0.08000000000000002))(x0.2)(π20.2)(e0.2)(e(1)0.2)+0.1986693307950612)(π20.2)2(e0.2)2(0.2e(1)+0.04000000000000001)2(π2x)(xe)2(xe(1))2(x0.2)2x2cos(π2)π4(π2e)2(π2e(1))2(π20.2)2+(xe)2(xe(1))2(x0.2)2x2(2(π2(π2e)(π2e(1))+π2(π2e)(π20.2)+π2(π2e(1))(π20.2)+(π2e)(π2e(1))(π20.2))(π2x)π2(π2e)(π2e(1))(π20.2)+1)sin(π2)π4(π2e)2(π2e(1))2(π20.2)2\displaystyle \frac{25.0 \, {\left(\pi^{2} - x\right)}^{2} {\left(x - e\right)}^{2} {\left(x - e^{\left(-1\right)}\right)}^{2} {\left(x - 0.2\right)}^{2} x}{\pi^{4}} + \frac{{\left(\pi^{2} - x\right)}^{2} {\left(x - e\right)}^{2} {\left(x - e^{\left(-1\right)}\right)} {\left(x - 0.2\right)}^{2} x^{2} \cos\left(e^{\left(-1\right)}\right) e^{2}}{{\left(\pi^{2} - e^{\left(-1\right)}\right)}^{2} {\left(e - e^{\left(-1\right)}\right)}^{2} {\left(e^{\left(-1\right)} - 0.2\right)}^{2}} + \frac{{\left(\pi^{2} - x\right)}^{2} {\left(x - e\right)} {\left(x - e^{\left(-1\right)}\right)}^{2} {\left(x - 0.2\right)}^{2} x^{2} \cos\left(e\right) e^{\left(-2\right)}}{{\left(\pi^{2} - e\right)}^{2} {\left(e - e^{\left(-1\right)}\right)}^{2} {\left(e - 0.2\right)}^{2}} - \frac{{\left(\pi^{2} - x\right)}^{2} {\left(x - e^{\left(-1\right)}\right)}^{2} {\left(x - 0.2\right)}^{2} x^{2} {\left(\frac{2 \, {\left({\left(\pi^{2} - e\right)} {\left(e - e^{\left(-1\right)}\right)} {\left(e - 0.2\right)} + {\left(\pi^{2} - e\right)} {\left(e - e^{\left(-1\right)}\right)} e + {\left(\pi^{2} - e\right)} {\left(e - 0.2\right)} e - {\left(e - e^{\left(-1\right)}\right)} {\left(e - 0.2\right)} e\right)} {\left(x - e\right)} e^{\left(-1\right)}}{{\left(\pi^{2} - e\right)} {\left(e - e^{\left(-1\right)}\right)} {\left(e - 0.2\right)}} - 1\right)} e^{\left(-2\right)} \sin\left(e\right)}{{\left(\pi^{2} - e\right)}^{2} {\left(e - e^{\left(-1\right)}\right)}^{2} {\left(e - 0.2\right)}^{2}} - \frac{{\left(\pi^{2} - x\right)}^{2} {\left(x - e\right)}^{2} {\left(x - 0.2\right)}^{2} x^{2} {\left(\frac{2 \, {\left({\left(\pi^{2} - e^{\left(-1\right)}\right)} {\left(e - e^{\left(-1\right)}\right)} {\left(e^{\left(-1\right)} - 0.2\right)} + {\left(\pi^{2} - e^{\left(-1\right)}\right)} {\left(e - e^{\left(-1\right)}\right)} e^{\left(-1\right)} - {\left(\pi^{2} - e^{\left(-1\right)}\right)} {\left(e^{\left(-1\right)} - 0.2\right)} e^{\left(-1\right)} - {\left(e - e^{\left(-1\right)}\right)} {\left(e^{\left(-1\right)} - 0.2\right)} e^{\left(-1\right)}\right)} {\left(x - e^{\left(-1\right)}\right)} e}{{\left(\pi^{2} - e^{\left(-1\right)}\right)} {\left(e - e^{\left(-1\right)}\right)} {\left(e^{\left(-1\right)} - 0.2\right)}} - 1\right)} e^{2} \sin\left(e^{\left(-1\right)}\right)}{{\left(\pi^{2} - e^{\left(-1\right)}\right)}^{2} {\left(e - e^{\left(-1\right)}\right)}^{2} {\left(e^{\left(-1\right)} - 0.2\right)}^{2}} + \frac{0.9800665778412416 \, {\left(\pi^{2} - x\right)}^{2} {\left(x - e\right)}^{2} {\left(x - e^{\left(-1\right)}\right)}^{2} {\left(x - 0.2\right)} x^{2}}{{\left(\pi^{2} - 0.2\right)}^{2} {\left(e - 0.2\right)}^{2} {\left(-0.2 \, e^{\left(-1\right)} + 0.04000000000000001\right)}^{2}} + \frac{{\left(\pi^{2} - x\right)}^{2} {\left(x - e\right)}^{2} {\left(x - e^{\left(-1\right)}\right)}^{2} x^{2} {\left(-\frac{0.1986693307950612 \, {\left(5.0 \, {\left(\pi^{2} - 0.2\right)} {\left(e - 0.2\right)} {\left(2 \, e^{\left(-1\right)} - 0.4\right)} - 5.0 \, {\left(\pi^{2} - 0.2\right)} {\left(0.4 \, e - 0.08000000000000002\right)} - 5.0 \, {\left(\pi^{2} - 0.2\right)} {\left(0.4 \, e^{\left(-1\right)} - 0.08000000000000002\right)} - 5.0 \, {\left(e - 0.2\right)} {\left(0.4 \, e^{\left(-1\right)} - 0.08000000000000002\right)}\right)} {\left(x - 0.2\right)}}{{\left(\pi^{2} - 0.2\right)} {\left(e - 0.2\right)} {\left(e^{\left(-1\right)} - 0.2\right)}} + 0.1986693307950612\right)}}{{\left(\pi^{2} - 0.2\right)}^{2} {\left(e - 0.2\right)}^{2} {\left(-0.2 \, e^{\left(-1\right)} + 0.04000000000000001\right)}^{2}} - \frac{{\left(\pi^{2} - x\right)} {\left(x - e\right)}^{2} {\left(x - e^{\left(-1\right)}\right)}^{2} {\left(x - 0.2\right)}^{2} x^{2} \cos\left(\pi^{2}\right)}{\pi^{4} {\left(\pi^{2} - e\right)}^{2} {\left(\pi^{2} - e^{\left(-1\right)}\right)}^{2} {\left(\pi^{2} - 0.2\right)}^{2}} + \frac{{\left(x - e\right)}^{2} {\left(x - e^{\left(-1\right)}\right)}^{2} {\left(x - 0.2\right)}^{2} x^{2} {\left(\frac{2 \, {\left(\pi^{2} {\left(\pi^{2} - e\right)} {\left(\pi^{2} - e^{\left(-1\right)}\right)} + \pi^{2} {\left(\pi^{2} - e\right)} {\left(\pi^{2} - 0.2\right)} + \pi^{2} {\left(\pi^{2} - e^{\left(-1\right)}\right)} {\left(\pi^{2} - 0.2\right)} + {\left(\pi^{2} - e\right)} {\left(\pi^{2} - e^{\left(-1\right)}\right)} {\left(\pi^{2} - 0.2\right)}\right)} {\left(\pi^{2} - x\right)}}{\pi^{2} {\left(\pi^{2} - e\right)} {\left(\pi^{2} - e^{\left(-1\right)}\right)} {\left(\pi^{2} - 0.2\right)}} + 1\right)} \sin\left(\pi^{2}\right)}{\pi^{4} {\left(\pi^{2} - e\right)}^{2} {\left(\pi^{2} - e^{\left(-1\right)}\right)}^{2} {\left(\pi^{2} - 0.2\right)}^{2}}

3) Splines

Let us finally present a series of problems related to splines and in particular (naturall) cubic splines.

Our guide here is Section 5.1.9 by Brisk Guide to Mathematics.


Exercise

Find and plot the spline for the points (k+sin(k2)/2,kcos(k2))(k+\sin(k^2)/2, k-\cos(k^2)) for k=0,1,,9k=0, 1,\ldots, 9.

Solution:

We present the solutin directly via Sage Math

pts =[(k+sin(k^2)/2, k-cos(k^2)) for k in range(10)] f=spline(pts) a=plot(f, 0, 9, color="green", figsize=8) show(points(pts)+a) show(f)
Image in a Jupyter notebook

[(0, -1), (1/2*sin(1) + 1, -cos(1) + 1), (1/2*sin(4) + 2, -cos(4) + 2), (1/2*sin(9) + 3, -cos(9) + 3), (1/2*sin(16) + 4, -cos(16) + 4), (1/2*sin(25) + 5, -cos(25) + 5), (1/2*sin(36) + 6, -cos(36) + 6), (1/2*sin(49) + 7, -cos(49) + 7), (1/2*sin(64) + 8, -cos(64) + 8), (1/2*sin(81) + 9, -cos(81) + 9)]\displaystyle \verb|[(0,|\verb| |\verb|-1),|\verb| |\verb|(1/2*sin(1)|\verb| |\verb|+|\verb| |\verb|1,|\verb| |\verb|-cos(1)|\verb| |\verb|+|\verb| |\verb|1),|\verb| |\verb|(1/2*sin(4)|\verb| |\verb|+|\verb| |\verb|2,|\verb| |\verb|-cos(4)|\verb| |\verb|+|\verb| |\verb|2),|\verb| |\verb|(1/2*sin(9)|\verb| |\verb|+|\verb| |\verb|3,|\verb| |\verb|-cos(9)|\verb| |\verb|+|\verb| |\verb|3),|\verb| |\verb|(1/2*sin(16)|\verb| |\verb|+|\verb| |\verb|4,|\verb| |\verb|-cos(16)|\verb| |\verb|+|\verb| |\verb|4),|\verb| |\verb|(1/2*sin(25)|\verb| |\verb|+|\verb| |\verb|5,|\verb| |\verb|-cos(25)|\verb| |\verb|+|\verb| |\verb|5),|\verb| |\verb|(1/2*sin(36)|\verb| |\verb|+|\verb| |\verb|6,|\verb| |\verb|-cos(36)|\verb| |\verb|+|\verb| |\verb|6),|\verb| |\verb|(1/2*sin(49)|\verb| |\verb|+|\verb| |\verb|7,|\verb| |\verb|-cos(49)|\verb| |\verb|+|\verb| |\verb|7),|\verb| |\verb|(1/2*sin(64)|\verb| |\verb|+|\verb| |\verb|8,|\verb| |\verb|-cos(64)|\verb| |\verb|+|\verb| |\verb|8),|\verb| |\verb|(1/2*sin(81)|\verb| |\verb|+|\verb| |\verb|9,|\verb| |\verb|-cos(81)|\verb| |\verb|+|\verb| |\verb|9)]|

Exercise

Plot the spline corresponding to the data (0,1)(0, 1), (1,2)(1, 2), (2,2.5)(2, 2.5), (3,1.5)(3, 1.5), (4,0.5)(4, 0.5 ), (5,0.9)(5, 0.9) and (6,0.5)(6, -0.5). Next compute its value at

the points 0.30.3, 2.42.4, 5.65.6. Can we compute also its value at the point 6.16.1? Then compute the area underneath the spline.

Solution:

pts = [(0,1),(1,2),(2,2.5),(3,1.5),(4,0.5),(5,0.9),(6,-0.5)] func=spline(pts) a=plot(func, 0, 6, color="green", figsize=4) show(points(pts)+a)
Image in a Jupyter notebook
func(0.3)
1.309765
func(2.4)
2.2493353846153847
func(5.6)
0.2534584615384622
func(6.1)
nan

Hence we cannot compute the value of the spline at points that are not included in between the given points.

Now, to compute the area underneath the spline we should apply an integration (look in Section 3.9 for more details and further applications of integration).

func.definite_integral(0, 6)
7.900961538461538

Here is another example:

v = [(i + RDF(i).sin()/2, i + RDF(i^2).cos()) for i in range(10)] s = spline(v) show(point(v) + plot(s,0,9, hue=.6))
Image in a Jupyter notebook

4.3. Limits of sequences and functions

Next we are interested in studying the convergence of given sequences or computing limits of functions. We begin with sequences.

Recall that a sequence is an enumerated collection of objects in which repetitions are allowed and order matters. Sequences can be viewed as functions NAR \mathbb{N}\to A\subseteq\mathbb{R}, where as usual N={1,2,3,}R\mathbb{N}=\{1, 2, 3, \ldots\}\subset\mathbb{R} denotes the set of natural numbers (we often include 00 as a natural), and AA is a subset of real numbers. We often denote a sequence by (an)nN(a_{n})_ {n\in\mathbb{N}}, where ana_{n} is its nnth-term. In the previous Chapter we met sequences defining by recursion, as the Fibonacci sequence, Fn=Fn1+Fn2. F_{n}=F_{n-1}+F_{n-2}\,.

However below we will treat more general sequences. An important property of a sequence is convergence, and this is what we are interested in below. If a sequence converges, it converges to a particular value known as the limit. If a sequence (an)nN(a_{n})_ {n\in\mathbb{N}}, converges to some limit, then it is called convergent, and moreover this limit is unique, denoted by limnan\lim_{n\to\infty}a_{n}. A sequence that does not converge is called divergent. An example of a divergent sequence is given by bn=(1)nb_n=(-1)^n.


Example (Limits of sequences)

Try to verify graphically{\textit{graphically}} via SageMath that the sequence an=(1)nna_{n}=\frac{(-1)^n}{n} converges to zero (for example, plot the first 100 terms of the sequece at hand).

Solution:

We have already seen in Chapter 1 how to plot sequences, via a method combining the Graphicsobject{\tt{Graphics object}} with the for{\tt{for}} command. For our case one can give the cell

A=Graphics() for i in srange (1, 100+1): A=A+points((i, (-1)^i/i)) A
Image in a Jupyter notebook

Hence (an)(a_n) should converge to 00. To comfirm this fact via Sage use the command limit{\tt{limit}} as follows (see also below)

var("n"); limit((-1)^n/n, n=oo) # oo is an alias of Infinity or infinity in Sage
0

Exercise

Prove that the seqence u(n):=2n2+4n+12n+1u(n):=\displaystyle\frac{2n^2+4n+1}{2n+1} tends to ++\infty as nn\to\infty.

Solution:

We need to show that

limnu(n)=+.\lim_{n\to\infty}u(n)=+\infty.

We may multiple both the numerators and denominator of u(n)u(n) by 1/n1/n. This gives the claim limn2n2+4n+12n+1=limn1n(2n2+4n+1)1n(2n+1)=limn2n+4+1n2+1n=+4+02+0=+,\begin{darray}{rcl} \lim_{n\to\infty}\frac{2n^2+4n+1}{2n+1}&=&\lim_{n\to\infty}\frac{\frac{1}{n}(2n^2+4n+1)}{\frac{1}{n}(2n+1)}\\\\ &=&\lim_{n\to\infty}\frac{2n+4+\frac{1}{n}}{2+\frac{1}{n}}\\\\ &=&\frac{\infty+4+0}{2+0}=+\infty\,, \end{darray} since 1n0\frac{1}{n}\to 0 as nn\to\infty. To verify this via Sage we may use the command limit{\tt{limit}} or its alias lim{\tt{lim}}, as follows (see also below)

n=var("n") u(n)=(2*n^2+4*n+1)/(2*n+1) lim(u(n), n=infinity)
+Infinity

Exercise

Prove that limna(n)=0\lim_{n\to\infty}a(n)=0, where a(n):=4n+12n2+4n+1a(n):=\displaystyle\frac{4n+1}{2n^2+4n+1}. Next plots its first 100 terms.

Solution:

n=var("n") a(n)=(4*n+1)/(2*n**2+4*n+1) lim(a(n), n=infinity)
0

Let us now plot the first 100 temrs of a(n)a(n) and illustrate the convergence to zero: Notice, on could type

plot(a(n), n, 1, 100)
Image in a Jupyter notebook

Recall however that a sequence is a function which has the natural numbers as domain, thus to get the correct graph of a(n)a(n) you should use the method presented above.

B=Graphics() for n in srange (1, 100+1): B=B+points((n, (4*n+1)/(2*n**2+4*n+1))) B
Image in a Jupyter notebook

Exercise

Compute the following limits in SageMath:

  1. limnA(n)\lim_{n\to\infty}A(n), where A(n):=(4n2+n)2nA(n):=\sqrt{(4n^2+n)}-2n.

  2. limnB(n)\lim_{n\to\infty}B(n), where B(n):=2n+2n2n2nB(n):=\displaystyle\frac{2^n+2^{-n}}{2^n-2^{-n}}.

  3. limnE(n)\lim_{n\to\infty}E(n), wherer E(n):=(1+1n)nE(n):=(1+\frac{1}{n})^n.

Solution:

The solutions can be obtained sing the command limit{\tt{limit}}, or its alias lim{\tt{lim}}.

n=var("n") A(n)=sqrt((4*n**2+n))-2*n lim(A(n), n=infinity)
1/4
n=var("n") B(n)=(2**n+2**(-n))/(2**n-2**(-n)) lim(B(n), n=infinity)
1
n=var("n") E(n)=(1+1/n)**n lim(E(n), n=infinity)
e

Hence limnE(n)=e\lim_{n\to\infty}E(n)=e, where ee is the Euler number (the base of the natural logarithm) (a classical result in calculus).

Limits of 1-variable functions.

The sequences are special examples of functions defined on the subset of natural numbers NR\mathbb{N}\subset\mathbb{R}. Next we treat limits of more general real-valued functions (with one variable).

In SageMath an appropriate command for computing limits of a given function f(x)f(x) is again based on the command limit{\tt{limit}} and has the form limit(f(x),x=a,direction){\tt{limit(f(x), x=a, direction)}} where x=aRx=a\in\mathbb{R}. The option ``direction'' is to study the limits to the left or to the right. This option can be omitted, but in case we need to use it can be introduced as dir="left"{\tt{dir="left"}} or dir="right"{\tt{dir="right"}}, or dir=""{\tt{dir="-"}} and dir="+"{\tt{dir="+"}}, respectively.

Remark.

By default, Sage calculates the limit to the left.

Example

Evaluate the limit limx04sin(x)x\lim_{x\to 0}\displaystyle\frac{4\sin(x)}{x}.

Solution

It is well-know that limx0sin(x)x=1\lim_{x\to 0}\displaystyle\frac{\sin(x)}{x}=1, so we expect that the given limit equals to 4. Indeed, in Sage we verify this as follows:

f(x)=4*sin(x)/x limit(f, x=0)
x |--> 4

Or just type

lim(4*sin(x)/x, x=0)
4

Let us also verify the claim by computing the left and right limits.

lim(4*sin(x)/x, x=0, dir="minus")
4
lim(4*sin(x)/x, x=0, dir="left") # as an alternative of the dir=minus is dir=left
4
lim(4*sin(x)/x, x=0, dir="plus")
4
lim(4*sin(x)/x, x=0, dir="right") # as an alternative of the dir=plus is dir=right
4

Since the left and right limit of f(x)=4sin(x)/xf(x)=4\sin(x)/x are both equal to 4, the claim follows.

Some remarks on the computation of limits

If f1,g1f_1, g_1 are functions on R\R with f1(x)0f_1(x)\ne 0 and g1(x)0g_1(x)\ne 0, for all xx, we can write f1(x)+f2(x)++fm(x)g1(x)+g2(x)++gn(x)=f1(x)(1+f2(x)f1(x)++fm(x)f1(x))g1(x)(1+g2(x)g1(x)++gn(x)g1(x)) \frac{f_1 (x) + f_2 (x) + \cdots + f_m (x)}{g_1 (x) + g_2 (x) + \cdots + g_n (x)} = \frac{f_1 (x)\bigl(1 + \frac{f_2 (x)}{f_1(x)} + \cdots + \frac{f_m(x)}{f_1(x)}\bigr)}{g_1 (x)\bigl(1 + \frac{g_2 (x)}{g_1(x)} + \cdots +\frac{g_n(x)}{g_1(x)}\bigr)} Thus, if limxx0fi(x)f1(x)=0\lim_{x\to x_0} \frac{f_i (x)}{f_1 (x)} = 0, for all i{2,,m}i \in \{2, \dots, m\}, and limxx0gi(x)g1(x)=0\lim_{x\to x_0} \frac{g_i (x)}{g_1 (x)}=0, for all i{2,,n}i \in \{2, \dots, n\}, then limxx0f1(x)+f2(x)++fm(x)g1(x)+g2(x)++gn(x)=limxx0f1(x)g1(x),\lim_{x\to x_0} \frac{f_1 (x) + f_2 (x) + \cdots + f_m (x)}{g_1 (x) + g_2 (x) + \cdots + g_n (x)}=\lim_{x\to x_0} \frac{f_1 (x)}{g_1 (x)}, supposing the limit on the right-hand side exists. This is an approach that we often apply to compute limits of functions given as fractions.

Squeeze Theorem

Another very useful result for the computation of limits is the so called Squeeze Theorem, which states the following:

Let IRI\subset\mathbb{R} be an interval containing the point x0x_0. Let g,f,hg, f, h be functions defined on II, except possibly at x0x_0 itslef. Suppose that for any xIx\in I we have that g(x)f(x)h(x),    and    limxx0g(x)=L=limxx0h(x). g(x)\leq f(x)\leq h(x), \ \ \ \ \text{and} \ \ \ \ \lim_{x\to x_0}g(x)=L=\lim_{x\to x_0}h(x). Then limxx0f(x)=L\lim_{x\to x_0}f(x)=L.

Exercise

Compute the following limits (note that some of them may not exist - specify them)

  1. limx2x+81\lim\frac{\sqrt{x}-2}{\sqrt{x+8}-1} when x0x\to 0, x3x\to 3, and x8x\to 8.

  2. limx8x32x+1933\lim_{x\to 8}\frac{\sqrt[3]{x}-2}{\sqrt[3]{x+19}-3}.

  3. limx0xx\lim_{x\to 0}\frac{x}{|x|}.

  4. limx0xx\lim_{x\to 0^-}\frac{x}{|x|}.

  5. limx0+xx\lim_{x\to 0^+}\frac{x}{|x|}.

  6. limx0sinxx\lim_{x\to 0}\frac{\sin x}{x} (and also the two-side limits, as above).

  7. limx0tanxx\lim_{x\to 0}\frac{\tan x}{x} (and also the two-side limits).

  8. limxsinxx\lim_{x\to\infty}\frac{\sin x}{x}.

  9. limx0sinxx\lim_{x\to 0}\frac{\sin x}{|x|} (and also the two-side limits).

  10. limx0sin(1x)\lim_{x\to 0}\sin(\frac{1}{x}).

  11. limx0xsin(1x)\lim_{x\to 0}x\sin(\frac{1}x).

  12. limx+ex\lim_{x\to+\infty}e^{-x}.

  13. limx+(1+1x)x\lim_{x\to+\infty}\left(1+\frac{1}{x}\right)^{x}.

  14. limx+lnxx\lim_{x\to+\infty}\frac{\ln x}{x}

  15. limx±xlnx\lim_{x\to\pm\infty}\frac{x}{\ln x}.

Solution:

  1. For the first case, one may apply the cell

b=lim((sqrt(x)-2)/(sqrt(x+8)-1), x=0) b.full_simplify ()
-2/(2*sqrt(2) - 1)

In case you like a decimal presentation, proceed as usual, i.e.,

N(b)
-1.09383632135605

For the case x3x\to 3 similarly we compute

N(lim((sqrt(x)-2)/(sqrt(x+8)-1), x=3))
-0.115663612660389

and for x8x\to 8 via the same method we obtain

N(lim((sqrt(x)-2)/(sqrt(x+8)-1), x=8))
0.276142374915397
  1. In the second case we get an indefinite of type 00\frac{0}{0}. Hence we can apply l'Hopital's rule (see below) to compute this limit by hand. By SageMath we get directly the result as follows:

lim((x**(1/3)-2)/((x+19)**(1/3)-3), x=8)
9/4

Let us present the rest limits in question.

lim(x/abs(x), x=0) # this limit doesn't exist
und

This answer means that the limit is undetermined, i.e., it doesn't exist. It may be useful to plot this function, to understand this result in a visual way.

plot(x/abs(x), x, -10, 10, exclude=[0])
Image in a Jupyter notebook
lim(x/abs(x), x=0, dir="left")
-1
lim(x/abs(x), x=0, dir="right")
1

Hence although the left and right limits exist, they are not equal and hece this gives another explanation

why the limit in case 3. does not exist.

  1. This is a very classical limit that everyone should know to prove also in a formal way (such one is based on the squeeze theorem).

lim(sin(x)/x, x=0, dir="right")
1
lim(sin(x)/x, x=0, dir="left")
1
lim(sin(x)/x, x=0)
1

7,

lim(tan(x)/x, x=0)
1
lim(tan(x)/x, x=0, dir="plus")
1
lim(tan(x)/x, x=0, dir="minus")
1
lim(sin(x)/x, x=oo)
0

We can illustrate this result by plotting the function f(x)=sinxxf(x)=\frac{\sin x}{x}, in the usual way

plot(sin(x)/x, x, -10*pi, 10*pi)
Image in a Jupyter notebook
lim(sin(x)/abs(x), x=0) # this limit doesnt exist
und
lim(sin(x)/abs(x), x=0, dir="+")
1
lim(sin(x)/abs(x), x=0, dir="-")
-1
limit((1+1/x)**x,x=infinity) # An indeterminate form!
e

We leave the rest few cases for practice.


4.4. Continuous functions

We now proceed with tasks on the continuity of functions. We recall first some basic definitions.

Definition. Let f:RRf: \mathbb{R} \to \mathbb{R} be a function, x0Rx_0 \in \mathbb{R}. ff is continuous at x0x_0 if limxx0f(x)=f(x0), \lim_{x \to x_0} f(x) = f(x_0), and ff is continuous on ARA \subset \mathbb{R} if for all aAa \in A, ff is continuous at aa.


Example

Consider the function f(x)f(x) defined as: f(x)={xif x<0x2if x0.f(x) = \begin{cases} -x & \text{if } x < 0 \\ x^2 & \text{if } x \ge 0. \end{cases} Show that ff is continuous.

Solution:

If a<0a < 0 then limxaf(x)=limxa(x)=a,\lim_{x \to a}f(x) = \lim_{x \to a}(-x) = -a, which coincides with f(a)f(a). Same argument applies when a>0a > 0. It remains to show that limx0f(x)=f(0)\lim_{x\to 0}f(x) = f(0).

(-x).limit(x=0, dir='-')
0
(x^2).limit(x=0, dir='+')
0

One-sided limits at 00 coincide: they are both equal to 00. Combining this with the relation f(0)=0f(0)=0,we deduce that ff is also continuous at 00.

Exercise

Consider the function f(x,a,b)f(x, a, b) defined as: f(x,a,b)={ax+bif x<0x2if x0.f(x, a, b) = \begin{cases} ax + b & \text{if } x < 0 \\ x^2 & \text{if }x \ge 0. \end{cases} Specify the values of the parameters aa and bb such that f(x,a,b)f(x, a, b) is continuous at x=0x = 0.

Solution:

var('a') var('b') # we already know that the right limit of x^2 at x=0 is 0 (a*x+b).limit(x=0, dir='-') == 0
b == 0

Exercise for practice

Consider the piecewise function defined below directly in Sage. Examine its continuity at x=0x=0.

f1(x)=2*e^x; f2(x)=x^2+2*x+2 f=piecewise([[RealSet.open_closed(0, 5), f1(x)], [RealSet.closed(-5, 0), f2(x)]]) p=plot(f(x), x, -5, 5) show(p)
Image in a Jupyter notebook

Exercise

Consider the function g(x,a)g(x, a) defined as: g(x,a)={x3+axif x<13x1if x1.g(x, a) = \begin{cases} x^3 + ax & \text{if } x < 1 \\ 3x - 1 & \text{if }x \ge 1. \end{cases} Specify the values of the parameter aa such that g(x,a)g(x, a) is continuous at x=1x = 1.

Solution:

We present the solution directly in Sage, without any comments (try to confirm the mathematical computation behind the code below)

var('a') eq = (x^3 + a*x).limit(x=1, dir='-') == (3*x-1).limit(x=1, dir='+') print(eq) eq.solve(a)
a + 1 == 2
[a == 1]

Intermediate Value Theorem

Let f:[a;b]Rf : [a; b] \to \mathbb{R} be a continuous function. Then ff takes any given value between f(a)f(a) and f(b)f(b) at some point within [a;b][a;b].

Corollary

If f(a)f(b)<0f(a) \cdot f(b) < 0 then there is x[a;b]x \in [a;b] such that f(x)=0f(x) = 0.

Bisection method

There is a simple algorithm to find an approximation of this xx that is called bisection method.

def bisect_method(f, a, b, eps): '''Bisection method. Returns c, (a,b), where c is the midpoint of the interval [a;b], and it holds |f(c)| < eps, f(a)*f(b) < 0. ''' eps = float(eps) while True: c = (a+b)/float(2) fa = f(a); fb = f(b); fc = f(c) if abs(fc) < eps: return c, (a,b) if fa*fc < 0: b = c elif fc*fb < 0: a = c else: raise ValueError("f must have a sign change in the interval (%s,%s)"%(a,b))

Example

Find 33\sqrt[3]{3} by approximating the root of f(x)=3x3f(x) = 3 - x^3 on the interval [1;2][1;2] via the method described above.

f(x) = 3 - x^3 (r, _) = bisect_method(f, 1, 2, 0.0000001) r
1.4422495663166046
r^3
2.999999975096381

4.5. Derivatives

We proceed with tasks on derivatives. We first recall some theoretical details on derivatives.

Theoretical notes on derivatives

Let ff be a real or complex function defined on an interval ARA ⊂ R and x0Ax_0 ∈ A. If the limit limxx0f(x)f(x0)xx0=a\lim_{x\to x_0}\frac{f(x) − f(x_0)}{x − x_0}= a exists, the function ff is said to be differentiable at x0x_0, provided aa is finite. The value of the derivative at x0x_0, namely aa, is denoted by f(x0)f'(x_0) or dfdx(x0)\frac{df}{dx} (x_0) or ddxf(x0).\frac{d}{dx} f(x_0). If aa is finite, the derivative is also sometimes called proper. If aa is infinite, it is improper. If x0x_0 is one of the boundary points of AA, we arrive at one-sided derivatives (i.e., left-sided derivative and right-sided derivative). If a function has a derivative at x0x_0, the function is said to be differentiable at x0x_0. A function which is differentiable at every point of a given interval is said to be differentiable on the interval.

Rules of differentiation:

  1. (cf)(x)=c(f(x))(cf)'(x) = c(f'(x)), for every real or complex number cc;

  2. (f+g)(x)=f(x)+g(x)(f+g)'(x) = f'(x)+g'(x)

  3. (fg)(x)=f(x)g(x)+g(x)f(x)(f\cdot g)'(x) = f'(x)g(x)+g'(x)f(x)

  4. (fg)(x)=f(x)g(x)g(x)f(x)g2(x)(\frac{f}{g})'(x) = \frac{f'(x)g(x)-g'(x)f(x)}{g^2(x)}

  5. h(f(x))=h(f(x))f(x)h(f(x))' = h(f(x))'\cdot f'(x)

We can solve derivatives in SageMath in at least three different ways:

  1. by using limits via the command limit{\tt{limit}};

  2. by using diff(){\tt{diff()}} or derivative(){\tt{derivative()}} functions;

  3. by using SAGE function method .derivative(){\tt{.derivative()}}.

We will explore the options that we can apply for any case via examples.

Exercise

Find the derivative of f(x)=x3f(x) = x^3 in SageMath by apply the difinition of ff', and next use the diff{\tt{diff}} function to verify the result of the first method.

Solution:

Let us first determine the derivative in question by using the definition and hence limits. Recall that f(x)f'(x) is defined as the limit f(x):=limh0f(x+h)f(x)hf'(x):=\lim_{h\to0} \frac{f(x+h)-f(x)}{h} and we say that f(x)f'(x) exists when this limit exists. In Sage this goes as follows:

f(x)=x^3 # define our function var('h') # create variable h limit((f(x+h)-f(x))/h, h=0) # find the limit
3*x^2

Great! We received our result, that is, f(x)=(x3)=3x2f'(x)=(x^3)'=3x^2. Let us now shorten the code with the built-in function diff(){\tt{diff()}}, as suggested in the statement.

diff(x^3,x) # first argument is the function and second one is variable to differentiate
3*x^2

Note that the same can be done with the derivative(){\tt{derivative()}} function, as follows:

derivative(x^3,x)
3*x^2

Finally, let us also use a third variant that Sage offers for the computation of the derivative of a function, based on the method derivative(){\tt{derivative()}}. This is essentially the same as above, but is written in the following form:

f(x)=x^3 f.derivative(x)
x |--> 3*x^2

This second version of the syntax is more like Python, and makes it clear that the derivative is an attribute, or method associated with the object f(x). In particular, one could also type f(x).derivative(x){\tt{f(x).derivative(x)}}, as it is shown below:

f(x).derivative(x)
3*x^2

Exercise (evaluation of derivatives at certain points)

Find the particular value of derivative of f(x)=x3f(x)=x^3 for f(3)f'(3) Here we can also use two methods

Hints:{\color{red} \textsf{Hints:}}

  1. Create a function dfdf and just input a value

  2. We can specify the variable after diff()(x=...){\tt{diff()(x=...)}}, see below.

Solution:

Let us apply first the first method:

f(x)=x^3 df(x) = f.derivative(x) # or diff(f) or derivative(f) df(3)
27

Let us also specify the details of another method, among those mentioned above, that one can apply to get the same answer:

f.derivative(x)(x=3) # or diff(f)(x=3) or derivative(f)(x=3)
27

Or we can type

f(x).derivative(x)(x=3)
27

In general, we prefer more the second method.

Example

Find the derivative of the function g(x)=5x3e(2x2)g(x)=5x^3\,e^(-2x^2) for all real numbers xx and next find g(1)g'(1). Moreover, present a numerical approxiamtion of g(1)g'(1).

g(x)=5*x^3*e^(-2*x^2) dg(x)=derivative(g(x),x) dg(1)
-5*e^(-2)

For the numerical approximation we can use the N(){\tt{N()}} function:

N(dg(1))
-0.676676416183064

Similarly for another value:

N(dg(2))
-0.0872202832546531

With Sage we can also handle derivatives when other variables are involved. In this case, all variables other than the variable of differentiation (the variable after the comma in the derivative command) are treated as constants (unspecified numbers). However, one should still introduce these varianbles as symbolic variables in Sage.

Exercise

Find the derivative of the function f(x)=ax4+bx3+cx+edx\displaystyle f(x)=ax^4+\sqrt{b}x^3+\frac{c}{x}+e^{dx}, where a,b,c,da, b, c, d are positive constants and x(0,+)x\in(0, +\infty). Next evaluate the derivative at x=1x=1

Solution:

var("a, b, c, d, x") f(x)=a*x^4+sqrt(b)*x^3+(c/x)+e^(d*x) df=diff(f(x), x) print("The derivative of f is given by", df)
The derivative of f is given by 4*a*x^3 + 3*sqrt(b)*x^2 + d*e^(d*x) - c/x^2

To compute the value f(1)f'(1) we can add the cell:

print("The deivative of f at x=1 is given by", diff(f(x), x)(x=1))
The deivative of f at x=1 is given by d*e^d + 4*a - c + 3*sqrt(b)

Application of first order derivatives (velocity of an object)

Recall that the velocity of a moving object is the derivative of its position function, and its acceleration is the derivative of its velocity function. Next we will use the units m/secm/sec for velocity and m/sec2m/sec^2 for acceleration (though one could instead use km/hkm/h for velocity and respectively for acceleration).

Exercise

If the position of a moving object in time t is given by the function s(t)=(2t3)29, s(t)=(2t-3)^2-9, determine the velocity and the acceleration of the object for general time tt and next for time t=3t=3.

Solution
var("t") s(t)=(2*t-3)^2-9 #the position function u(t)=diff(s(t), t).factor() #the velocity print("The velocity of the object at general time t is the function given by", u(t), "m/sec.")
The velocity of the object at general time t is the function given by 8*t - 12 m/sec.
var("t") s(t)=(2*t-3)^2-9 u(t)=diff(s(t), t).factor() a(t)=diff(u(t), t) #the acceleration print("The acceleration of the object at general time t is the derivative of the velocity, given by", a(t), "m/sec^2.") print("The velocity of the object at time t=3 is given by", u(3), "m/sec.")
The acceleration of the object at general time t is the derivative of the velocity, given by 8 m/sec^2. The velocity of the object at time t=3 is given by 12 m/sec.

Hence the object is moving with constant acceleratiob 8m/sec^2 (assuming that we compute the velocity in m/sec).

Higher order derivatives

The previous problem it could solved by a different approach, relying on second order derivatives (as the acceleration is the second derivatibe of the position function). Usually, given a function f:RRf : \mathbb{R}\to\mathbb{R} we denote the second derivative by f(x)=ddff(x)f''(x)=\frac{d}{df}f'(x). In general for higher order derivatives we use the notation f(3)=ddxf(x),,f(n)(x)=ddxf(n1)(x). f^{(3)}=\frac{d}{dx}f''(x)\,, \quad \ldots, \quad f^{(n)}(x)=\frac{d}{dx}f^{(n-1)}(x)\,. We call f(n)(x)f^{(n)}(x) the nnth derivative of ff at xx. See the 6th Chapter of the BG book for the formal definition of higher-order derivatives as limits.

Exercise

Find the 2nd order, 3rd order and 5th order derivatives of the function f(x)=x5e5f(x)=x^5e^5 with xRx\in\mathbb{R}. Next evaluate the second-order derivative at x=0x=0, at x=1x=1 and at x=1x=-1.

Solution:

One can use our previous functions and methods diff(){\tt{diff()}}, derivative(){\tt{derivative()}} or, f.derivative(){\tt{f.derivative()}} and add a third argument (or second if the variable is not specified), which will show the order of the derivative that need to compute.

f(x)=x^5*e^5 # define function diff(f, x, 2) #or derivative(f,x,2)) or f(x).derivative(2) diff(f,x,3) # or derivative(f,x,3)) or f(x).derivative(3) diff(f,x,5) # or derivative(f,x,3)) or f(x).derivative(3) show("The second order derivative of f is given by:", diff(f,x,2)) # for a fancy output we use show() function show("The third order derivative of f is given by:", diff(f,x,3)) show("The fifth order derivative of f is given by:", diff(f,x,5))

The second order derivative of f is given by:x  20x3e5\displaystyle \verb|The|\verb| |\verb|second|\verb| |\verb|order|\verb| |\verb|derivative|\verb| |\verb|of|\verb| |\verb|f|\verb| |\verb|is|\verb| |\verb|given|\verb| |\verb|by:| x \ {\mapsto}\ 20 \, x^{3} e^{5}

The third order derivative of f is given by:x  60x2e5\displaystyle \verb|The|\verb| |\verb|third|\verb| |\verb|order|\verb| |\verb|derivative|\verb| |\verb|of|\verb| |\verb|f|\verb| |\verb|is|\verb| |\verb|given|\verb| |\verb|by:| x \ {\mapsto}\ 60 \, x^{2} e^{5}

The fifth order derivative of f is given by:x  120e5\displaystyle \verb|The|\verb| |\verb|fifth|\verb| |\verb|order|\verb| |\verb|derivative|\verb| |\verb|of|\verb| |\verb|f|\verb| |\verb|is|\verb| |\verb|given|\verb| |\verb|by:| x \ {\mapsto}\ 120 \, e^{5}

To evaluate a higher-order derivative at a certain point we can use the same approach as for the evaluation of the derivative of a function ff.

f(x)=x^5*e^5 diff(f, x, 2)(x=0)
0

Or we can type:

print("The second order derivative of f at x=0 equals to", diff(f, x, 2)(x=0))
The second order derivative of f at x=0 equals to 0
print("The second order derivative of f at x=1 equals to", diff(f, x, 2)(x=1))
The second order derivative of f at x=1 equals to 20*e^5

In case we need a numerical approximation of f(1)f''(1) we can add:

print("The second order derivative of f at x=1 equals to", N(diff(f, x, 2)(x=1)))
The second order derivative of f at x=1 equals to 2968.26318205153

Similarly we can compute f(1)f''(-1):

print("The second order derivative of f at x=-1 equals to", N(diff(f, x, 2)(x=-1)))
The second order derivative of f at x=-1 equals to -2968.26318205153

Exercise

Find the acceleration of an object at time general time tt and at time t=1/2t=1/2 sec, when its position function is given by s(t)=t(1+t2)+t4s(t)=-\frac{t}{(1+t^2)}+t^4, with t0t\geq 0.

Solution:

var("t") s(t)=-(t/(1+t^2))+t^4 a(t)=diff(s, t, 2) show(a(t))

12t28t3(t2+1)3+6t(t2+1)2\displaystyle 12 \, t^{2} - \frac{8 \, t^{3}}{{\left(t^{2} + 1\right)}^{3}} + \frac{6 \, t}{{\left(t^{2} + 1\right)}^{2}}

print("The acceleration of the object at time t=1/2 is", a(1/2), "m/sec^2.")
The acceleration of the object at time t=1/2 is 551/125 m/sec^2.

Or we could directly type:

var("t") s(t)=-(t/(1+t^2))+t^4 print("The acceleration of the object at time t=1/2 is", diff(s, t, 2)(t=1/2), "m/sec^2.")
The acceleration of the object at time t=1/2 is 551/125 m/sec^2.

A few details on partial derivatives

Partial derivatives are derivatives that can be used when you have a function of multiple variables. Unlike ordinary derivatives, which describe how a function changes with respect to one variable, partial derivatives measure how a function changes with respect to one variable while keeping all other variables constant. For instance, consider the function f(x,y)=x2y+3xy2f(x, y)=x^2y+3xy^2. Then we have the following two first-order partial derivatives:

fx=differentiation with respect to x=2xy+3y2.\frac{\partial f}{\partial x}=\text{differentiation with respect to x}=2xy+3y^2.fy=differentiation with respect to y=x2+6xy.\frac{\partial f}{\partial y}=\text{differentiation with respect to y}=x^2+6xy.

For more details see Chapter 8 in the BG book.

Exercise

Consider the function f(x,y)=x5exy5f(x,y)=x^5e^xy^5. Find the 1st and 3rd order partial derivatives with respect to x,yx, y.

Solution:
var('y') # declare y variable f(x,y) = x^5*e^x*y^5 # define function show('1st order diff by x:', diff(f,x)) # differentiate by x show('1st order diff by y:', diff(f,y)) # differentiate by y print('\n') show('3rd order diff by x:', diff(f,x,3)) # differentiate by x show('3rd order diff by y:', diff(f,y,3)) # differentiate by y print('\n')

1st order diff by x:(x,y)  x5y5ex+5x4y5ex\displaystyle \verb|1st|\verb| |\verb|order|\verb| |\verb|diff|\verb| |\verb|by|\verb| |\verb|x:| \left( x, y \right) \ {\mapsto} \ x^{5} y^{5} e^{x} + 5 \, x^{4} y^{5} e^{x}

1st order diff by y:(x,y)  5x5y4ex\displaystyle \verb|1st|\verb| |\verb|order|\verb| |\verb|diff|\verb| |\verb|by|\verb| |\verb|y:| \left( x, y \right) \ {\mapsto} \ 5 \, x^{5} y^{4} e^{x}

3rd order diff by x:(x,y)  x5y5ex+15x4y5ex+60x3y5ex+60x2y5ex\displaystyle \verb|3rd|\verb| |\verb|order|\verb| |\verb|diff|\verb| |\verb|by|\verb| |\verb|x:| \left( x, y \right) \ {\mapsto} \ x^{5} y^{5} e^{x} + 15 \, x^{4} y^{5} e^{x} + 60 \, x^{3} y^{5} e^{x} + 60 \, x^{2} y^{5} e^{x}

3rd order diff by y:(x,y)  60x5y2ex\displaystyle \verb|3rd|\verb| |\verb|order|\verb| |\verb|diff|\verb| |\verb|by|\verb| |\verb|y:| \left( x, y \right) \ {\mapsto} \ 60 \, x^{5} y^{2} e^{x}

Exercise

Find the first order and second order partial derivatives of f(x,y)=x5y5f(x, y)=x^5y^5.

show(diff(x^5*y^5, x)) show(diff(x^5*y^5, y)) show(diff(x^5*y^5,[x,y])) show(diff(x^5*y^5,[y,x]))

5x4y5\displaystyle 5 \, x^{4} y^{5}

5x5y4\displaystyle 5 \, x^{5} y^{4}

25x4y4\displaystyle 25 \, x^{4} y^{4}

25x4y4\displaystyle 25 \, x^{4} y^{4}

Geometrical meaning of the derivative

Derivatives allow us to get instantaneous rate of change or slope of a function at a given point. More specifically, the derivative at a point on a curve represents the slope of the tangent line to the curve at that point. Geometrically, the tangent line touches the curve at only one point and is the best linear approximation of the curve near that point. The derivative gives us information about the direction and steepness of the curve at any point, allowing us to analyze the behavior of the curve and make predictions about its future behavior. In essence, the derivative provides us with a powerful tool for understanding and analyzing the geometry of functions.

The geometric meaning of the derivative allows to approximate

f(x+h)=f(x)+ϕ(x+h)hf(x)+f(x)hf(x + h) = f(x) + \phi(x + h)h ≃ f(x) + f′(x)h

A real function is called increasing at x0x_0 of its domain, if for all points x of some neighbourhood of a point x0,f(x)>f(x0)x_0, f(x) > f(x_0) if x>x0x > x_0 and f(x)<f(x0)f(x) < f(x_0) if x<x0x < x_0. A real function is increasing on an interval AA if f(x)f(y)>0f(x)−f(y) > 0 for all x>y,x,yAx > y, x, y \in A. Similarly, a function is said to be decreasing at a point x0x_0 if and and only if there is a neighbourhood of the point x0x_0 such that f(x)<f(x0)f(x) < f(x_0) for all x>x0x > x_0, while f(x)>f(x0)f(x) > f(x_0) for all x<x0x < x_0 from this neighbourhood. A function is decreasing on an interval A if f(x)f(y)<0f(x) − f(y) < 0 for all x>y,x,yAx > y, x, y \in A. Thus a function having a non-zero finite derivative at a point is either increasing or decreasing at that point, according to the sign of the derivative.

Exercise

Find the derivative of f(x)=x2sin(x)f(x)=x^2sin(x), and plot the tangent lines for x0=3x_0=3 and x1=2.3x_1=-2.3.

Solution:

f(x) = x**2*sin(x) df(x) = f.differentiate() x0 = 3 x1 = -2.3 tang_line_1 = f(x0) + df(x0)*(x-x0) tang_line_2 = f(x1) + df(x1)*(x-x1) (plot(f, xmin=-4, xmax=5, color='black') + # plot original function plot(tang_line_1, xmin=1, xmax=5, color='red') + # plot first tangent line plot(tang_line_2, xmin=-5, xmax=0.5, color='red') + # plot first tangent line point((x0,f(x0)), color='blue', size=20) + # plot first point of intersection point((x1,f(x1)), color='blue', size=20)) # plot second point of intersection
Image in a Jupyter notebook

Exercise

Find the local extrema for f(x)=5x3+2x23xf(x)=5x^3+2x^2-3x. Find out whether there any maxima or minima.

Solution:

f(x) = 5*x**3 + 2*x**2 - 3*x - 5#x**2+(2*x**2)+3 # First order derivative df(x) = f.differentiate() # Solve for df(x)==0 results = solve(df(x)==0, x, domain='real') show(results) # Plot the function plots = plot(f, xmin=-1.5, xmax=1.5, color='black') for result in results: # Find second order derivative d2f = df.diff() # or f.diff(x,2) # Find out whether the point is maximum or minumum print("\nSecond order derivative of f(x) in ({0}) =".format(result), d2f(result.right())) if d2f(result.right()) > 0: print(result, 'is a local minimum') elif d2f(result.right()) < 0: print(result, 'is a local maximum') else: print(result,'needs more analysis') # Plot tangent lines and points tang_line = f(result.right()) + df(result.right())*(x-result.right()) plots += plot(tang_line, color='blue') + point((result.right(), f(result.right())), color='red', size=40) plots

[x=(13),x=(35)]\displaystyle \left[x = \left(\frac{1}{3}\right), x = \left(-\frac{3}{5}\right)\right]

Second order derivative of f(x) in (x == (1/3)) = 14 x == (1/3) is a local minimum Second order derivative of f(x) in (x == (-3/5)) = -14 x == (-3/5) is a local maximum
Image in a Jupyter notebook

Exercise (tangent lines)

Find the tangent line to the function g(x)=x4exg(x)=x^4\, e^x for x=2x=2. Next create an illustration including the graph of gg and the tangent line at question.

Solution:

g(x)=x^4*e^(x) dg(x)=derivative(g(x),x) tangent(x)=g(2)+dg(2)*(x-2) show("The tangent line of the function g(x) at x=2 is the line given by:", tangent(x)) plot(g,xmin=0,xmax=4)+plot(tangent,xmin=0,xmax=4,color='red')+point((2,g(2)),color='black',size=25)

The tangent line of the function g(x) at x=2 is the line given by:48(x2)e2+16e2\displaystyle \verb|The|\verb| |\verb|tangent|\verb| |\verb|line|\verb| |\verb|of|\verb| |\verb|the|\verb| |\verb|function|\verb| |\verb|g(x)|\verb| |\verb|at|\verb| |\verb|x=2|\verb| |\verb|is|\verb| |\verb|the|\verb| |\verb|line|\verb| |\verb|given|\verb| |\verb|by:| 48 \, {\left(x - 2\right)} e^{2} + 16 \, e^{2}

Image in a Jupyter notebook

Exercise

Find the second derivative of h(x)=sin(x)+cos(x)h(x) = sin(x) + cos(x)

Solution:

Let us use Sage's symbolic calculus capabilities to find the first and second derivatives of the given function:

# Define the function h(x) h(x) = sin(x) + cos(x) # Find the first derivative of h(x) dh = diff(h(x), x) # Find the second derivative of h(x) d2h = diff(dh, x) # Simplify the second derivative using algebraic manipulation show(d2h.simplify_full())

cos(x)sin(x)\displaystyle -\cos\left(x\right) - \sin\left(x\right)

Derivatives of polynomial functions

Polynomials are some of the easiest and most commonly used class of functions that we encounter in Mathematics. They can easily and efficiently be evaluated at every specific value and have very nice properties: they are continuous everywhere, differentiable everywhere, and integrable on any bounded interval. Finding the derivatives or antiderivatives of a polynomial are very easy tasks too.

Exercise

Find the derivative of the polynomial function given below at x=2x=2: f(x)=x34x2+2x7.f(x) = x^3 - 4x^2 + 2x - 7.

Solution:

We can use SAGE's symbolic calculus capabilities to find the derivative of the function and then evaluate it at x = 2.

# Define the function f(x) f(x) = x^3 - 4*x^2 + 2*x - 7 # Find the derivative of f(x) df = diff(f(x), x) # Evaluate the derivative at x = 2 df.subs(x=2)
-2

Exercise

Find the derivative of g(x)=(x2+3x+1)/(x+2)g(x) = (x^2 + 3x + 1)/(x + 2)

Solution:

We can use the quotient rule to find the derivative of g(x) and simplify the expression.

# Define the function g(x) g(x) = (x^2 + 3*x + 1)/(x + 2) # Find the derivative of g(x) dg = diff(g(x), x) # Simplify the derivative using algebraic manipulation show(dg.simplify_full())

x2+4x+5x2+4x+4\displaystyle \frac{x^{2} + 4 \, x + 5}{x^{2} + 4 \, x + 4}

4.6. Elementary theorems on derivatives

We now proceed with some elementary properties of derivatives and in particular some basic theorems related to them.

Derivative of the inverse function

If ff is a real-valued function differentiable at y0y_0, such that the inverse f1(x)f^{−1}(x) exists on a neighbourhood of the value x0=f(y0)x_0 = f(y_0) and f(y0)0f′(y_0) \neq0, then (f1)(x0)=1f(f1(x0))=1f(y0)(f^{-1})′(x_0) = \frac{1}{f′(f^{-1}(x_0))} = \frac{1}{f′(y_0)}

Rolle’s theorem

Assume that the function f:RRf : R → R is continuous on a closed bounded interval [a, b] and differentiable inside this interval. If f(a)=f(b)f(a) = f(b), then there is a number c(a,b)c ∈ (a, b) such that f(c)=0f′(c) = 0

Lagrange’s mean value theorem

Assume the function f:RRf : R → R is continuous on an interval [a, b] and differentiable at all points inside this interval. Then there is a number c(a,b)c ∈ (a, b) such that f(c)=(f(b)f(a))/(ba)f′(c) = (f(b) − f(a))/(b − a)

Cauchy’s mean value theorem

Let functions y=f(t)y = f(t) and x=g(t)x = g(t) be continuous on an interval [a, b] and differentiable inside this interval. Further, let g(t)0g′(t) \neq 0 for all t(a,b)t ∈ (a, b) and g(b)g(a)g(b) \neq g(a). Then there is a point c(a,b)c ∈ (a, b) such that f(b)f(a)g(b)g(a)=f(c)g(c)\frac{f(b)-f(a)}{g(b)-g(a)} = \frac{f′(c)}{g′(c)}

The l'Hopital's rule

Finally let us recall the l'Hopital's rule, which is related to differentiation.

Assume that ff and gg are differentiable functions on an open interval IRI\subset\mathbb{R}, except possibly at a point cIc\in I. Suppose that limxcf(x)=limxcg(x)=0\lim_{x\to c}f(x)=\lim_{x\to c}g(x)=0, or limxcf(x)=limxcg(x)=±\lim_{x\to c}f(x)=\lim_{x\to c}g(x)=\pm\infty, and g(x)0g'(x)\neq 0 for all xIx\in I with xcx\neq c. Moreover, assume that the limit limxcf(x)g(x)\lim_{x\to c}\frac{f'(x)}{g'(x)} exists. Then the l'Hopital's rule states that limxcf(x)g(x)=limxcf(x)g(x). \lim_{x\to c}\frac{f(x)}{g(x)}=\lim_{x\to c}\frac{f'(x)}{g'(x)}\,.


Exercise (critical points)

Find the critical points of f(x)=x214x2+8f(x)=\frac{x^2-1}{4x^2+8}. Next plot ff for 2x2-2\leq x\leq 2.

Solution:

# Define the variable x = var('x') # Define the P(x) f = (x^2+1)/(4*x^2+8) show(f.factor()) # Step 1: Plot the polynomial to visualize its behavior plotf=f.plot(-2, 2) show(plotf) # Step 2: Compute the first derivative of f(x) f_prime = diff(f, x) show(f_prime) # Step 3: Solve f'(x) = 0 to find critical points (local maxima and minima) critical_points = solve(f_prime == 0, x) show(critical_points)

x2+14(x2+2)\displaystyle \frac{x^{2} + 1}{4 \, {\left(x^{2} + 2\right)}}

Image in a Jupyter notebook

x2(x2+2)(x2+1)x2(x2+2)2\displaystyle \frac{x}{2 \, {\left(x^{2} + 2\right)}} - \frac{{\left(x^{2} + 1\right)} x}{2 \, {\left(x^{2} + 2\right)}^{2}}

[x=0]\displaystyle \left[x = 0\right]

Exercise (critical points)

Find the critical points of the function P(x)=4x33x+2P(x)=4x^3-3x+2. Next plot in one figure the given function PP and its (first) derivative and mark the c Use green color for the graph of the derivative.

Solution:

# Define the variable and the function f(x) x = var('x') f = 4*x^3 - 3*x + 2 # Compute the derivative of f(x) f_prime = diff(f, x) show("The derivative of f is given by:", f_prime) criticalp = solve(f_prime== 0, x) show("The critical points of f are the points:", criticalp) # Plot f(x) and f'(x) p1 = plot(f, (x, -2, 2), color='blue', legend_label=r"$f(x)$") p2 = plot(f_prime, (x, -2, 2), color='green', legend_label=r"$f'(x)$") # Mark the critical points critical_points = [(-1/2, f.subs(x=-1/2)), (1/2, f.subs(x=1/2))] p3 = point(critical_points, color='red', size=30, legend_label="Critical points") # Combine the plots combined_plot = p1 + p2 + p3 # Show the plot with labels and legend combined_plot.show()

The derivative of f is given by:12x23\displaystyle \verb|The|\verb| |\verb|derivative|\verb| |\verb|of|\verb| |\verb|f|\verb| |\verb|is|\verb| |\verb|given|\verb| |\verb|by:| 12 \, x^{2} - 3

The critical points of f are the points:[x=(12),x=(12)]\displaystyle \verb|The|\verb| |\verb|critical|\verb| |\verb|points|\verb| |\verb|of|\verb| |\verb|f|\verb| |\verb|are|\verb| |\verb|the|\verb| |\verb|points:| \left[x = \left(-\frac{1}{2}\right), x = \left(\frac{1}{2}\right)\right]

Image in a Jupyter notebook

Note that as the 2 critical points of ff are solutions of the equation f(x)=0f'(x)=0, we can attain them as follows:

xa = criticalp [0].rhs() show("The first critical point of f is the point:", xa) xb = criticalp [1].rhs() show("The second critical point of f is the point:", xb)

The first critical point of f is the point:12\displaystyle \verb|The|\verb| |\verb|first|\verb| |\verb|critical|\verb| |\verb|point|\verb| |\verb|of|\verb| |\verb|f|\verb| |\verb|is|\verb| |\verb|the|\verb| |\verb|point:| -\frac{1}{2}

The second critical point of f is the point:12\displaystyle \verb|The|\verb| |\verb|second|\verb| |\verb|critical|\verb| |\verb|point|\verb| |\verb|of|\verb| |\verb|f|\verb| |\verb|is|\verb| |\verb|the|\verb| |\verb|point:| \frac{1}{2}

This gives us an alternative to rewrite the program for marking these critical points as follows

# Plot f(x) and f'(x) p1 = plot(f, (x, -2, 2), color='blue', legend_label=r"$f(x)$") p2 = plot(f_prime, (x, -2, 2), color='green', legend_label=r"$f'(x)$") # Mark the critical points critical_points = [(xa, f.subs(x=xa)), (xb, f.subs(x=xb))] p3 = point(critical_points, color='red', size=30, legend_label="Critical points") # Combine the plots combined_plot = p1 + p2 + p3 # Show the plot with labels and legend combined_plot.show()
Image in a Jupyter notebook

Exercise (Rolle's theorem)

Verify Rolle's theorem for the function f(x)=x24x+3f(x)=x^2-4x+3 in the interval [1,3][1, 3].

Solution

# Define the function f(x) x = var('x') f = x^2 - 4*x + 3 plotf=plot(f, 0, 4, legend_label=r"$f(x)$") # Check if f(1) = f(3) f_at_1 = f.subs(x=1) show(f_at_1) f_at_3 = f.subs(x=3) show(f_at_3) # Derivative of f(x) f_prime = diff(f, x) plotf_prime=plot(f_prime, 0, 4, color="darkblue", thickness=2, legend_label=r"$f'(x)$") show(plotf+plotf_prime) # Find c such that f'(c) = 0 by solving f'(x) = 0 c_value = solve(f_prime == 0, x) show(c_value)

0\displaystyle 0

0\displaystyle 0

Image in a Jupyter notebook

[x=2]\displaystyle \left[x = 2\right]

4.7. Series

Recall that a series is the sum of the terms of a sequence, usually denoted by

n=1an,\sum_{n=1}^{\infty} a_n,

where ana_n is the nnth term of the sequence. If the sequence Sn=k=1nakS_n = \sum_{k=1}^n a_k converges to a limit LL we write n=1an=L\sum_{n=1}^{\infty} a_n = L. If the sequence SnS_n does not converge to a limit, then the series n=1an\sum_{n=1}^{\infty} a_n is said to diverge.

A geometric series is a series of the form n=0arn\sum_{n=0}^{\infty} ar^n, where aa and rr are constants. The sum of a geometric series is given by:

n=0arn=a1r, if r<1.\sum_{n=0}^\infty a r^n = \frac{a}{1-r}, \text{ if } |r| < 1.

A telescoping series is a series in which most of the terms cancel out:

n=1(anan1)=limnana0\sum_{n=1}^\infty (a_n - a_{n-1}) = \lim_{n \to \infty} a_n - a_0

The sum of a series can be computed in SageMath by using the sum(){\tt{sum()}} function. The syntax is sum(f,x,a,b){\tt{sum(f, x, a, b)}}, where f is the function to be summed, x is the variable of summation, a is the starting index of summation, and b is the ending index of summation.

Example

Consider a series k=0xk, \sum_{k=0}^\infty x^k, that converges if x<1|x| < 1 and is equal to 11x. \frac{1}{1 - x}.

var('x') var('n') assume(abs(x) < 1) print(sum(x^n, n, 0, oo))
-1/(x - 1)

Let us now compute a partial sum

reset('x') var('x') var('m') # diverges if abs(x) >= 1 print(sum(x^n, n, 0, m))
(x^(m + 1) - 1)/(x - 1)

Exercise

Use SageMath to compute the sum of the series n=11n2\sum_{n=1}^{\infty} \frac{1}{n^2}.

Solution:

var('n') sum(1/n^2, n, 1, oo)
1/6*pi^2

Exercise

Compute the sum of the series n=0(1)n2n+1\sum_{n=0}^{\infty} \frac{(-1)^n}{2n+1}.

Solution:

var("n") sum((-1)^n/(2*n+1), n, 0, oo)
1/4*pi

Exercise

Find the value of n=11n(n+1)\sum_{n=1}^{\infty} \frac{1}{n(n+1)} using partial fraction decomposition. Verify in this way that the sum converges and compare the result to the value of the telescoping sum n=1(1n1n+1)\sum_{n=1}^{\infty} \left(\frac{1}{n}-\frac{1}{n+1}\right).

Solution:

var("n") sum(1/(n*(n+1)), n, 1, oo)
1

On the other hand we see that:

sum(1/n-1/(n+1), n, 1, oo) #hence the claim
1

Exercise

Compute the sum of the series n=0xnn!\sum_{n=0}^{\infty} \frac{x^n}{n!}.

Solution:

sum(x^n/factorial(n), n, 0, oo)
e^x

We proceed by analyzing the use of power series in Sage. Below is an example of computing the square root of 2 using power series.

Example

Compute the square root of 2 using power series.

Solution:

K.<t> = PowerSeriesRing(QQ, 20)
# B is a power series of sqrt(1 + t) up to power 14 B = sqrt(1+t).truncate(15); B
-185725/33554432*t^14 + 52003/8388608*t^13 - 29393/4194304*t^12 + 4199/524288*t^11 - 2431/262144*t^10 + 715/65536*t^9 - 429/32768*t^8 + 33/2048*t^7 - 21/1024*t^6 + 7/256*t^5 - 5/128*t^4 + 1/16*t^3 - 1/8*t^2 + 1/2*t + 1
# our approximation for sqrt(2) sqrt2approx = B(1).n() sqrt2approx
1.41159650683403
# its square is close to 2 (sqrt2approx^2).n()
1.99260469810604

This example shows how to approximate π\pi ising power series.

# C is the series for acos(t) C = sin(t).reverse(); C
t + 1/6*t^3 + 3/40*t^5 + 5/112*t^7 + 35/1152*t^9 + 63/2816*t^11 + 231/13312*t^13 + 143/10240*t^15 + 6435/557056*t^17 + 12155/1245184*t^19 + O(t^20)
# we know that acos(0.5) = pi/6, so we should get pi 6*C.truncate(15)(0.5)
3.14158942531912

An illustration of a telescoping series

Series of the form n=1+an\sum_{n=1}^{+\infty}a_{n} with an=bnbn+1a_{n}=b_{n}-b_{n+1} are referred to as telescoping series. For example, a telescoping series is given below:

112+123+134++1k(k+1)+=(112)+(1213)+(1314)++(1k1k+1)+=112+1213+131k+1k1k+1+1k+1=1\frac{1}{1\cdot 2} + \frac{1}{2\cdot 3} + \frac{1}{3\cdot 4} + \dots + \frac{1}{k\cdot (k+1)} + \dots = (1 - \frac{1}{2}) + (\frac{1}{2} - \frac{1}{3}) + (\frac{1}{3} - \frac{1}{4}) + \dots + (\frac{1}{k} - \frac{1}{k+1}) + \dots = 1 - \frac{1}{2} + \frac{1}{2} - \frac{1}{3} + \frac{1}{3} - \dots - \frac{1}{k} + \frac{1}{k} - \frac{1}{k+1} + \frac{1}{k+1} - \dots = 1

To illustrate this series we can proceed as follows:

var('n'); var('k') a = 1/(k*(k+1))
B=Graphics() for n in srange(1, 100+1): B=B+points((n, sum(1/(k*(k+1)), k, 1, n))) B
Image in a Jupyter notebook

4.8. Taylor polynomials (Taylor expansions)

We now want to approximate functions, by Taylor polynomials.

Taylor polynomials are mathematical approximations of a function using a series of polynomials. Specifically, they are a way of representing a function as an infinite sum of terms, where each term is a polynomial multiplied by a power of the independent variable.

The Taylor polynomial of a function f(x)f(x) centered at x=ax=a is given by:

f(x)=f(a)+f(a)(xa)+f(a)/2!(xa)2+f(a)/3!(xa)3+...+f(n)(a)/n!(xa)n+Rn(x)f(x) = f(a) + f'(a)(x-a) + f''(a)/2!(x-a)^2 + f'''(a)/3!(x-a)^3 + ... + f^{(n)}(a)/n!(x-a)^n + R_n(x)

where

  • f(n)(a)f^{(n)}(a) denotes the nthnth derivative of f(x)f(x) evaluated at x=ax=a;

  • Rn(x)R_n(x) is the remainder term of the nthnth-degree Taylor polynomial.

The first few terms of the Taylor polynomial are given by:

  • First-degree Taylor polynomial (also known as the tangent line approximation): f(x)f(a)+f(a)(xa)f(x) ≈ f(a) + f'(a)(x - a)

  • Second-degree Taylor polynomial: f(x)f(a)+f(a)(xa)+f(a)/2!(xa)2f(x) ≈ f(a) + f'(a)(x - a) + f''(a)/2!(x - a)^2

  • Third-degree Taylor polynomial: f(x)f(a)+f(a)(xa)+f(a)/2!(xa)2+f(a)/3!(xa)3f(x) ≈ f(a) + f'(a)(x - a) + f''(a)/2!(x - a)^2 + f'''(a)/3!(x - a)^3

  • Fourth-degree Taylor polynomial: f(x)f(a)+f(a)(xa)+f(a)/2!(xa)2+f(a)/3!(xa)3+f(4)(a)/4!(xa)4f(x) ≈ f(a) + f'(a)(x - a) + f''(a)/2!(x - a)^2 + f'''(a)/3!(x - a)^3 + f^(4)(a)/4!(x - a)^4

and so on. As the degree of the polynomial increases, the approximation gets closer to the actual value of the function.

In SageMath we can find the Taylor polynomial of a given function f(x)f(x) in two ways:

  • apply the taylor(f(x),x,a,n){\tt{taylor(f(x), x, a, n)}} function;

  • apply the .taylor(x,a,n){\tt{.taylor(x,a,n) }} SAGE function method.

Let us present some examples.

Exercise

For the function sin(x)\sin(x) find the Taylor polynomial of degree at most 5, around 00.

Solution:

show("The Taylor polynomial of f, of degree at most 5, is the polynomial:", taylor(sin(x), x, 0, 5))

The Taylor polynomial of f, of degree at most 5, is the polynomial:1120x516x3+x\displaystyle \verb|The|\verb| |\verb|Taylor|\verb| |\verb|polynomial|\verb| |\verb|of|\verb| |\verb|f,|\verb| |\verb|of|\verb| |\verb|degree|\verb| |\verb|at|\verb| |\verb|most|\verb| |\verb|5,|\verb| |\verb|is|\verb| |\verb|the|\verb| |\verb|polynomial:| \frac{1}{120} \, x^{5} - \frac{1}{6} \, x^{3} + x

Or we could directly type:

show(sin(x).taylor(x, 0, 5))

1120x516x3+x\displaystyle \frac{1}{120} \, x^{5} - \frac{1}{6} \, x^{3} + x

Remark

To get a better view of how the approximation works, one can plot both the given functions as follows (below we will described another procedure for doing this)

( plot(sin(x), xmax=5, xmin=-5, legend_label=r'$sin(x)$') + plot(taylor(sin(x), x, 0, 5), xmax=5, xmin=-5, color='red', legend_label=r'$\frac{1}{120}x^5-\frac{1}{6}x^3+x$') )
Image in a Jupyter notebook

As you can see our Tailor polynomial of degree 5 approximates the function pretty well to some extent

Exercise

For the exponential function f(x)=exf(x)=e^x find the Taylor polynomial of degree at most 10, around 0.

Solution:

show("The Taylor polynomial in question has the form:", taylor(e^x, x, 0, 10))

The Taylor polynomial in question has the form:13628800x10+1362880x9+140320x8+15040x7+1720x6+1120x5+124x4+16x3+12x2+x+1\displaystyle \verb|The|\verb| |\verb|Taylor|\verb| |\verb|polynomial|\verb| |\verb|in|\verb| |\verb|question|\verb| |\verb|has|\verb| |\verb|the|\verb| |\verb|form:| \frac{1}{3628800} \, x^{10} + \frac{1}{362880} \, x^{9} + \frac{1}{40320} \, x^{8} + \frac{1}{5040} \, x^{7} + \frac{1}{720} \, x^{6} + \frac{1}{120} \, x^{5} + \frac{1}{24} \, x^{4} + \frac{1}{6} \, x^{3} + \frac{1}{2} \, x^{2} + x + 1

Exercise

For the function f(x)=x2+log(x)f(x)=x^2+\log(x) find the Taylor polynomial of degree 2 and 8, for x=3. Moreover, plot the approximations.

Solution:

f(x) = x^2+log(x) taylor_2_deg = taylor(f(x), x, 3, 2) taylor_8_deg = taylor(f(x), x, 3, 8) show(taylor_2_deg) show(taylor_8_deg)

1718(x3)2+193x+log(3)10\displaystyle \frac{17}{18} \, {\left(x - 3\right)}^{2} + \frac{19}{3} \, x + \log\left(3\right) - 10

152488(x3)8+115309(x3)714374(x3)6+11215(x3)51324(x3)4+181(x3)3+1718(x3)2+193x+log(3)10\displaystyle -\frac{1}{52488} \, {\left(x - 3\right)}^{8} + \frac{1}{15309} \, {\left(x - 3\right)}^{7} - \frac{1}{4374} \, {\left(x - 3\right)}^{6} + \frac{1}{1215} \, {\left(x - 3\right)}^{5} - \frac{1}{324} \, {\left(x - 3\right)}^{4} + \frac{1}{81} \, {\left(x - 3\right)}^{3} + \frac{17}{18} \, {\left(x - 3\right)}^{2} + \frac{19}{3} \, x + \log\left(3\right) - 10

We can check the accuracy of approximation near the point aa

print('Original function value at x=3:', round(f(x=3.1),10), '\n2nd degree Taylor aproximation:', round(taylor_2_deg(x=3.1), 10), '\n8th degree Taylor aproximation:', round(taylor_8_deg(x=3.1), 10) )
Original function value at x=3: 10.7414021115 2nd degree Taylor aproximation: 10.7413900664 8th degree Taylor aproximation: 10.7414021115

As we can see, the higher the degree of approximation, the more precise are the values near the center aa. Let us now plot the approximations:

( plot(f, xmax=10, xmin=-3, legend_label=r'$x^3log(x)$') + plot(taylor_2_deg, xmax=10, xmin=-3, color='green', alpha=0.7, legend_label='2-degree Taylor poly') + plot(taylor_8_deg, xmax=10, xmin=-3, color='red', legend_label='8-degree Taylor poly') )
verbose 0 (3899: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 47 points. verbose 0 (3899: plot.py, generate_plot_points) Last error message: 'Unable to compute f(-0.01964549734289004)'
Image in a Jupyter notebook

Exercise (optional) - Construction of a Taylor series function in Python

The following three exercises are optional, as we dont use built-in functions from Sage to study Taylor series.

In these examples we are totally based on Python so one may skip them.

Construct your own Taylor series Sage function for f(x)=sin(x)f(x)=\sin(x). Note that this function should satisfy the following:

  • It takes as an input the point x_0 around which we want to consider the Taylor expansion.

  • It takes as an input the function f for which the Taylor series is computed.

  • It takes as an input the order at which we want to stop the Taylor expansion.

  • The resulting Taylor series is a function of one free variable, i.e., it can be evaluated at some point x_value

Solution:

A possible solution goes as follows:

# Import necessary module import matplotlib.pyplot as plt import numpy as np # Define the variable and function x = var('x') f = function('f')(x) f = sin(x) # Define a function to compute the Taylor polynomial def taylor_series(x_value, x0, f, order): approximation = 0 for n in range(order + 1): term = diff(f, x, n).subs(x=x0) * (x_value - x0)**n / factorial(n) approximation += term return approximation

We can pick some concrete inputs to see that our Sage function "taylor_series{\tt{taylor\_series}}" for Taylor polynomial works as it should. For example, Taylor polynomial of degree 5 of sin(x) at 0 is

show(taylor_series(x, 0, f, 5))

1120x516x3+x\displaystyle \frac{1}{120} \, x^{5} - \frac{1}{6} \, x^{3} + x

Notice that Taylor polynomial of degree 6 of sin(x) at 0 is the same as the degree 5 Taylor polynomial

show(taylor_series(x, 0, f, 6))

1120x516x3+x\displaystyle \frac{1}{120} \, x^{5} - \frac{1}{6} \, x^{3} + x

This is becaue sin(x) is an odd function, and thus the Taylor series contains only summands with odd exponents. Indeed, observe the change if we shift to degree 7

show(taylor_series(x, 0, f, 7))

15040x7+1120x516x3+x\displaystyle -\frac{1}{5040} \, x^{7} + \frac{1}{120} \, x^{5} - \frac{1}{6} \, x^{3} + x

Exercise (optional)

Compare the graphs of the Taylor polynomials of the sin function, centered at the same point x0=0x_0=0.

Solution:

(Note that the Sage cell given below will work only if the previous cell, where the function "taylor_series{\tt{taylor\_series}}" is defined, has already been compiled).

# Define the x range x_values = np.linspace(-8, 8, 1000) # Plot the function and the Taylor polynomial approximations fig, ax = plt.subplots(figsize=(10, 6)) # Hide the right and top spines ax.spines['right'].set_visible(False) ax.spines['top'].set_visible(False) # Only show ticks on the left and bottom spines ax.yaxis.set_ticks_position('left') ax.xaxis.set_ticks_position('bottom') # Move bottom spine (x-axis) to y=0 ax.spines['bottom'].set_position(('data',0)) ax.plot(x_values, np.sin(x_values), label='sin(x)') # Calculate the Taylor polynomial approximations and plots for order in range(1,15, 2): taylor_approximation = [taylor_series(val, 0, f, order) for val in x_values] ax.plot(x_values, taylor_approximation, label=f'order {order}') # Set the x axis label ax.set_xlabel('x') # Set the y-axis limits ax.set_ylim([-8, 8]) # Show the legend ax.legend() # Show the plot plt.show()
Image in a Jupyter notebook

Exercise (optional)

Repeat the above procedure and display a complete code for the computation and visualization of Taylor polynomials of the exponential function around the origin. In particular, chose to display the graphs of Taylor polynomials of orders 4,5,,94,5, \ldots, 9. (Note that exe^x is neither odd, nor an even function, so there is no reason to expect some vanishing of terms in the Taylor expansion, as we saw above for the case of sin(x)\sin(x)).

Solution:

# Calculate the Taylor polynomials of degrees k = 1,2,...5 for k in range(1,6): print('order', k, 'Taylor polynomial:') show(taylor_series(x, 0, exp(x), k))
order 1 Taylor polynomial:

x+1\displaystyle x + 1

order 2 Taylor polynomial:

12x2+x+1\displaystyle \frac{1}{2} \, x^{2} + x + 1

order 3 Taylor polynomial:

16x3+12x2+x+1\displaystyle \frac{1}{6} \, x^{3} + \frac{1}{2} \, x^{2} + x + 1

order 4 Taylor polynomial:

124x4+16x3+12x2+x+1\displaystyle \frac{1}{24} \, x^{4} + \frac{1}{6} \, x^{3} + \frac{1}{2} \, x^{2} + x + 1

order 5 Taylor polynomial:

1120x5+124x4+16x3+12x2+x+1\displaystyle \frac{1}{120} \, x^{5} + \frac{1}{24} \, x^{4} + \frac{1}{6} \, x^{3} + \frac{1}{2} \, x^{2} + x + 1

Let us now display the code for the computation and visualization of Taylor polynomials of order 4,5,,94, 5, \ldots, 9 of exe^x around 00.

# Import necessary module import matplotlib.pyplot as plt import numpy as np # Define the variable and function x = var('x') f = function('f')(x) f = np.exp(x) # Define a function to compute the Taylor polynomial def taylor_series(x_value, x0, f, order): approximation = 0 for n in range(order + 1): term = diff(f, x, n).subs(x=x0) * (x_value - x0)**n / factorial(n) approximation += term return approximation # Define the x range x_values = np.linspace(-8, 8, 1000) # Plot the function and the Taylor polynomial approximations fig, ax = plt.subplots(figsize=(10, 6)) # Hide the right and top spines ax.spines['right'].set_visible(False) ax.spines['top'].set_visible(False) # Only show ticks on the left and bottom spines ax.yaxis.set_ticks_position('left') ax.xaxis.set_ticks_position('bottom') # Move bottom spine (x-axis) to y=0 ax.spines['bottom'].set_position(('data',0)) ax.plot(x_values, np.exp(x_values), label='exp(x)', linewidth=2.5) # Calculate the Taylor polynomial approximations and plot them for order in range(4,10): taylor_approximation = [taylor_series(val, 0, f, order) for val in x_values] ax.plot(x_values, taylor_approximation, label=f'order {order}') # Set the x axis label ax.set_xlabel('x') # Set the y-axis limits ax.set_ylim([-6, 6]) # Show the legend ax.legend() # Show the plot plt.show()
Image in a Jupyter notebook

4.9. Integration

We now proceed with learning how to compute simple integrals via Sage. First we recall some theoretical background.

Theoretical remarks on integration

An integral is a mathematical concept that expresses the area under a curve. It is a way to calculate the length, area, or volume of an irregular object or shape. There are two types of integrals, definite and indefinite.

A definite integral is used to find the area under a curve between two specific points. It is written as:

abf(x)dx\int_a^b f(x) \, dx

where f(x)f(x) is the function being integrated, aa is the lower limit and bb is the upper limit of integration.

An indefinite integral, also called an antiderivative, is the set of all functions that would have f(x)f(x) as their derivative. It is written as:

f(x)dx\int f(x) \, dx

Here are some common integrals:

  1. (a)dx=ax+C\int (a) dx = ax + C, where aa is a coefficient, and CC is a constant of integration.

  2. xndx=xn+1n+1+C\int x^n dx = \frac{x^{n+1}}{n+1} + C, where n1n \neq -1.

  3. dxx=logx+C\int \frac{dx}{x} = \log{|x|} + C.

  4. axdx=axlnx+C\int a^x dx = \frac{a^x}{\ln{x}} + C.

  5. exdx=ex+C\int e^x dx = e^x + C.

  6. cos(x)dx=sin(x)+C\int \cos(x) dx = \sin(x) + C.

  7. sin(x)dx=cos(x)+C\int \sin(x) dx = -\cos(x) + C.

Now, in order to compute/find the definite integral abf(x)dx \int_a^b f(x) dx we can use:

  1. integral(f(x), x, a, b) function;

  2. f.integral(x,a,b) or f.integrate(x,a,b) methods.

To find the indefinite integral we can use the same:

  1. integral(f(x), x) function;

  2. f.integral(x) or f.integrate(x) methods.

but we need to omit the aa and bb parameters, lets have a look

Remark

Note that Sage does not add the constant of integration (the CC constant posed above)

Example

Conside the function f(x)=x37x+1f(x)=x^3-7x+1.

  1. Use Sage to plot the graph of ff for 2x2-2\leq x\leq 2 and the area bounded by the graph of ff an the xx-axis (andf the vertical lines x=±2x=\pm 2).

  2. Compute the integral 22f(x)dx\displaystyle\int_{-2}^{2}f(x)dx and next provide a geometric interpretation.

Solution:

f(x)=x^3-7*x+1 a=-2 b=2 plot(f(x), (x, a, b), color="darkgreen", figsize=4)
Image in a Jupyter notebook

To mark the region bounded by the graph of ff, the xx-axis and the lines x=±1x=\pm 1, we can instead type:

plot(f(x), (x, a, b), color="darkgreen", figsize=4, fill= True, fillcolor="grey")
Image in a Jupyter notebook
integral(f(x), x, a, b)
4
General geometric interpretation:

The integral abf(x)dx\displaystyle\int_{a}^{b}f(x)dx, when it exists, represents the net area between the curve of the function f(x)f(x), and the xx-axis, over the interval [a,b][a, b]

For our case the integral represents the signed area between the curve of f(x)=x35x+1f(x)=x^3-5x+1, the xx-axis and the vertical lines x=±2x=\pm 2. Here, by the term signed area we mean that if the curve is above the xx-axis, then the area contributes positively to the integral. If the curve is below the xx-axis, then the area contributes negatively.

The total area is "net" in the sense that the areas above the xx-axis are positive, and those below the xx-axis are negative. Thus, if part of the curve lies above the xx-axis and part lies below, the positive and negative areas are subtracted from each other.

Example (Riemann sums)

For the theory on Riemann sums we refer to Chapter 6 in the BG book. Below we will use the function f(x)=x37x+1f(x)=x^3-7x+1 to illustrate them and also explain how we can use them to approximate an integral.

f(x)=x^3-7*x+1 a=-2 b=2 n=10 ### r=0 for a left Riemann sum, r=1 for a right Riemman sum, r=1/2 for a middle on r=0 delta=(b-a)/n; rdelta=r*delta; xk=a; L=[]; S=0 for k in range(n): L=L + [(xk, 0)] y=f(xk+rdelta) S=S + y L=L + [(xk, y)] xk= xk + delta L=L + [(xk, y)] L=L + [(xk, 0)] G=plot(f(x), (x, a, b), color="darkblue", thickness=2) G=G+plot(f(x), (x, a, b), color="green", thickness=1, fill=True, fillcolor="grey") G=G+polygon(L, edgecolor="red", color="lightblue") G.show(aspect_ratio="automatic", figsize=6)
Image in a Jupyter notebook
f(x)=x^3-7*x+1 a=-2 b=2 n=10 ### r=0 for a left Riemann sum, r=1 for a right Riemman sum, r=1/2 for a middle on r=1 delta=(b-a)/n; rdelta=r*delta; xk=a; L=[]; S=0 for k in range(n): L=L + [(xk, 0)] y=f(xk+rdelta) S=S + y L=L + [(xk, y)] xk= xk + delta L=L + [(xk, y)] L=L + [(xk, 0)] G=plot(f(x), (x, a, b), color="darkblue", thickness=2) G=G+plot(f(x), (x, a, b), color="green", thickness=1, fill=True, fillcolor="grey") G=G+polygon(L, edgecolor="red", color="lightblue") G.show(aspect_ratio="automatic", figsize=6)
Image in a Jupyter notebook
f(x)=x^3-7*x+1 a=-2 b=2 n=10 ### r=0 for a left Riemann sum, r=1 for a right Riemman sum, r=1/2 for a middle on r=1/2 delta=(b-a)/n; rdelta=r*delta; xk=a; L=[]; S=0 for k in range(n): L=L + [(xk, 0)] y=f(xk+rdelta) S=S + y L=L + [(xk, y)] xk= xk + delta L=L + [(xk, y)] L=L + [(xk, 0)] G=plot(f(x), (x, a, b), color="darkblue", thickness=2) G=G+plot(f(x), (x, a, b), color="green", thickness=1, fill=True, fillcolor="grey") G=G+polygon(L, edgecolor="red", color="lightblue") G.show(aspect_ratio="automatic", figsize=6)
Image in a Jupyter notebook

Observe for all the three case (left, right and middle Rieman sums) we can get better approximations by increasing nn.

We can use this method to estimate the integral (or area) as follows:

f(x)=x^3-7*x+1 a=-2 b=2 n=20 ### r=0 for a left Riemann sum, r=1 for a right Riemman sum, r=1/2 for a middle on r=0 I=integral(f(x), x, a, b).n() delta=(b-a)/n; rdelta=r*delta; xk=a; L=[]; S=0 for k in range(n): L= L + [(xk, 0)] y= f(xk+rdelta) S= S + y L= L + [(xk, y)] xk= xk + delta L= L + [(xk, y)] S= delta*S.n() L= L + [(xk, 0)] show("The actuall integral is given by:", I) show("The approximate left Riemann Sum is given by:", S) G=plot(f(x), (x, a, b), color="darkblue", thickness=2) G=G+plot(f(x), (x, a, b), color="green", thickness=1, fill=True, fillcolor="grey") G=G+polygon(L, edgecolor="red", color="lightblue") G.show(aspect_ratio="automatic", figsize=6)

The actuall integral is given by:4.00000000000000\displaystyle \verb|The|\verb| |\verb|actuall|\verb| |\verb|integral|\verb| |\verb|is|\verb| |\verb|given|\verb| |\verb|by:| 4.00000000000000

The approximate left Riemann Sum is given by:5.20000000000000\displaystyle \verb|The|\verb| |\verb|approximate|\verb| |\verb|left|\verb| |\verb|Riemann|\verb| |\verb|Sum|\verb| |\verb|is|\verb| |\verb|given|\verb| |\verb|by:| 5.20000000000000

Image in a Jupyter notebook

By increasing nn we can get a much better approximation: (similarly are treated the right or middle Riemann sums).

f(x)=x^3-7*x+1 a=-2 b=2 n=200 ### r=0 for a left Riemann sum, r=1 for a right Riemman sum, r=1/2 for a middle on r=0 I=integral(f(x), x, a, b).n() delta=(b-a)/n; rdelta=r*delta; xk=a; L=[]; S=0 for k in range(n): L= L + [(xk, 0)] y= f(xk+rdelta) S= S + y L= L + [(xk, y)] xk= xk + delta L= L + [(xk, y)] S= delta*S.n() L= L + [(xk, 0)] show("The actuall integral is given by:", I) show("The approximate left Riemann Sum is given by:", S) G=plot(f(x), (x, a, b), color="darkblue", thickness=2) G=G+plot(f(x), (x, a, b), color="green", thickness=1, fill=True, fillcolor="grey") G=G+polygon(L, edgecolor="red", color="lightblue") G.show(aspect_ratio="automatic", figsize=6)

The actuall integral is given by:4.00000000000000\displaystyle \verb|The|\verb| |\verb|actuall|\verb| |\verb|integral|\verb| |\verb|is|\verb| |\verb|given|\verb| |\verb|by:| 4.00000000000000

The approximate left Riemann Sum is given by:4.12000000000000\displaystyle \verb|The|\verb| |\verb|approximate|\verb| |\verb|left|\verb| |\verb|Riemann|\verb| |\verb|Sum|\verb| |\verb|is|\verb| |\verb|given|\verb| |\verb|by:| 4.12000000000000

Image in a Jupyter notebook

Example (computation of different types of integrals via SageMath)

f = x^2 # define the function # Indefinite integral print('Get indefinite integral:') show('First way via "integral()" function: ', integral(f, x)) show('Second way via "f.integral()" method: ', f.integral(x)) show('Third way via "f.integrate()" method: ', f.integrate(x)) print('\n') # Definite integral print('Get definite integral:') a = 0 # define a (lower bound) b = 5 # define b (upper bound) show('First way via "integral()" function: ', integral(f, x, a, b)) show('Second way via "f.integral()" method: ', f.integral(x, a, b)) show('Third way via "f.integrate()" method: ', f.integrate(x, a, b))
Get indefinite integral:

First way via "integral()" function:13x3\displaystyle \verb|First|\verb| |\verb|way|\verb| |\verb|via|\verb| |\verb|"integral()"|\verb| |\verb|function:| \frac{1}{3} \, x^{3}

Second way via "f.integral()" method:13x3\displaystyle \verb|Second|\verb| |\verb|way|\verb| |\verb|via|\verb| |\verb|"f.integral()"|\verb| |\verb|method:| \frac{1}{3} \, x^{3}

Third way via "f.integrate()" method:13x3\displaystyle \verb|Third|\verb| |\verb|way|\verb| |\verb|via|\verb| |\verb|"f.integrate()"|\verb| |\verb|method:| \frac{1}{3} \, x^{3}

Get definite integral:

First way via "integral()" function:1253\displaystyle \verb|First|\verb| |\verb|way|\verb| |\verb|via|\verb| |\verb|"integral()"|\verb| |\verb|function:| \frac{125}{3}

Second way via "f.integral()" method:1253\displaystyle \verb|Second|\verb| |\verb|way|\verb| |\verb|via|\verb| |\verb|"f.integral()"|\verb| |\verb|method:| \frac{125}{3}

Third way via "f.integrate()" method:1253\displaystyle \verb|Third|\verb| |\verb|way|\verb| |\verb|via|\verb| |\verb|"f.integrate()"|\verb| |\verb|method:| \frac{125}{3}

Another example:

integral((x-1)/(x^2+9),x)
-1/3*arctan(1/3*x) + 1/2*log(x^2 + 9)

If we want to have a Latex output style, we can add:

show(_)

13arctan(13x)+12log(x2+9)\displaystyle -\frac{1}{3} \, \arctan\left(\frac{1}{3} \, x\right) + \frac{1}{2} \, \log\left(x^{2} + 9\right)

Exercise

Find the antiderivative of the function f(x)=x3+2x2+3x+4f(x) = x^3 + 2x^2 + 3x + 4

Solution

integral(x^3 + 2*x^2 + 3*x + 4, x)
1/4*x^4 + 2/3*x^3 + 3/2*x^2 + 4*x

Exercise

Evaluate the definite integral of the function f(x)=1x2+4x+5f(x) = \frac{1}{x^2 + 4x + 5} from x=0x=0 to x=1x=1

Solution:

f = (x^2 + 4*x + 5) show('Indef integral = ', f.integral(x)) show('Answer:', f.integrate(x, 0, 1 ))

Indef integral =13x3+2x2+5x\displaystyle \verb|Indef|\verb| |\verb|integral|\verb| |\verb|=| \frac{1}{3} \, x^{3} + 2 \, x^{2} + 5 \, x

Answer:223\displaystyle \verb|Answer:| \frac{22}{3}

Exercise

Find the integral of the function f(x)=2x3+3x24x+1x2+1.f(x) = \frac{2x^3 + 3x^2 - 4x + 1}{x^2 + 1}.

Solution:

integral((2*x^3 + 3*x^2 - 4*x + 1) / (x^2 + 1), x)
x^2 + 3*x - 2*arctan(x) - 3*log(x^2 + 1)

Exercise

Calculate the definite integral of the function f(x)=ln(x)xf(x) = \frac{\operatorname{ln}(x)}{x} from x=1x=1 to x=ex=e

Solution:

f(x) = ln(x) / x f.integral(x, 1, e)
1/2

Exercise

Check the integration by parts formula abf(x)g(x)dx=f(b)g(b)f(a)g(a)abf(x)g(x)dx\int_a^b f(x) g'(x) dx = f(b)g(b) - f(a)g(a) - \int_a^b f'(x) g(x) dx for f(x)=sin(x),g(x)=x2f(x) = sin(x), g(x) = x^2 and a=0,b=π2a = 0, b = \frac{\pi}{2}.

Solution:

f = sin(x) g = x^2 a = 0 b = pi/2 LHS = integral(f*diff(g,x), x, a , b) RHS = f(x=b)*g(x=b) - f(x=a)*g(x=a) - integral(diff(f,x)*g, x, a , b) print(f'Left hand side = {LHS}') print(f'Right hand side = {RHS}')
Left hand side = 2 Right hand side = 2

It is important to understand that Sage can also handle definite integrals involving variables, as in the case of derivatives. Let us treat such an example.

Exercise

Compute the integral I=01(ax3+bx2+cx+d)dx. I=\displaystyle\int_{0}^{1}(ax^3+bx^2+cx+d)\,dx\,.

Solution:

This is an application of symbolic variables as first we should introduce the variables a,b,c,da, b, c, d.

var("a, b, c, d") I=integral(a*x^3+b*x^2+c*x+d, x, 0, 1) show("The integral I is given in tersm of the parameteres a, b, c, d, as follows I=", I)

The integral I is given in tersm of the parameteres a, b, c, d, as follows I=14a+13b+12c+d\displaystyle \verb|The|\verb| |\verb|integral|\verb| |\verb|I|\verb| |\verb|is|\verb| |\verb|given|\verb| |\verb|in|\verb| |\verb|tersm|\verb| |\verb|of|\verb| |\verb|the|\verb| |\verb|parameteres|\verb| |\verb|a,|\verb| |\verb|b,|\verb| |\verb|c,|\verb| |\verb|d,|\verb| |\verb|as|\verb| |\verb|follows|\verb| |\verb|I=| \frac{1}{4} \, a + \frac{1}{3} \, b + \frac{1}{2} \, c + d

**Fundamental Theorem of Integral Calculus - applications **

Let us now discuss a few tasks related to the fundamental theorem of calculus, see Chapter 6 on the BG book for details.

Here we only recall that given some continuous function f:[a,b]Rf : [a, b]\to\mathbb{R}, defined on a finite interval [a,b][a, b], we can introduce the function F(x)=axf(t)dt,x[a,b]. F(x)=\int_{a}^{x}f(t)dt\,, \quad x\in[a, b]\,.

Its derivative at xx is given by F(x)=f(x)F'(x)=f(x).

Exercise

Introduce in SageMath the fucntion

F(x)=0xet2(t2+1)dt.F(x)=\int_{0}^{x}\frac{e^{t^2}}{(t^2+1)}\,dt\,.

Next derive the derivative of FF and compute F(1)F'(1).

Solution:

var("x, t") f(t)=e^(t^2)/(t^2 + 1) F(x) = integral(f(t), t, 0, x) show("The function F is given by F(x)=", F(x)) Fprime(x)=diff(F(x), x) show("The derivative of F is given by F'(x)=", Fprime(x)) show("The value of F'(1) is given by:", Fprime(1))

The function F is given by F(x)=0xe(t2)t2+1dt\displaystyle \verb|The|\verb| |\verb|function|\verb| |\verb|F|\verb| |\verb|is|\verb| |\verb|given|\verb| |\verb|by|\verb| |\verb|F(x)=| \int_{0}^{x} \frac{e^{\left(t^{2}\right)}}{t^{2} + 1}\,{d t}

The derivative of F is given by F'(x)=e(x2)x2+1\displaystyle \verb|The|\verb| |\verb|derivative|\verb| |\verb|of|\verb| |\verb|F|\verb| |\verb|is|\verb| |\verb|given|\verb| |\verb|by|\verb| |\verb|F'(x)=| \frac{e^{\left(x^{2}\right)}}{x^{2} + 1}

The value of F'(1) is given by:12e\displaystyle \verb|The|\verb| |\verb|value|\verb| |\verb|of|\verb| |\verb|F'(1)|\verb| |\verb|is|\verb| |\verb|given|\verb| |\verb|by:| \frac{1}{2} \, e

Be aware that in SageMathCell in order to obtain the result posed above we shoud include in our program the command assume(x>0){\tt{assume(x>0)}}.

In fact below we describe a bit quicker way to introduce FF.

# Define the variables x and t var('x, t') assume(x>0) # Define the function f(t) f = exp(t^2) / (t^2 + 1) # Define the integral function F(x) F = integral(f, t, 0, x) # Display the function F(x) show(F)

0xe(t2)t2+1dt\displaystyle \int_{0}^{x} \frac{e^{\left(t^{2}\right)}}{t^{2} + 1}\,{d t}

Or finally one could type

var("x, t") f(t)=e^(t^2)/(t^2 + 1) F=f.integral(t, 0, x); show(F)

0xe(t2)t2+1dt\displaystyle \int_{0}^{x} \frac{e^{\left(t^{2}\right)}}{t^{2} + 1}\,{d t}

Fprime=F.derivative(x); show(Fprime)

e(x2)x2+1\displaystyle \frac{e^{\left(x^{2}\right)}}{x^{2} + 1}

However, in this case we cannot type Fprime(1){\tt{Fprime(1)}} as we will get an error (try it!).

Using this approach, to evaluate F(1)F'(1) we should type the following:

show(F.derivative(x)(x=1))

12e\displaystyle \frac{1}{2} \, e

Exercise for practice

  1. Use Sage to evaluate the derivatives G(1)G'(1) and G(1)G''(1), where

G(x)=1xte3tdt.G(x)=\int_{1}^{x} te^{3t}\,dt\,.

Next confirm Sage's result by hand.

  1. Let f(t)=tt2f(t)=t-t^2 with t0t\geq 0, and F(x)=0xf(t)dtF(x)=\displaystyle\int_{0}^{x} f(t)\,dt. Find the positive real xx where F(x)F(x) starts decreasing.

Remark

Some functions do not have elementary antiderivatives. For example, consider the integral I=01ex2dx. I=\int_{0}^{1}e^{-x^2}\,dx\,.

show(integral(e^(-x^2),x,0,1))

12πerf(1)\displaystyle \frac{1}{2} \, \sqrt{\pi} \operatorname{erf}\left(1\right)

Let us ask Sage for the error function erf{\tt{erf}}

erf?
Type: LazyImport String form: erf File: /ext/sage/9.7/src/sage/misc/lazy_import.pyx Docstring: The error function. The error function is defined for real values as \operatorname{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2} dt. This function is also defined for complex values, via analytic continuation. EXAMPLES: We can evaluate numerically: sage: erf(2) erf(2) sage: erf(2).n() 0.995322265018953 sage: erf(2).n(100) 0.99532226501895273416206925637 sage: erf(ComplexField(100)(2+3j)) -20.829461427614568389103088452 + 8.6873182714701631444280787545*I Basic symbolic properties are handled by Sage and Maxima: sage: x = var("x") sage: diff(erf(x),x) 2*e^(-x^2)/sqrt(pi) sage: integrate(erf(x),x) x*erf(x) + e^(-x^2)/sqrt(pi) ALGORITHM: Sage implements numerical evaluation of the error function via the "erf()" function from mpmath. Symbolics are handled by Sage and Maxima. REFERENCES: * https://en.wikipedia.org/wiki/Error_function * http://mpmath.googlecode.com/svn/trunk/doc/build/functions/expin tegrals.html#error-functions Class docstring: EXAMPLES: sage: from sage.misc.lazy_import import LazyImport sage: my_integer = LazyImport('sage.rings.all', 'Integer') sage: my_integer(4) 4 sage: my_integer('101', base=2) 5 sage: my_integer(3/2) Traceback (most recent call last): ... TypeError: no conversion of this rational to integer Init docstring: EXAMPLES: sage: from sage.misc.lazy_import import LazyImport sage: lazy_ZZ = LazyImport('sage.rings.all', 'ZZ') sage: type(lazy_ZZ) <class 'sage.misc.lazy_import.LazyImport'> sage: lazy_ZZ._get_object() is ZZ True sage: type(lazy_ZZ) <class 'sage.misc.lazy_import.LazyImport'> Call docstring: Calling self calls the wrapped object. EXAMPLES: sage: from sage.misc.lazy_import import LazyImport sage: my_isprime = LazyImport('sage.all', 'is_prime') sage: is_prime(12) == my_isprime(12) True sage: is_prime(13) == my_isprime(13) True

Hence, according to Sage the error function is defined as erf(x)=2π0xet2dt. \operatorname{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2} dt.

Let us evaluate this fucntion at some points:

erf(2)
erf(2)
erf(2).n()
0.995322265018953

Finally, we could directly aapproximate the initial intregral based on the N(){\tt{N()}} function.

In the final chapter we will see that Sage provides another built-in method to approximate integrals numerically.

N(integral(e^(-x^2),x,0,1))
0.746824132812427

Another example where we get an expression in terms of the error function is for example the following integral:

sin(x2)dx.\int\sin(x^2)\, dx\,.
show(integral(sin(x^2), x))

116π((i+1)2erf((12i+12)2x)+(i1)2erf((12i12)2x)(i1)2erf(ix)+(i+1)2erf((1)14x))\displaystyle \frac{1}{16} \, \sqrt{\pi} {\left(\left(i + 1\right) \, \sqrt{2} \operatorname{erf}\left(\left(\frac{1}{2} i + \frac{1}{2}\right) \, \sqrt{2} x\right) + \left(i - 1\right) \, \sqrt{2} \operatorname{erf}\left(\left(\frac{1}{2} i - \frac{1}{2}\right) \, \sqrt{2} x\right) - \left(i - 1\right) \, \sqrt{2} \operatorname{erf}\left(\sqrt{-i} x\right) + \left(i + 1\right) \, \sqrt{2} \operatorname{erf}\left(\left(-1\right)^{\frac{1}{4}} x\right)\right)}

I=integral(sin(x^2), x) N(I) #this will produce an error, try it! I.n() #similarly we will get an error

An example of numerical integration in SageMath

Compute the integral

I=12x+1+x2x.I=\int_{1}^{2} \frac{ \sqrt{x+\sqrt{1+x^2}}}{x} \,.
f(x)=sqrt(x+sqrt(1+x^2))/x fplot=f.plot(x, 1, 100) show(fplot) fint=integral(f(x), x) show(fint) # we will get a very complicated expression depending on Gamma function
Image in a Jupyter notebook

xΓ(14)Γ(14)23F2(14,14,1412,34;1x2)8πΓ(34)\displaystyle \frac{\sqrt{x} \Gamma\left(\frac{1}{4}\right) \Gamma\left(-\frac{1}{4}\right)^{2} \,_3F_2\left(\begin{matrix} -\frac{1}{4},-\frac{1}{4},\frac{1}{4} \\ \frac{1}{2},\frac{3}{4} \end{matrix} ; -\frac{1}{x^{2}} \right)}{8 \, \pi \Gamma\left(\frac{3}{4}\right)}

Let us introduce and have a look at the so called Gamma function.

var("x, t") assume(x>0) Gamma(x)=integral(t^(x-1)*e^(-t), t, 0, oo) show(Gamma(x)) Gamma(0)

Γ(x)\displaystyle \Gamma\left(x\right)

Infinity
Gamma(0)==gamma(0)
Infinity == Infinity

Hence our Gamma function at least at x=0x=0 gives the same result with the built-in function in Sage that represents the Γ\Gamma function.

Let us try to plot our function for certain xx:

plot(Gamma, x, -10, 10, ymin=-10, ymax=10, detect_poles="True")
Image in a Jupyter notebook

We can plot in a better way the Gamma function by importing in it!

# Import the Gamma function from sage.functions.gamma import gamma # Define the range for the plot x_range = (-5, 5) # Avoid zero because Gamma(x) has a singularity at x=0 # Plot the Gamma function plot(gamma(x), x_range, title="Gamma Function", color='blue', ymin=-10, ymax=10, detect_poles="True")
Image in a Jupyter notebook

Type gamma?{\tt{gamma?}} to learn more information for the Gamma function in SageMath and see also Chapter 6 on the BG book.

gamma(0)
Infinity

Let us now return to our initial problem. One way to obtain a meaninfull result is to apply the following method:

(integral(f(x), x, 1, 2)).n()
1.23692911795883

Let us now use Sage via the built-in function numerical_integral{\tt{numerical\_integral}}. This approximates the value of an integral on an interval.

numerical_integral(f(x), 1, 2) # this will provide an approximation of the integral between 0 and 1
(1.2369291179588335, 1.3732671865871831e-14)

This provides a good approximation, as the error (the second component in Sage's output above is very small). However, observe that:

(integral(f(x), x, 0, 1)).n()
2.72269476898150

But:

numerical_integral(f(x), 0, 1)
(69.99281746999927, 13.112269516493129)

So in this case the numericalintegral{\tt{numerical_integral}} method provides a very big error and the estimation is wrong, and not acceptable. We will discuss mathematical techniques of numerical integration and how we can implement them using Sage in our final chapter, the final week of the course.

Improper integrals

Compute the integral

I=0sin(x)xdx.I=\int_{0}^{\infty}\frac{\sin(x)}{x}\, dx\,.
I=integral(sin(x)/x, x, 0, infinity) show("The integral in question equals", ":", I)

The integral in question equals:12π\displaystyle \verb|The|\verb| |\verb|integral|\verb| |\verb|in|\verb| |\verb|question|\verb| |\verb|equals| \verb|:| \frac{1}{2} \, \pi

numerical_integral(sin(x)/x, 0, oo)
(1.902348806239928, 1.6980673032591547)

Let us compare this result with the previous result:

n(pi/2)
1.57079632679490

Thus, the numerical_integral{\tt{numerical\_integral}} method is not very accurate in this case, and has a very high error, close to 1.7.