# Lecture 20: Computing with Functions¶

List comprehensions in Python take the form
`[f(x) for x in L]`

where `L`

is some list and `f`

some function.
The result is a new list that contains all function values of `f(x)`

for all elements `x`

in the list `L`

. In SageMath we can apply list
comprehensions to build expressions of arbitrary size and shape.

## List Comprehensions¶

We can use a list comprehension
to create a sequence of 20 variables,
starting from `x01`

, `x02`

, .., `x20`

.

```
v = [var('x' + '%02d' % k) for k in range(1,21)]
v
```

Observe the string concatenation with the proper formatting
of the integer index. This leads to the list
```
[x01, x02, x03, x04, x05, x06, x07, x08, x09, x10, x11, x12,
x13, x14, x15, x16, x17, x18, x19, x20]
```

.
Now we raise every variable to the power three.

```
p = [x^3 for x in v]
p
```

and we computed
```
[x01^3, x02^3, x03^3, x04^3, x05^3, x06^3, x07^3, x08^3, x09^3, x10^3,
x11^3, x12^3, x13^3, x14^3, x15^3, x16^3, x17^3, x18^3, x19^3, x20^3]
```

.
We add up the powers into a sum to make an expression.

```
e = sum(p)
print(e, 'has type', type(e))
```

The sum leads to the expression
```
x01^3 + x02^3 + x03^3 + x04^3 + x05^3 + x06^3 + x07^3 + x08^3 + x09^3
+ x10^3 + x11^3 + x12^3 + x13^3 + x14^3 + x15^3 + x16^3 + x17^3 + x18^3
+ x19^3 + x20^3
```

.
We can return to the list representation via the `operands()`

method.

```
ope = e.operands()
print(ope, 'has type', type(ope))
```

and this gives a list
```
[x01^3, x02^3, x03^3, x04^3, x05^3, x06^3, x07^3, x08^3, x09^3, x10^3,
x11^3, x12^3, x13^3, x14^3, x15^3, x16^3, x17^3, x18^3, x19^3, x20^3]
```

.
Raising the k-th operand to the power k goes as follows.

```
r = [ope[k]^k for k in range(len(ope))]
r
```

and then we see ```
[1, x02^3, x03^6, x04^9, x05^12, x06^15, x07^18, x08^21,
x09^24, x10^27, x11^30, x12^33, x13^36, x14^39, x15^42, x16^45, x17^48,
x18^51, x19^54, x20^57]
```

.
Note that all operands are expressions.

```
print(r[2], 'has type', type(r[2]))
print('its degree in' , v[2], ':', r[2].degree(v[2])
```

With the proper selection of the variable in the argument
of `degree()`

we get the right degree, which is `6`

for `r[2]`

.
We can compute the degrees in their variables.

```
[(r[k]).degree(v[k]) for k in range(len(ope))]
```

which leads to ```
[0, 3, 6, 9, 12, 15, 18, 21, 24, 27, 30,
33, 36, 39, 42, 45, 48, 51, 54, 57]
```

.

Suppose we want to select those operands of degree less than 13. Recall the unit_step() function.

```
print(r[2], ':', unit_step(13 - r[2].degree(v[2])))
print(r[8], ':', unit_step(13 - r[8].degree(v[8])))
```

and we see `x03^6 : 1`

and `x09^24 : 0`

.
Now we can remove all operands with degree higher than 13.

```
filter = lambda x, k: x*unit_step(13 - x.degree(v[k]))
s = [filter(r[k],k) for k in range(len(r))]
s
```

and we see ```
[1, x02^3, x03^6, x04^9, x05^12,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
```

and then build the expression again with `sum()`

.

```
t = sum(s)
t
```

and we obtain `x05^12 + x04^9 + x03^6 + x02^3 + 1`

.

## Composing Functions¶

Consider the trajectory of a projectile in space modeled by a parabola, subject to the following constraints. At time t = 0 it is launched from the ground and at time t = 45 it hits the ground 120 miles further. Assuming constant horizontal speed we create a function f(t) to give the altitude of the projectile in function of time. First we model the shape of the trajectory.

```
y(x) = x*(120 - x)
plot(y(x), (x, 0, 120), aspect_ratio=1/100, thickness=3, color='red')
```

The parameters `thickness`

and `color`

embellish the plot,
adjusting the default value for the `aspect_ratio`

to `1/100`

alters the display of the shape of the trajectory,
so we do not have to worry about units.
The plot is displayed in Fig. 24.

Now we want a function in time t. The assumption of constant horizontal speed implies that x is just a rescaling of t. This means that when t = 45, x must be 120.

```
x(t) = 120/45*t
print('x at 0 :', x(0))
print('x at 45 :', x(45))
```

and indeed, we see `x at 0 : 0`

and `x at 45 : 120`

printed.
Then the altitude is given as the composition of y after x.

```
f(t) = y(x(t))
print('altitude at 0 :', f(0))
print('altitude at 22.5 :', f(22.5))
print('altitude at 45 :', f(45))
```

as altitudes at 0, 22.5, and 45, we respectively see `0`

,
`3600.00000000000`

and `0`

.
With the composition of functions we have separated the geometric
shape of the trajectory and its evolution in time.
Suppose we want to halve the speed.

```
h(t) = t/2
hf(t) = f(h(t))
print('altitude at 45 :', hf(45))
print('altitude at 90 :', hf(90))
```

and the prints confirm that the projectile reaches its peak at 45 and falls to the ground at 90.

We can also make a 3d-plot of the space curve.

```
parametric_plot3d([x(t), 0, f(t)], (t, 0, 45), thickness=5, color='red')
```

and this shows the following image:

## Functions Returning Functions¶

We can define functions that return functions. Suppose we are interested in solving equations with Newton’s method.

```
def newton_step(f, x):
"""
Returns the function that will perform one Newton step
to solve the equation f(x) = 0.
"""
n(x) = x - f(x)/diff(f,x)
return n
```

Let us try it out on the square function, so we can compute the square root of two.

```
sqr(x) = x^2 - 2
our_sqrt = newton_step(sqr, x)
our_sqrt
```

Then we see `x |--> x - 1/2*(x^2 - 2)/x`

.
Observe that if we start at an integer value,
we get rational approximations for the square root of two.

```
our_sqrt(our_sqrt(our_sqrt(2)))
```

Starting at `2`

, three iterations of Newton’s method yields `577/408`

.
For a repeated application of a function, we multiply strings.

```
s = 'our_sqrt('*5 + '2' + ')'*5
s
```

We will then execute
`our_sqrt(our_sqrt(our_sqrt(our_sqrt(our_sqrt(2)))))`

and we evaluate the expression.

```
v = eval(s)
v
```

This gives `886731088897/627013566048`

.
To compare with an actual approximation for sqrt(2).

```
print(v.n(digits=30))
print(sqrt(2).n(digits=30))
```

and we see

```
1.41421356237309504880168962350
1.41421356237309504880168872421
```

## Assignments¶

Type

`L = [randint(0,99) for i in range(10)]`

to generate a list of 10 random numbers between 0 and 99. Give the SageMath commands for the following operations.divide every element in the list by 100;

convert the list to a list of floating-point numbers of a precision of 3 decimal places;

select all elements in the list that are larger than 0.5;

compute the sum of the elements in the list.

The command

`is_prime()`

applied to a number returns`True`

if the number is prime and`False`

otherwise. Use`is_prime()`

to make a list of all primes less than 1000. How many 3-digit primes are there?Do

`P.<x> = PolynomialRing(RR, sparse=False)`

and then`q = P.random_element(degree=10)`

.Select from

`q`

all terms with positive coefficients and make a new polynomial with those terms.For a positive integer \(n\), and a list \(X\) of \(n\) variables \(X\) =

`[`

`x1`

,`x2`

, .. ,`xn`

`]`

, define the function\[H(n,X,k,c) = 2 n X[k] - c X[k] \left( 1 + \sum_{i=0}^{n-1} \frac{i+1}{i+3} X[i] \right) - 2 n,\]where \(c\) is some variable parameter. Show the definition of your function is correct for a list of 9 variables.

Let \(n\) and \(k\) be positive natural numbers with \(k < n\) and consider the polynomial \({\displaystyle \sum_{i=0}^{n-1} \prod_{j=0}^{k-1} z_{(i+j) \mbox{ mod } n}}.\)

Define a function that takes on input a list of variables (where \(n\) is the length of the list) and a value for \(k\). On output is the expression as defined by the polynomial above. Hint: use

`prod`

to make the product of a list of expressions.Show the result of your function on a sequence of 10 variables and for \(k = 3\).

Execute the

`newton_step`

function on the`our_sqr`

function, defined again as`x^2 - 2`

. Newton’s method has the property of converging quadratically, that is: in each step the correct number of decimal places doubles. Starting at`2`

, perform 10 steps with Newton’s method and verify the progess of the accuracy of the rational approximations for the square root of 2? Do you observe quadratic convergence?Define a Python function

`makeHorner`

which takes on input the coefficients of a polynomial and the symbol of the variable in the polynomial. For example, if`C`

is`[c0, c1, c2, c3, c4]`

and the variable is`x`

, then`makeHorner(C, x)`

returns(((c4*x + c3)*x + c2)*x + c1)*x + c0

Your function

`makeHorner`

should directly build the polynomial and not apply the`horner()`

to an expanded polynomial.