In [1]:
%matplotlib inline
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
In [2]:
# You can change the default figure size to be a bit larger if you want,
# uncomment the next line for that:
plt.rc('figure', figsize=(10, 6))
The function plot_taylor_approximations included here was written by Fernando Perez and was part of work on the original IPython project. Although attribution seems to have been lost over time, we gratefully acknowledge FP and thank him for this code!
In [3]:
def plot_taylor_approximations(func, x0=None, orders=(2, 4), xrange=(0,1), yrange=None, npts=200):
"""Plot the Taylor series approximations to a function at various orders.
Parameters
----------
func : a sympy function
x0 : float
Origin of the Taylor series expansion. If not given, x0=xrange[0].
orders : list
List of integers with the orders of Taylor series to show. Default is (2, 4).
xrange : 2-tuple or array.
Either an (xmin, xmax) tuple indicating the x range for the plot (default is (0, 1)),
or the actual array of values to use.
yrange : 2-tuple
(ymin, ymax) tuple indicating the y range for the plot. If not given,
the full range of values will be automatically used.
npts : int
Number of points to sample the x range with. Default is 200.
"""
if not callable(func):
raise ValueError('func must be callable')
if isinstance(xrange, (list, tuple)):
x = np.linspace(float(xrange[0]), float(xrange[1]), npts)
else:
x = xrange
if x0 is None: x0 = x[0]
xs = sp.Symbol('x')
# Make a numpy-callable form of the original function for plotting
fx = func(xs)
f = sp.lambdify(xs, fx, modules=['numpy'])
# We could use latex(fx) instead of str(), but matploblib gets confused
# with some of the (valid) latex constructs sympy emits. So we play it safe.
plt.plot(x, f(x), label=str(fx), lw=2)
# Build the Taylor approximations, plotting as we go
apps = {}
for order in orders:
app = fx.series(xs, x0, n=order).removeO()
apps[order] = app
# Must be careful here: if the approximation is a constant, we can't
# blindly use lambdify as it won't do the right thing. In that case,
# evaluate the number as a float and fill the y array with that value.
if isinstance(app, sp.numbers.Number):
y = np.zeros_like(x)
y.fill(app.evalf())
else:
fa = sp.lambdify(xs, app, modules=['numpy'])
y = fa(x)
tex = sp.latex(app).replace('$', '')
plt.plot(x, y, label=r'$n=%s:\, %s$' % (order, tex) )
# Plot refinements
if yrange is not None:
plt.ylim(*yrange)
plt.grid()
plt.legend(loc='best').get_frame().set_alpha(0.8)
# For an expression made from elementary functions, we must first make it into
# a callable function, the simplest way is to use the Python lambda construct.
# plot_taylor_approximations(lambda x: 1/sp.cos(x), 0, [2,4,6,8], (0, 2*sp.pi), (-5,5))
In [4]:
plot_taylor_approximations(sp.sin, 0, [2, 4, 6, 8], (0, 2*sp.pi), (-2,2))
It is common in physics and engineering to represent transcendental functions and other nonlinear expressions using a few terms from a Taylor series. This series provides a fast and efficient way to compute quantities such as $\mathrm{sin}(x)$ or $e^x$ to a prescribed error. Learning how to calculate the series representation of these functions will provide practical experience with the Taylor series and help the student understand the results of Python methods designed to accelerate and simplify computations. The series can be written generally as:
$$f(x) = f{\left (0 \right )} + x \left. \frac{d}{d x} f{\left (x \right )} \right|_{\substack{ x=0 }} + \frac{x^{2}}{2} \left. \frac{d^{2}}{d x^{2}} f{\left (x \right )} \right|_{\substack{ x=0 }} + \frac{x^{3}}{6} \left. \frac{d^{3}}{d x^{3}} f{\left (x \right )} \right|_{\substack{ x=0 }} + \frac{x^{4}}{24} \left. \frac{d^{4}}{d x^{4}} f{\left (x \right )} \right|_{\substack{ x=0 }} + \frac{x^{5}}{120} \left. \frac{d^{5}}{d x^{5}} f{\left (x \right )} \right|_{\substack{ x=0 }} + \mathcal{O}\left(x^{6}\right)$$Of equal importance the Taylor series permits discrete representation of derivatives and is a common way to perform numerical integration of partial and ordinary differential equations. Expansion of a general function $f(x)$ about a point, coupled with algebraic manipulation, will produce expressions that can be used to approximate derivative quantities. Although any order of derivative can be computed, this lesson will focus on first and second derivatives that will be encountered in the diffusion equation.
You will practice the following skills:
Optional challenge: A list is one of the fundamental data structures within Python. Numpy (a Python library) and other parts of Python libraries use vectorized computations. From Wikipedia, vectorization is "a style of computer programming where operations are applied to whole arrays instead of individual elements."
With this in mind, we certainly can iterate over our list of points and apply the function that you will soon write in an element by element fashion, however, it is a more common practice in Python and other modern languages to write vectorized code. If this is your first exposure to vectorized computation, I recommend two initial strategies: write out your algorithms and use "classic" flow control and iteration to compute the results. From that point you will more easily see the strategy you should use to write vectorized code.
Using the discrete forms of the first and second derivatives (based on central differences) can you devise a vectorized operation that computes the derivative without looping in Python?
A high quality communication provides an organized, logically progressing, blend of narrative, equations, and code that teaches the reader a particular topic or idea. You will be assessed on:
If your notebook is just computer code your assignment will be marked incomplete.
Ideas relating to sequences, series, and power series are used in the formulation of integral calculus and in the construction of polynomial representations of functions. The limit of functions will also be investigated as boundary conditions for differential equations. For this reason understanding concepts related to sequences and series are important to review.
A sequence is an ordered list of numbers. A list such as the following represents a sequence:
$$a_1, a_2, a_3, a_4, \dots, a_n, \dots $$The sequence maps one value $a_n$ for every integer $n$. It is typical to provide a formula for construction of the nth term in the sequence. While ad-hoc strategies could be used to develop sequences using SymPy and lists in Python, SymPy has a sequence class that can be used. A short demonstration is provided next:
In [5]:
import sympy as sp
x, y, z, t, a = sp.symbols('x y z t a')
k, m, n = sp.symbols('k m n', integer=True)
f, g, h = sp.symbols('f g h', cls=sp.Function)
sp.var('a1:6')
sp.init_printing()
It is important to read about SymPy symbols at this time.
We can generate a sequence using SeqFormula.
In [6]:
a1 = sp.SeqFormula(n**2, (n,0,5))
list(a1)
Out[6]:
if we want the limit of the sequence at infinity:
$$[0, 1, 4, 9, \ldots ]$$we can use limit_seq:
In [7]:
sp.limit_seq(a1.formula, n)
Out[7]:
In [8]:
# Your code here.
In [9]:
a2 = sp.Sum(1/2**n, (n,0,1))
a2
Out[9]:
In [10]:
a2.doit()
Out[10]:
In [11]:
a4 = sp.Sum(n**2, (n,0,5))
a4
Out[11]:
In [12]:
a5 = sp.Sum(k**2, (k, 1, m))
a5
Out[12]:
In [13]:
a4.doit()
Out[13]:
In [14]:
a5.doit()
Out[14]:
A power series is of the form:
$$ \sum_{n=0}^{\infty} M_{n} x^{n} = M_0 + M_1 x + M_2 x^2 + \cdots $$
In [15]:
M = sp.IndexedBase('M')
In [16]:
sp.Sum(M[n]*x**n, (n,0,m))
Out[16]:
We can define the series about the point $a$ as follows:
In [17]:
sp.Sum(M[n]*(x-a)**n, (n,0,m))
Out[17]:
SymPy has a function that can take SymPy expressions and represent them as power series:
In [18]:
sp.series(f(x), x, x0=0)
Out[18]:
In [19]:
# Your code here.
Below we present a derivation of Taylor's series and small algebraic argument for series representations of functions. In contrast to the ability to use sympy functions without any deeper understanding, these presentations are intended to give you insight into the origin of the series representation and the factors present within each term. While the algebraic presentation isn't a general case, the essential elements of a general polynomial representation are visible.
The function $f(x)$ can be expanded into an infinite series or a finite series plus an error term. Assume that the function has a continuous nth derivative over the interval $a \le x \le b$. Integrate the nth derivative n times:
$$\int_a^x f^n(x) dx = f^{(n-1)}(x) - f^{(n-1)}(a)$$The power on the function $f$ in the equation above indicates the order of the derivative. Do this n times and then solve for f(x) to recover Taylor's series. One of the key features in this derivation is that the integral is definite. This derivation is outlined on Wolfram’s Mathworld.
As a second exercise, assume that we wish to expand sin(x) about x=0. First, assume that the series exists and can be written as a power series with unknown coefficients. As a first step, differentiate the series and the function we are expanding. Next, let the value of x go to the value of the expansion point and it will be possible to evaluate the coefficients in turn:
We can choose an expansion point (e.g. $x = 0$) and differentiate to get a set of simultaneous equations permitting determination of the coefficients. The computer algebra system can help us with this activity:
In [20]:
import sympy as sp
sp.init_printing()
In [21]:
x, A, B, C, D, E = sp.symbols('x, A, B, C, D, E')
To help us get our work done we can use sympy's diff function. Testing this function with a known result, we can write:
In [22]:
sp.diff(sp.sin(x),x)
Out[22]:
A list comprehension is used to organize the results. In each iteration the exact function and the power series are differentiated and stored as an element of a list. The list can be inspected and a set of simultaneous equations can be written down and solved to determine the values of the coefficients. Casting the list as a sympy Matrix object clarifies the correspondance between entries in the list.
In [23]:
orderOfDifferentiation = 1
powerSeries = A+B*x+C*x**2+D*x**3+E*x**4
# Differentiate, element by element, the list [sp.sin(x),powerSeries]
[sp.diff(a,x,orderOfDifferentiation) for a in [sp.sin(x),powerSeries]]
Out[23]:
A list comprehension can be used to organize and extend the results further. We can wrap the list above into another list that changes the order of differentiation each iteration.
In [24]:
maximumOrder = 5
funcAndSeries = [[sp.diff(a,x,order) for a in [sp.sin(x),powerSeries]] for order in range(maximumOrder)]
funcAndSeries
Out[24]:
Casting the results as a sympy Matrix object the list is more easily viewed in the Jupyter notebook:
In [25]:
sp.Matrix(funcAndSeries)
Out[25]:
In [26]:
# Your code here if you feel you need it.
Your markdown here.
In [27]:
from sympy import init_printing, symbols, Function
init_printing()
x, c = symbols("x,c")
f = Function("f")
f(x).series(x, x0=c, n=3)
Out[27]:
One of the major uses of Taylor's series in computation is for the evaluation of derivatives. Take note of the fact that the derivatives of a function appear in the evaluation of the series.
It may be straightforward to compute the derivative of some functions. For example:
$$f(x) = x^2$$$$f'(x) = 2x$$In numerical computing situations there is no analytical solution to the problem being solved and therefore no function to integrate or differentiate. The approximate solution is available as a list of discrete points in the domain of the problem's independent variables (e.g. space, time). The values could be represented as a list of numbers:
$$\{f(x_0), f(x_1), f(x_2), ...\}$$The neighboring points $f(x_0)$ and $f(x_1)$ are seperated by a distance $\Delta x = x_1 - x_0$ in the independent variable. Although this will not be apparent from the values, it is implicit in the structure of the data. Taylor's series can be used to compute approximate derivatives of the unknown function directly from the list of points in this situation. We are going to compute a series expansion for an unknown function $f(x)$ in the vicinity of the point $c$ and then examine the relationship between that function and it's derivative quantities at a point $c\pm h$. The goal of the activity is to see if we can find expressions for the derivatives using the data point of interest ($c$) and its neighbors ($c \pm h$). We are going to use the idea of forward and backward differences.
In the figure below we indicate the following:
Imagine that we take the above series expansion and use it to compute the value of the function near the point $c$. Let us evaluate this series by adding and subtracting to the independent varable the quantity $h$. To accomplish this we write down the series expansion for our function about the point $c$, then we let the independent variable $x \rightarrow c+h$ and $c-h$.
In [28]:
x, h, c = sp.symbols("x,h,c")
f = sp.Function("f")
# the .subs() method replaces occurences of 'x' with something else
taylorExpansionPlus = f(x).series(x, x0=c, n=3).removeO().subs(x,c+h)
taylorExpansionMinus = f(x).series(x, x0=c, n=3).removeO().subs(x,c-h)
In [29]:
taylorExpansionPlus
Out[29]:
Meaning that:
$$ f(c+h) = \frac{h^{2}}{2} \left. \frac{d^{2}}{d \xi_{1}^{2}} f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=c }} + h \left. \frac{d}{d \xi_{1}} f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=c }} + f{\left (c \right )} $$
In [30]:
taylorExpansionMinus
Out[30]:
Meaning that:
$$ f(c-h) = \frac{h^{2}}{2} \left. \frac{d^{2}}{d \xi_{1}^{2}} f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=c }} - h \left. \frac{d}{d \xi_{1}} f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=c }} + f{\left (c \right )} $$Inspection of the results shows that the signs on the terms containing the first derivative are different between the two expressions. We can use this to our advantage in solving for the derivative terms explicitly. Note that each grouped expression is equal to zero as is the default in sympy.
Find the first derivative in this expression:
In [31]:
(taylorExpansionMinus-f(c-h))-(taylorExpansionPlus-f(c+h))
Out[31]:
Remember that sympy expressions are zero by default. So this is true:
Find the second derivative in this expression:
In [32]:
(taylorExpansionMinus-f(c-h))+(taylorExpansionPlus-f(c+h))
Out[32]:
Similarly:
$$ h^{2} \left. \frac{d^{2}}{d \xi_{1}^{2}} f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=c }} + 2 f{\left (c \right )} - f{\left (c - h \right )} - f{\left (c + h \right )} = 0 $$Solve each of the equations above for the derivative quantites and you will have recovered the formula for discrete derivatives.
Use these expressions to complete your assignment as outlined above.
In [ ]: