| Hosted by CoCalc | Download
Kernel: SageMath 9.1

Lecture 4 : Calculus and power series

References:

Summary:
In this lecture, we start by recalling how new functions can be defined in Python. Then we show how to plot functions in SageMath and discuss concepts from calculus (or analysis) like limits, derivatives and integrals. We finish by discussing power series, which are a generalization of the Taylor series from analysis.

Python warmup: functions

We have already seen many functions implemented in SageMath (such as log, span or solve). But it is also possible (and often very useful) to write functions ourselves.

def parabola(x): x_squared = x^2 return x_squared

This creates a new function parabola, which we can immediately call on some value:

parabola(4)

The general code pattern for defining a function foo with nn inputs (or arguments) looks as follows:

def foo(arg1, arg2, ..., argn): some_code() more_code() return result

As was the case with if and for, we have a block of indented code (the body of the function) in which we can do what we want (e.g. create new variables, call other functions, use if and for). Most functions foo should in the end give back some result, which we do using the return keyword.

Let's see a function with more than one argument:

def bigger_than(x,y): if x > y: return True else: return False
bigger_than(5,6)

Exercise

Define the following functions and test them in some examples.

  • the function sign returning the sign of a real number

  • the function conjugate, taking n×nn \times n-matrices A,SA,S and returning SAS1S A S^{-1}

Solution (uncomment to see)

We can also have functions with optional arguments. These are arguments which can be given, but don't have to be given when calling the function. However, in exchange you need to provide a default value for each of them. Let's generalize the function parabola from above to implement the function xax2+cx \mapsto a \cdot x^2 + c:

def parabola(x, a=1, c=0): x_squared = x^2 return a*x_squared + c
parabola(2, 5, 1) # returns 5 * 2^2 +1
parabola(3)
parabola(3,-1)

As we see above, if a function gets a number of arguments which is smaller than the total number of fixed and optional arguments, it simply fills the optional arguments in the order that was specified in the function definition.

If we want to explicitly specify the argument c above, but leave a to take the default value, we can do so by explicitly specifying the name of the argument that we want:

parabola(3, c=-1)

The general pattern looks like this:

def foo(arg1, arg2, ..., argn, opt_arg1 = default1, ..., opt_argm = defaultm): some_code() more_code() return result

Here, it is important that the non-optional arguments come before the optional ones.

Exercise

Implement a function fibonacci computing the Fibonacci-Numbers FnF_n defined by F0=0,F1=1,Fn=Fn1+Fn2 for n2. F_0 = 0, \quad F_1 = 1, \quad F_{n} = F_{n-1} + F_{n-2} \text{ for }n \geq 2. Add the option to specify F0,F1F_0, F_1 by hand.

Solution (uncomment to see)

Calculus

SageMath can do many computations related to the infinitesimal behaviour of functions, such as limits, derivatives and integrals. In each case, the function should be represented by a symbolic expression, as seen in the last lecture. Since it is useful for illustrating the later computations (and in general), we'll start however with some information how to plot functions (a.k.a. producing pretty pictures).

Plotting

The most basic task here is to plot a function in one variable on an interval. For instance, we can give the function as a symbolic expression in x:

plot(sin(x),x,-5,5)

Alternatively, we can program a function as seen in the first part of the lecture:

def g(x): if x <= 0: return cos(x) else: return floor(x) # floor is the round-down function, sending 3.14 to 3, etc.
plot(g,x,-5,5)

Most of the time it's more useful to use the symbolic expressions, though.

A curiosity: if you type plot(g(x),x,-5,5) it does not do what you want. If you understand what goes wrong there, you already have a pretty good grasp on symbolic expressions!

Exercise

  • Plot your favorite function on an interval of your choice.

  • Plot your second-favorite function on the interval (3,3)(-3,3), but in green.
    Hint: Since unfortunately I forgot to say how to change the color, you will have to find it out yourself!

Solution (uncomment to see)

As you can guess, there are lots more optional arguments for the plot function, which you can again find in the documentation.

plot(sign, x, -3, 3, ymin=-2, ymax=2, axes=False, fill='axis')

The function plot actually returns a Graphics-object and you can literally add multiple such objects to plot multiple functions:

f = plot(sin(x),x,-6,6) g = plot(sin(2*x),x,-6,6,color='red') f+g

Three more variants of plots:

var('y') implicit_plot(x^2+y^2-1, (-3,3), (-2,2)) # plots the solution set of x^2 + y^2 -1 = 0 inside [-3,3] x [-2,2]
var('t') parametric_plot((3*cos(t), 2*sin(t)), (t,0,2*pi)) # plots the curve t -> (3 cos(t), 2 sin(t)) for t in [0, 2 pi]
plot3d(sin(x) + e^cos(y),(x,-3,3),(y,-3,3))

Exercise

Consider the following system of equations: 7x211xy+14y2+4x15y51=0,49x2+59xy94x43y+159=0 7 \, x^{2} - 11 \, x y + 14 \, y^{2} + 4 \, x - 15 \, y - 51 = 0\,,\\ -49 \, x^{2} + 59 \, x y - 94 \, x - 43 \, y + 159 = 0

  • Plot the two solution sets to the two equations in separate colors (in the same coordinate system).

  • Solve this system of equations (using the methods from last lecture) and verify graphically that you found all solutions in the first part of the exercise.

To spare you the typing, here are the left-hand sides of the equations:

var('x,y') F1 = 7*x^2 - 11*x*y + 14*y^2 + 4*x - 15*y - 51 F2 = -49*x^2 + 59*x*y - 94*x - 43*y + 159

Solution (uncomment to see)

Limits

For a function given by a symbolic expression, we can compute its limit at a given point as follows:

limit((x^2 - 1)/(x - 1), x=1)

Exercise

Compute the limit limx0sinh(x)x \lim_{x \to 0} \frac{\sinh(x)}{x} and check with a plot that this is plausible.

Solution (uncomment to see)

We can also deal with divergent functions, one-sided limits and limits as xx goes to ±\pm \infty:

limit(1/x, x=0, dir='+')
limit((x+1)/x, x=Infinity)

Derivatives and integrals

Of course SageMath can differentiate (symbolic) functions, and for many functions it can also integrate them. Here are three ways to obtain the derivative:

f = x * exp(x) print(diff(f,x)) print(derivative(f,x)) print(f.derivative(x))

To get the derivative at a specific point, we can evaluate the symbolic expression there, as seen in the last lecture:

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

Exercise

Write a function plot_tangent which takes as input

  • a symbolic expression f of a function in the variable x

  • a point x0

and outputs a plot of the function f on the interval [x03,x0+3][x_0-3, x_0+3] together with its tangent line at x0.
Bonus: If you are done early, feel free to include optional parameters to change the interval of plotting, the colour of the tangent line etc.

Solution (uncomment to see)

Of course we can also compute partial derivatives for functions depending on multiple parameters:

var('a,b') derivative(a^2*b^3,a)

Here is an example of how to do indefinite integration:

g = integral(sin(x)*cos(x),x); g

The syntax for a definite integral like 02πg(x)dx \int_0^{2 \pi} g(x) dx is as follows:

integral(g, (x,0,2*pi))

To apply the things we learned, here is a task from an actual Analysis II-exam at ETH Zurich (see Exercise 3).

Exercise

Let a,b>0a, b > 0. Let the function f:R2Rf : \mathbb{R}^2 \to \mathbb{R} be defined by f(x,y):=(ax2+by2)e(x2+y2). f(x,y) := (a x^2 + b y^2) \cdot e^{-(x^2 + y^2)}. Find the global maxima and minima of ff for

  • a>ba>b.

  • a=b=1a=b=1.

Hint: The function ff slopes down to zero towards infinity.

In case you have forgotten how these extreme value problems work, you can of course look it up on wikipedia.

Solution (uncomment to see)

Finally, we can also do some calculus involving vectors. Here is the parametrization r:[0,2π]R3r : [0, 2 \pi] \to \mathbb{R}^3 of a spiral:

var('t') r=vector([cos(t), sin(t), t]); r
parametric_plot3d(r, (t, 0, 2*pi))
v = derivative(r,t); v

The arc length of the path rr given by (r)=02πr˙(t)dt. \ell(r) = \int_{0}^{2 \pi} | \dot r(t)| dt.

integrate(v.norm(),t,0, 2*pi)

There are many more interesting functions, and we don't have time to go through all of them. In particular, see this tutorial on vector calculus in SageMath. There, you learn how to define a vector field, compute its divergence and curl, etc.

Power series

Power series are a very useful part of mathematics. Typically one first encounters them in the form of the Taylor series k=0f(k)(0)k!tk \sum_{k=0}^\infty \frac{f^{(k)}(0)}{k!} t^k of a smooth function ff around the point t0=0t_0=0. However, it is also possible to do calculations with more general series k=0aktk. \sum_{k=0}^\infty a_k t^k. The easiest way to compute with them in SageMath is to create a PowerSeriesRing, with some formal variable t. Here you have to specify a base ring, which is the ring containing the possible coefficients aka_k above.

R.<t> = PowerSeriesRing(QQ); R

Then, to compute the Taylor series of the function ff at 00, we can just plug the variable t of our power series ring into the function f.

P = exp(t); P

By default, this computes it to order t19t^{19}, but this can be changed as follows:

R.<t> = PowerSeriesRing(QQ, default_prec=7) P = exp(t); P

Note that default_prec=n means the last coefficient that is computed is the one for tn1t^{n-1}, which can be somewhat confusing ...

The nice thing about power series is that we can add and multiply them, and the Taylor series of fgf \cdot g is equal to the product of the series for ff and gg:

Q = sin(t); Q
P*Q
taylor(exp(x) * sin(x), x, 0, 6)

To get the coefficient of tkt^k of such a power series PP, we can type P[k]:

P[3]

A cool application of power series are so-called generating functions. Given a sequence (ak)k0(a_k)_{k \geq 0} of real numbers, we say that a function f(t)f(t) is a generating function for this sequence, if k0aktk\sum_{k \geq 0} a_k t^k is the Taylor series of ff at t0=0t_0 = 0.

Exercise

The Wikipedia page of the Fibonacci numbers claims that the function f(t)=t1tt2 f(t) = \frac{t}{1-t-t^2} is a generating function for these numbers. Check below if the first coefficients of the Taylor series of ff are correct, and use this method to compute the 4242nd Fibonacci number.

Solution (uncomment to see)

Assignments

Exercise (short)

Solve Exercise 2 b) from the Analysis II-exam above.
Remark: I couldn't get SageMath to solve Exercise 2 a) since it involves some third roots, which are kind of tricky.

Solution (uncomment to see)

Exercise

A partition of a natural number n0n \geq 0 is one representation of nn as an ordered sum of positive integers. For instance, the number n=5n=5 has 77 different partitions: 5,4+1,3+2,3+1+1,2+2+1,2+1+1+1,1+1+1+1+1 5,\quad 4+1,\quad 3+2,\quad 3+1+1,\quad 2+2+1,\quad 2+1+1+1,\quad 1+1+1+1+1 If p(n)p(n) denotes the number of partitions of nn, then the infinite product j=111tj=11t11t211t3 \prod_{j=1}^\infty \frac{1}{1-t^j} = \frac{1}{1-t} \cdot \frac{1}{1-t^2} \cdot \frac{1}{1-t^3} \cdots is a generating function of p(n)p(n).

  • (theoretical) Use the formula for the geometric series to show that the finite product j=1m11tj=11t11tm \prod_{j=1}^m \frac{1}{1-t^j} = \frac{1}{1-t} \cdots \frac{1}{1-t^m} is a generating series for the number p(n,m)p(n,m) of partitions of nn with maximal part at most mm.
    In particular, since for mnm \geq n we have p(n,m)=p(n)p(n,m)=p(n), this explains why the infinite product above is a generating series for p(n)p(n).

  • Use this knowledge to write a function partition_number(n, m=None) to compute the numbers p(n,m)p(n,m).
    Remark: In the above case, it is useful to set the optional parameter m to the Python value None as default. Then, at the beginning of the function, you should check whether m is None and if yes give it a reasonable value.

  • Check your function in a few cases (e.g. n=5n=5 using the enumeration of partitions of 55 above).

Solution (uncomment to see)