# Function Representation and Manipulation


```
%matplotlib inline
```


```
import numpy as np
import matplotlib
matplotlib.rcParams.update({'font.size': 14})
import matplotlib.pyplot as plt
```

From a mathematical point of view, a central point in numerical methods is how we represent a general function $f(x)$. As $x$ is a real number the function encodes an (uncountably) infinite amount of information: its value $f(x)$ for every point $x$. A computer can only usefully store and manipulate a finite amount of information, so it would appear impossible to represent even the simplest function of one variable.

Of course, most functions can be represented straightforwardly using little information. We often see functions that can be represented symbolically - $f(x) = x^2 + x$ is a simple example - but numerical methods are precisely for cases where symbolic methods are not applicable. Instead we must consider how a function is represented, or approximated, in terms of a finite amount of information.

Two standard representations can be used to make the central point. The first representation is in terms of polynomials, using Taylor series:

\begin{equation}
  f(x ; x_0) = \sum_{n=0}^{\infty} \frac{c_n}{n!} (x - x_0)^n.
\end{equation}

Here $x_0$ is a parameter about which the Taylor series is centred.

The second representation is in terms of trigonometric functions, using Fourier series. For simplicity we consider only the Fourier cosine series

\begin{equation}
  f(x ; L) = \frac{1}{2} a_0 + \sum_{n=1}^{\infty} a_n \cos \left( \frac{n \pi x}{L} \right).
\end{equation}

The Fourier series representation is $2L$-periodic.

Given a set of coefficients ($\{c_n\}$ for the Taylor series or $\{a_n\}$ for the Fourier series) we define a function.


```
an = np.zeros((10, 1))
cn = np.zeros_like(an)
for n in range(len(an)):
    an[n] = 1.0 / (n + 1)**2
    cn[n] = an[n]
x = np.linspace(-1.0, 1.0, 100)
taylor = an[0]*np.ones_like(x)
fourier = cn[0]*np.ones_like(x)
for n in range(1, len(an)):
    taylor += an[n] * x**n
    fourier += cn[n] * np.cos(n * np.pi * x)
fig = plt.figure(figsize = (8,6))
ax = fig.add_subplot(111)
ax.plot(x,taylor,label = 'Taylor')
ax.plot(x,fourier, label = 'Fourier')
plt.legend()
ax.set_xlabel('$x$')
plt.show()
```

Given a function, the standard analytical methods will give you the coefficients. This can rarely be directly applied in numerical methods - as noted above, it is usual that you're dealing with a function that isn't known in symbolic form.

Above, when we defined the function $f(x)$ from the coefficients, we were noting that we could evaluate the representation for any point $x$ given 

1. the form of the representation (e.g., as a Taylor series about $x_0=0$), and
2. the value of the coefficients.

We can (nearly always) go in the opposite direction: given the value of the function at enough points $x_n$ we can compute the coefficients of the representation.

For example, let us assume that 

1. the function will be represented as a polynomial of order 2, i.e. a Taylor series about $x_0=0$ such that all coefficients $c_n$ vanish for $n > 2$, and
2. the value of the function is known at three points $x_1, x_2, x_3$.

This is precisely the information needed to compute the coefficients $c_0, c_1, c_2$ and hence write the function as

\begin{equation}
  f(x) \approx c_0 + c_1 x + \frac{c_2}{2} x^2.
\end{equation}

Explicitly, we note that the function representation matches the known value of the function at the three points only if

\begin{equation}
  \begin{pmatrix} 1 & x_1 & \tfrac{1}{2} x_1^2 \\ 1 & x_2 & \tfrac{1}{2} x_2^2 \\ 1 & x_3 & \tfrac{1}{2} x_3^2 \end{pmatrix} \begin{pmatrix} c_0 \\ c_1 \\ c_2 \end{pmatrix} = \begin{pmatrix} f(x_1) \\ f(x_2) \\ f(x_3) \end{pmatrix}
\end{equation}

which, given the information we know, is a linear system for the coefficients.

Note in particular that the number of points at which the function is known must match the number of coefficients in the function representation.

The previous example hints at one of the standard methods of numerically representing functions. First, we take a representation that we know, given an infinite amount of information, converges to the function we want (in some sense). Second, assume that taking a *finite* number of terms is sufficient: e.g., truncate the sum in the Taylor series expansion after a finite number of terms. Then, using a finite amount of information about the function (such at its value at an appropriate number of points) we compute the coefficients of the representation. We then manipulate the representation *as if* it were the true function.

## From data to derivative

A standard example of manipulating "functions" would be the following. The *value* of the function $f$ is known at three points $x_{1,2,3}$ with $x_1 < x_2 < x_3$. Find the derivative of the function at some location in $[x_1, x_3]$.

One method is to use finite differences. However, it is instructive to go through the steps in terms of function representations.

Here we will use the Taylor series representation above. That is, we assume that $f(x) \simeq g(x)$, where $g(x)$ is the *unique* polynomial of degree 2 that matches the known values of $f$, i.e. $g(x_i) = f(x_i)$ for $i = 1, 2, 3$, and

$$
  g(x) = c_0 + c_1 x + \frac{c_2}{2} x^2
$$

as above. For simplicity we'll assume an evenly spaced grid of points, i.e. $x_1 = x_2 - h$ and $x_3 = x_2 + h$.


```
import sympy
from sympy.matrices import *
x1,x2,x3=sympy.symbols(("x1","x2","x3"))
h = sympy.symbols("h")
M = Matrix([[1, x1, x1**2/2],[1, x2, x2**2/2],[1, x3, x3**2/2]])
M = M.subs(x1, x2-h)
M = M.subs(x3, x2+h)
f1,f2,f3=sympy.symbols(("f1","f2","f3"))
b = Matrix([f1, f2, f3])
c = M.solve(b)
sympy.pretty_print(c.expand())
```

This gives us the $c_i$ coefficients. Now, as $f(x) \simeq g(x)$ we *assume* that $f'(x) \simeq g'(x)$. As we have the explicit form for $g$, we immediately know

$$
  \begin{align}
    && g(x) & = c_0 + c_1 x + \frac{c_2}{2} x^2 \\
    \implies && g'(x) & = c_1 + c_2 x \\
    && & = \frac{1}{h} \left( \frac{f_3 - f_1}{2} + \frac{x - x_2}{h} \left( f_1 - 2 f_2 + f_3 \right) \right)
  \end{align}
$$

We immediately see that 

$$
  f'(x_2) \simeq \frac{f_3 - f_1}{2 h}
$$

which is the standard central differencing formula.

Let us try it for a simple example - $f(x) = \exp(x)$ on $[0, 1]$. We see that the derivative is reasonable, but varies linearly with $x$. This is to be expected - $g(x)$ is a quadratic, so its derivative must be a linear polynomial.


```
x = np.array([0.0, 0.5, 1.0])
xp = np.linspace(0.0, 1.0) # Used for plotting
h = 0.5
f = np.exp(x) # Samples of f
df = np.exp(xp) # The derivative of f is again the exponential; this is for plotting
dg = 1/h * ((f[2] - f[0])/2.0 + (xp - x[1]) / h * (f[0] - 2.0*f[1] + f[2]))
fig = plt.figure(figsize = (8,6))
ax = fig.add_subplot(111)
ax.plot(xp, df, 'b-', label=r"$\exp(x)$")
ax.plot(xp, dg, 'r--', label=r"$g'(x)$")
ax.set_xlabel(r"$x$")
plt.legend(loc=2);
```

We should note that there is a distinction between what we have done above and the standard approach to finite differencing. In finite differencing we evaluate the result at a *single* point. Here, as we have the functional form of the approximation $g(x)$ at all points in the interval, we can evaluate the result *anywhere*. 

## From data to integral

If we wanted to integrate the function $f(x)$ over the interval, a simple approximation would be to integrate the approximate function $g(x)$. This is the basis for Simpson's rule, in the case of a second order polynomial interpolating function. In the general case, where the interpolating function is a (piecewise) polynomial of some order, it is the basis of the Newton-Cotes methods.

By using the form of $g(x)$ above we have

$$
  \begin{align}
    && g(x) & = c_0 + c_1 x + \frac{c_2}{2} x^2 \\
    \implies && \int_0^{2h} g(x) \, dx & = 2 h c_0 + 2 h^2 c_1 + \frac{4}{3} h^3 c_2
  \end{align}
$$

which we can evaluate (as $x_2 = h$, using the results for $c_i$ above)

$$
\begin{equation}
  \int_0^{2 h} g(x) \, dx =  \frac{h}{3} \left( f_1 + 4 f_2 + f_3 \right).
\end{equation}
$$

This is precisely Simpson's rule.

## From equations to algorithms

Suppose we want to solve a differential equation $y' = f(x, y)$.

We can *assume* that the solution that we want can be approximated by a simple form; again, let us assume that in the interval $[x_n, x_{n+1}]$ it can be approximated by a linear polynomial. Therefore $y(x) \simeq Y(x)$ where

$$\begin{equation}
  Y(x) = c_0 + c_1 (x - x_{n}).
\end{equation}$$

This immediately tells us that $y' \simeq Y' = c_1$, which is constant across the entire interval $[x_n, x_{n+1}]$.

Now, if we assume that we *know* the value of the solution at the start of the interval, $y_n = y(x_n)$, then we can insist that

$$\begin{equation}
  Y(x_n) = c_0 = y_n.
\end{equation}$$

As we also *know*, given the value of $y_n$, that

$$\begin{equation}
  y'(x_n) = f(x_n, y_n) \simeq Y'(x_n) = c_1,
\end{equation}$$

we have fixed all the parameters determining $Y(x)$.

Therefore we can *evaluate*

$$\begin{equation}
  y_{n+1} = y(x_{n+1}) \simeq Y(x_{n+1}) = c_0 + c_1 (x_{n+1} - x_n) = y_n + (x_{n+1} - x_n) f(x_n, y_n).
\end{equation}$$

If we denote $h = (x_{n+1} - x_n)$ then we see the standard formula for *Euler's method* for solving initial value problems

$$\begin{equation}
  y_{n+1} =  y_n + h f(x_n, y_n).
\end{equation}$$

8fa717b1b5302d6bfdf8183210c45eb0c4e3be30
ipynb
Jupyter Notebook
Lectures/Function Representation and Manipulation.ipynb
josh-gree/NumericalMethods
03cb91114b3f5eb1b56916920ad180d371fe5283
CC-BY-3.0
2015-02-12T19:51:52.000Z
2022-03-26T15:34:11.000Z
Lectures/Function Representation and Manipulation.ipynb
josh-gree/NumericalMethods
03cb91114b3f5eb1b56916920ad180d371fe5283
CC-BY-3.0
2017-05-24T19:49:52.000Z
2018-01-23T21:40:42.000Z
Lectures/Function Representation and Manipulation.ipynb
josh-gree/NumericalMethods
03cb91114b3f5eb1b56916920ad180d371fe5283
CC-BY-3.0
2015-01-05T13:30:47.000Z
2022-02-15T09:59:39.000Z
Qwen/Qwen-72B
1. YES
2. YES