**Path:**Introduction to computer algebra and mathematical programming/04 Analysis and power series.ipynb

**Views:**

^{42}

**Visibility:**Unlisted (only visible to those who know the link)

**Image:**ubuntu2004

**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 $n$ 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 numberthe function

`conjugate`

, taking $n \times n$-matrices $A,S$ and returning $S 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 $x \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 $F_n$ defined by $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 $F_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)$, 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: $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 $\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 $x$ 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 $[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 $\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 > 0$. Let the function $f : \mathbb{R}^2 \to \mathbb{R}$ be defined by $f(x,y) := (a x^2 + b y^2) \cdot e^{-(x^2 + y^2)}.$ Find the global maxima and minima of $f$ for

$a>b$.

$a=b=1$.

*Hint:* The function $f$ 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 \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 $r$ given by $\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 $\sum_{k=0}^\infty \frac{f^{(k)}(0)}{k!} t^k$ of a smooth function $f$ around the point $t_0=0$. However, it is also possible to do calculations with more general series $\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 $a_k$ above.

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

Then, to compute the Taylor series of the function $f$ at $0$, 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 $t^{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 $t^{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 $f \cdot g$ is equal to the product of the series for $f$ and $g$:

Q = sin(t); Q

P*Q

taylor(exp(x) * sin(x), x, 0, 6)

To get the coefficient of $t^k$ of such a power series $P$, we can type `P[k]`

:

P[3]

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

#### Exercise

The Wikipedia page of the Fibonacci numbers claims that the function $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 $f$ are correct, and use this method to compute the $42$nd 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 $n \geq 0$ is one representation of $n$ as an ordered sum of positive integers. For instance, the number $n=5$ has $7$ different partitions: $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)$ denotes the number of partitions of $n$, then the infinite product $\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)$.

(theoretical) Use the formula for the geometric series to show that the finite product $\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)$ of partitions of $n$ with maximal part at most $m$.

In particular, since for $m \geq n$ we have $p(n,m)=p(n)$, this explains why the infinite product above is a generating series for $p(n)$.Use this knowledge to write a function

`partition_number(n, m=None)`

to compute the numbers $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=5$ using the enumeration of partitions of $5$ above).

**Solution** (uncomment to see)