- \n",
"
- Overview and Learning Objectives \n", "
- Power Series \n", "
- Exponentiation Identity and Euler's Formula \n", "
- Differential Equations \n", "
- Roots of Unity \n", "
- Exercise 1: Approximation of Exponential Function via Power Series \n", "
- Exercise 2: Gaussian Function \n", "
- Exercise 3: Spiral Generation \n", "

\n",
"## Overview and Learning Objectives

\n",
"\n",
"The exponential function is one of the most important functions in mathematics. In everyday life, we encounter this function when a phenomenon (e.g., the spread of a viral infection) can be modeled by an initial value and a growth rate. The exponential function has several remarkable mathematical properties, which lead to different ways on how to approach and define this concept. In this unit, we introduce the **exponential function** by its **power series**. This definition allows for expanding the definition of a real (defined for real numbers in $\\mathbb{R}$) to a complex exponential function (defined for complex numbers in $\\mathbb{C}$). The complex version of the exponential function will play a central role for defining and understanding the **Fourier transform** covered in Unit 9. We will then go through important properties such as the exponentiation identity and Euler's Formula, which sheds a different, more natural light on the trigonometric identities of the sine and the cosine functions. Furthermore, we discuss the exponential function from the perspective of differential equations. This also leads to numerical methods for approximating the exponential function's values (methods that are much more efficient than using the power series). Finally, we introduce the notion of **roots of unity**, which are the roots of a specific polynomial (being of the form $z^N-1$ for some $N\\in\\mathbb{N}$) and can be expressed in terms of the exponential function. These roots of unity are the building blocks of the **discrete Fourier transform** (DFT) and the **FFT algorithm**—topics we cover in Unit 9. While discussing the exponential function, another goal of this unit is to further deepen your skills in Python programming by applying the concepts learned in previous units. In Exercise 1 we ask you to implement and compare two algorithms to approximate the exponential function for some given argument. Then, in Exercise 2, you will write a Python program to compute and plot the Gaussian function given the exponential function. Finally, you will apply in Exercise 3 the exponential function to create spirals with different properties, which will deepen your understanding of the relationship between the complex exponential function and angles. \n",
" \n",
"

"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" \n",
"## Power Series\n",
"\n",
"One encounters the **real exponential function** $\\exp:\\mathbb{R}\\to \\mathbb{R}$ in the context of many mathematical applications, and the function can be characterized in many different ways. Historically, the exponential function was studied already by **Johann Bernoulli** in the $17^\\mathrm{th}$ century when considering **interest rates**: Assume that an amount of $1$ earns an interest $x$ at an annual rate compounded monthly. Then the interest earned each month is $\\frac{x}{12}$ times the current value, so that each month the total value is multiplied by $\\left(1+\\frac{x}{12}\\right)$ and the value at the end of the year is $\\left(1+\\frac{x}{12}\\right)^{12}$. In case the interest is compounded every day, it becomes $\\left(1+\\frac{x}{365}\\right)^{365}$. Letting the time intervals grow per year by making them shorter leads to the limit definition of the exponential function\n",
"\n",
"$$\\exp(x) = \\mathrm{lim}_{n\\to\\infty} \\left(1+\\frac{x}{n}\\right)^{n},$$\n",
"\n",
"which was first given by **Leonhard Euler**. The constant $e:=\\exp(1)\\approx 2.71828 \\ldots$ is also known as **Euler's number**. By expanding the $n$-fold product in the above definition, one can show that the exponential function can also be expressed by the following power series:\n",
"\n",
"$$\\exp(x) := \\sum_{n=0}^{\\infty} \\frac{x^n}{n!} = 1 + x + \\frac{x^2}{1 \\cdot 2} + \\frac{x^3}{1 \\cdot 2 \\cdot 3} + \\dots$$\n",
"\n",
"with $x\\in\\mathbb{R}$. Replacing in the power series the real-valued variable $x\\in\\mathbb{R}$ by a complex-valued variable $z\\in\\mathbb{C}$, one still obtains the **complex exponential function** $\\exp:\\mathbb{C}\\to \\mathbb{C}$ given by \n",
"\n",
"$$\\exp(z) := \\sum_{n=0}^{\\infty} \\frac{z^n}{n!} = 1 + z + \\frac{z^2}{1 \\cdot 2} + \\frac{z^3}{1 \\cdot 2 \\cdot 3} + \\dots$$\n",
"\n",
"In the following plot, we visualize the real part, the imaginary part, as well as the absolute value of the complex exponential function over the complex plane. It can be seen that the absolute value $|\\exp(z)|$ only depends on the real part $x$ of the complex argument $z=x+iy$, while increasing exponentially with increasing $x$. Furthermore, the real part $\\mathrm{Re}(\\exp(z))$ and imaginary part $\\mathrm{Im}(\\exp(z))$ show periodic oscillations over $y$ for a fixed $x$. This behavior becomes clear from Euler's formula and other trigonometric identities that hold for the exponential function. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from matplotlib import pyplot as plt\n",
"%matplotlib inline\n",
"\n",
"X, Y = np.meshgrid(np.arange(-2, 2, 0.1), np.arange(-12, 12, 0.1))\n",
"Z = X + Y*1j\n",
"f_exp = np.exp(Z)\n",
"\n",
"plt.figure(figsize=(12, 4))\n",
"extent = [-2, 2, -12, 12]\n",
"plt.subplot(1, 3, 1)\n",
"plt.imshow(np.real(f_exp), aspect='auto', cmap='seismic', origin='lower', extent=extent)\n",
"plt.title('Real part Re(exp(z))')\n",
"plt.xlabel('x = Re(z)')\n",
"plt.ylabel('y = Im(z)')\n",
"plt.colorbar()\n",
"plt.subplot(1, 3, 2)\n",
"plt.imshow(np.imag(f_exp), aspect='auto', cmap='seismic', origin='lower', extent=extent)\n",
"plt.title('Imaginary part Im(exp(z))')\n",
"plt.xlabel('x = Re(z)')\n",
"plt.ylabel('y = Im(z)')\n",
"plt.colorbar()\n",
"plt.subplot(1, 3, 3)\n",
"plt.imshow(np.abs(f_exp), aspect='auto', cmap='gray_r', origin='lower', extent=extent)\n",
"plt.title('Absolute value |exp(z)|')\n",
"plt.xlabel('x = Re(z)')\n",
"plt.ylabel('y = Im(z)')\n",
"plt.colorbar()\n",
"plt.tight_layout()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" \n",
"## Exponentiation Identity and Euler's Formula \n",
"\n",
"Based on the power series definition, one may prove two famous formulas of the exponential function that explain many of its properties. The first formula is knowns as **exponentiation identity** and says that \n",
"\n",
"$$\n",
" \\exp(z_1 + z_2) = \\exp(z_1)\\cdot \\exp(z_2)\n",
"$$\n",
"\n",
"for any complex numbers $z_1, z_2\\in\\mathbb{C}$. In particular, this property explains the exponential increase for real arguments. For example, \n",
"\n",
"$$\n",
" \\exp(n) = \\exp(1+1+\\ldots +1) = \\exp(1)^n = e^n\n",
"$$\n",
"\n",
"for $n\\in\\mathbb{N}$. The second famous formula, which is known as **Euler's formula**, relates the values of the exponential function at purely imaginary arguments to trigonometric functions. It states that for the complex (and purely imaginary) number $c = i\\gamma$ with some real-valued $\\gamma\\in\\mathbb{R}$ one has the identity \n",
"\n",
"$$\\mathrm{exp}(i\\gamma) = \\cos(\\gamma) + i\\sin(\\gamma) .$$\n",
"\n",
"Actually, starting with the real sine and cosine functions, one often defines $\\mathrm{exp}(i\\gamma)$ by means of the Euler formula (rather than using the power series). This explains the periodic behavior of the real and imaginary part of $\\exp$ along the imaginary (vertical) axis as shown in the previous figure. The real value $\\gamma$ can be thought of as an angle (given in radians). The following visualization shows how the values $\\mathrm{exp}(i\\gamma)$ change when increasing the angle $\\gamma$ from $0$ to $2\\pi$:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from matplotlib import ticker \n",
"%matplotlib inline\n",
"\n",
"cmap = plt.cm.get_cmap('hsv') # hsv is nice because it defines a circular color map\n",
"\n",
"N = 64\n",
"\n",
"fig = plt.figure(figsize=(12, 4))\n",
"ax1 = fig.add_subplot(1, 3, 1, projection='polar')\n",
"ax2 = fig.add_subplot(1, 3, 2)\n",
"ax3 = fig.add_subplot(1, 3, 3)\n",
"\n",
"for i in range(N):\n",
" gamma = 2 * np.pi * i / N\n",
" c = np.exp(1j * gamma)\n",
" color = cmap(i / N)\n",
" ax1.plot([0, np.angle(c)], [0, np.abs(c)], color=color)\n",
" ax1.plot(np.angle(c), np.abs(c), 'o', color=color)\n",
" ax2.plot(gamma, np.real(c), 'o', color=color)\n",
" ax3.plot(gamma, np.imag(c), 'o', color=color)\n",
" \n",
"ax2.grid()\n",
"ax2.set_xlabel('$\\gamma$ [radians]')\n",
"ax2.set_ylabel('$\\mathrm{Re}(\\exp(i \\gamma))$')\n",
"ax2.xaxis.set_major_formatter(ticker.FormatStrFormatter('$%s$')) \n",
"\n",
"ax3.grid()\n",
"ax3.set_xlabel('$\\gamma$ [radians]')\n",
"ax3.set_ylabel('$\\mathrm{Im}(\\exp(i \\gamma))$')\n",
"ax3.xaxis.set_major_formatter(ticker.FormatStrFormatter('$%s$')) \n",
"plt.tight_layout()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We have seen from Euler's formula that the complex values $\\mathrm{exp}(i\\gamma)$ lie on the unit circle of the complex plane for all $\\gamma\\in\\mathbb{R}$. Furthermore, due to periodicity, it suffices to consider $\\gamma\\in[0,2\\pi)$. In fact, $\\gamma$ encodes the angle (in radians) of the complex number $c = \\mathrm{exp}(i\\gamma)$, while $|c|=1$. From the exponentiation identity and Euler's formula, one can derive the following properties of the exponential function:\n",
"\n",
"\\begin{eqnarray}\n",
"\\exp(i\\gamma) & = & \\exp(i(\\gamma+2\\pi)) \\\\\n",
"|\\exp(i\\gamma)| & = & 1 \\\\\n",
"\\overline{\\exp(i\\gamma)} & = & \\exp(-i\\gamma) \\\\\n",
"\\exp(i(\\gamma_1+\\gamma_2)) & = & \\exp(i\\gamma_1) \\exp(i\\gamma_2) \\\\\n",
"\\end{eqnarray}\n",
"\n",
"Plugging in Euler's formula in the last identity, on obtains the trigonometric identities for the sine and cosine:\n",
"\n",
"\\begin{eqnarray}\n",
" \\cos(\\gamma_1+\\gamma_2) &=& \\cos(\\gamma_1)\\cos(\\gamma_2) - \\sin(\\gamma_1)\\sin(\\gamma_2)\\\\\n",
" \\sin(\\gamma_1+\\gamma_2) &=& \\cos(\\gamma_1)\\sin(\\gamma_2) + \\cos(\\gamma_2)\\sin(\\gamma_1)\n",
"\\end{eqnarray}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" \n",
"## Differential Equations\n",
"\n",
"The exponential function can be characterized by another important property in terms of differential equations. Let us consider the differential equation $\\frac{df}{dx}(x)=f(x)$ with initial condition $f(0)=1$. The [**Picard–LindelĂ¶f Theorem**](https://en.wikipedia.org/wiki/Picard%E2%80%93Lindel%C3%B6f_theorem) implies that there exists a unique solution. Using the power series definition, one can easily check that the exponential function indeed fulfills these properties. The differential equation also holds for the complex exponential function. In particular, one obtains the following equations:\n",
"\n",
"\\begin{eqnarray}\n",
"\\frac{d\\exp(x)}{dx} & = & \\exp(x)\\\\\n",
"\\frac{d\\exp(z)}{dz} & = & \\exp(z)\\\\\n",
"\\frac{d\\exp(i\\gamma)}{d\\gamma} & = & i\\exp(i\\gamma)\n",
"\\end{eqnarray}\n",
"\n",
"There are many numerical methods for approximating solutions of differential equations with initial conditions. The easiest method is known as **Euler's method**, where one starts with the initial value and then uses a tangent line over a short step size to estimate the function's value at the next step. In the next code cell, we implement this procedure for the case of the real exponential function. The figure shows the approximative solution compared to the NumPy function `np.exp`. To better understand the quality of the solution of the interval considered, we apply `plt.semilogy` for plotting."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def exp_approx_Euler(x_min=0, x_max=2, x_delta=0.01, f_0=1):\n",
" \"\"\"Approximation of exponential function using Euler's method\n",
"\n",
" Notebook: PCP_07_exp.ipynb\n",
"\n",
" Args:\n",
" x_min: Start of input interval (Default value = 0)\n",
" x_max: End of input interval (Default value = 2)\n",
" x_delta: Step size (Default value = 0.01)\n",
" f_0: Initial condition (Default value = 1)\n",
"\n",
" Returns:\n",
" f: Signal\n",
" x: Sampled input interval\n",
" \"\"\"\n",
" x = np.arange(x_min, x_max+x_delta, x_delta)\n",
" N = len(x)\n",
" f = np.zeros(N)\n",
" f[0] = f_0\n",
" for n in range(1, N):\n",
" f[n] = f[n-1] + f[n-1]*x_delta\n",
" return f, x\n",
"\n",
"\n",
"plt.figure(figsize=(12, 4))\n",
"x_max = 3\n",
"x_delta = 0.1\n",
"f, x = exp_approx_Euler(x_min=0, x_max=x_max, x_delta=x_delta, f_0=1)\n",
"\n",
"plt.subplot(1, 3, 1)\n",
"plt.plot(x, f, 'r')\n",
"plt.plot(x, np.exp(x), 'k')\n",
"plt.legend(['Approximation','exp'])\n",
"plt.xlim([0, x_max])\n",
"plt.grid()\n",
"plt.title('Approximation with $\\Delta$ = %.1f' % x_delta)\n",
"\n",
"plt.subplot(1, 3, 2)\n",
"plt.semilogy(x, f, 'r')\n",
"plt.semilogy(x, np.exp(x), 'k')\n",
"plt.legend(['Approximation','exp'])\n",
"plt.xlim([0, x_max])\n",
"plt.grid(which='both')\n",
"plt.title('Approximation with $\\Delta$ = %.1f' % x_delta)\n",
"\n",
"\n",
"x_delta = 0.01\n",
"f, x = exp_approx_Euler(x_min=0, x_max=x_max, x_delta=x_delta, f_0=1)\n",
"plt.subplot(1, 3, 3)\n",
"plt.semilogy(x, f, 'r')\n",
"plt.semilogy(x, np.exp(x), 'k')\n",
"plt.legend(['Approximation','exp'])\n",
"plt.xlim([0, x_max])\n",
"plt.grid(which='both')\n",
"plt.title('Approximation with $\\Delta$ = %.2f' % x_delta)\n",
"plt.tight_layout()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" \n",
"## Roots of Unity\n",
"\n",
"Let $N \\in \\mathbb{N}_{>0}$ be a positive integer. A complex number $\\rho \\in \\mathbb{C}$ is called an $N^\\mathrm{th}$ **root of unity** if $\\rho^N = 1$. It is not hard to see that there are exactly $N$ distinct $N^\\mathrm{th}$ roots of unity, which are exactly the $N$ different roots of the polynomial $z^N-1$ (see also Exercise 2 of Unit 6). Additionally, if $\\rho^n \\neq 1$ for all $n\\in [1:N-1]$, the root $\\rho$ is called a **primitive** $N^\\mathrm{th}$ root of unity. With the properties mentioned above, it is easy to see that $\\rho_N:=\\exp(2 \\pi i / N)$ is such a **primitive** $N^\\mathrm{th}$ root of unity. Furthermore, all $N^\\mathrm{th}$ roots of unity can be generated by considering the powers of $\\rho_N$:\n",
"\n",
"$$1=\\rho_N^0, \\quad \\rho_N^1, \\quad \\rho_N^2, \\quad ...,\\quad \\rho_N^{N-1}$$\n",
"\n",
"The following plot shows all roots of unity for different integers $N \\in \\mathbb{N}_{>0}$. The primitive roots are indicated in red. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from math import gcd\n",
"\n",
"def plot_vector(c, color='k', start=0, linestyle='-'):\n",
" \"\"\"Plotting complex number as vector\n",
"\n",
" Notebook: PCP_07_exp.ipynb\n",
"\n",
" Args:\n",
" c: Complex number\n",
" color: Vector color (Default value = 'k')\n",
" start: Start of vector (Default value = 0)\n",
" linestyle: Line Style of vector (Default value = '-')\n",
" \"\"\"\n",
" return plt.arrow(np.real(start), np.imag(start), np.real(c), np.imag(c),\n",
" linestyle=linestyle, head_width=0.05,\n",
" fc=color, ec=color, overhang=0.3, length_includes_head=True)\n",
"\n",
"\n",
"def plot_root_unity(N, ax):\n",
" \"\"\"Plotting N-th root of unity into figure with axis\n",
"\n",
" Notebook: PCP_07_exp.ipynb\n",
"\n",
" Args:\n",
" N: Root number\n",
" ax: Axis handle\n",
" \"\"\"\n",
" root_unity = np.exp(2j * np.pi / N)\n",
" root_unity_power = 1\n",
"\n",
" ax.grid()\n",
" ax.set_xlim([-1.4, 1.4])\n",
" ax.set_ylim([-1.4, 1.4])\n",
" ax.set_xlabel('$\\mathrm{Re}$')\n",
" ax.set_ylabel('$\\mathrm{Im}$')\n",
" ax.set_title('Roots of unity for $N=%d$' % N)\n",
"\n",
" for n in range(0, N):\n",
" colorPlot = 'r' if gcd(n, N) == 1 else 'k'\n",
" plot_vector(root_unity_power, color=colorPlot)\n",
" ax.text(np.real(1.2*root_unity_power), np.imag(1.2*root_unity_power),\n",
" r'$\\rho_{%s}^{%s}$' % (N, n), size='14',\n",
" color=colorPlot, ha='center', va='center')\n",
" root_unity_power *= root_unity\n",
"\n",
" circle_unit = plt.Circle((0, 0), 1, color='lightgray', fill=0)\n",
" ax.add_artist(circle_unit)\n",
"\n",
"plt.figure(figsize=(12, 4))\n",
"ax = plt.subplot(1, 3, 1)\n",
"plot_root_unity(N=8, ax=ax) \n",
"ax = plt.subplot(1, 3, 2)\n",
"plot_root_unity(N=11, ax=ax)\n",
"ax = plt.subplot(1, 3, 3)\n",
"plot_root_unity(N=12, ax=ax)\n",
"plt.tight_layout()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercises and Results"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import libpcp.exp\n",
"show_result = True"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 1: Approximation of Exponential Function via Power Series**

\n", "Implement a function

"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#\n", "Implement a function

`exp_power_series`

with input arguments $z\\in\\mathbb{C}$ and $N\\in\\mathbb{N}$, which outputs an approximation $\\sum_{n=0}^{N} \\frac{z^n}{n!}$ of $\\exp(z)$. Similarly, implement a function `exp_limit_compound`

that approximates $\\exp(z)$ via $\\left(1+\\frac{z}{N}\\right)^{N}$. Test the two functions for various input arguments and compare the results with the NumPy function `np.exp`

. In particular, compare the approximation quality of `exp_power_series`

and `exp_limit_compound`

for increasing $N$ (fixing a complex number $z$).\n",
"\n",
"**Exercise 2: Gaussian Function**

\n", "Gaussian functions are often used in statistics to represent the probability density function of a normally distributed random variable with expected value $\\mu$ and variance $\\sigma^2$. For these parameters, the Gaussian function $g:\\mathbb{R}\\to \\mathbb{R}$ is defined by\n", "\n", "$$\n", " g(x):= \\frac{1}{\\sigma\\sqrt{2\\pi}} \\exp\\left(-\\frac{1}{2} \\left(\\frac{x-\\mu}{\\sigma} \\right)^2 \\right).\n", "$$\n", "\n", "Using

"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#\n", "Gaussian functions are often used in statistics to represent the probability density function of a normally distributed random variable with expected value $\\mu$ and variance $\\sigma^2$. For these parameters, the Gaussian function $g:\\mathbb{R}\\to \\mathbb{R}$ is defined by\n", "\n", "$$\n", " g(x):= \\frac{1}{\\sigma\\sqrt{2\\pi}} \\exp\\left(-\\frac{1}{2} \\left(\\frac{x-\\mu}{\\sigma} \\right)^2 \\right).\n", "$$\n", "\n", "Using

`np.exp`

, implement a function `compute_gaussian_1D`

that inputs a NumPy array `X`

(as well as input arguments for $\\mu$ and $\\sigma$) and evaluates the Gaussian function `X`

in a point-wise fashion. Plot the result for an input vector and various choices of $\\mu$ and $\\sigma$.\n",
"\n",
"**Exercise 3: Spiral Generation**

\n", "Implement a function

"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#\n", "Implement a function

`generate_spiral`

that generates a spiral of increasing radius (Hint: Make use of the exponential function). The function should have the following arguments:\n",
"- \n",
"
`rad_start`

: Radius to start with \n",
" `rad_end`

: Radius to stop with \n",
" `num_rot`

: Number of rotations \n",
" `angle_start`

: Angle to start with (given in degrees) \n",
" `N`

: Number of data points to represent the spiral \n",
"