Taylor Series

Quarto
Python
desp
Author
Published

October 24, 2022

\[\begin{equation} f(x) = \sum_{n=0}^{\infty}\frac{f^{(n)}(a)}{n!}(x-a)^{n} \end{equation}\]

Let \(f(x)\) be a function that is \(n+1\) times differentiable on an interval \(I\) containing \(a\) and let \(P_n(x)\) be the \(n\)th degree Taylor polynomial for \(f(x)\) about \(a\). Then, there exists a number \(c\) between \(a\) and \(x\) such that: \[\begin{equation} f(x)=P_n(x)+R_n(x), \end{equation}\] where the remainder \(R_n(x)\) is given by: \[\begin{equation} R_n(x)=\frac{f^{(n+1)}(c)}{(n+1)!}(x-a)^{n+1}. \end{equation}\]

Define the function to be approximated

Code
import torch
import math
import matplotlib.pyplot as plt
import numpy as np
# Define the sine function to be approximated
def f(x):
    return torch.sin(x)

x = torch.linspace(-3.14, 3.14, 100)
y = f(x)
plt.plot(x, y)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Sine Function')
plt.show()

First order Taylor approximation for f(x) at x = 0

Code
x = torch.tensor([0.], requires_grad=True)
y = f(x)
approx = y + torch.autograd.grad(y, x, create_graph=True)[0] * x
x_vals = torch.linspace(-np.pi, np.pi, 100)
y_vals = f(x_vals)
approx_vals = (approx.detach() + torch.autograd.grad(approx, x, create_graph=True)[0] * x_vals).detach()
plt.plot(x_vals.numpy(), y_vals.numpy(), label='sin(x)')
plt.plot(x_vals.numpy(), approx_vals.numpy(), label='approx')
plt.legend()
plt.show()

Find the nth order Taylor approximation for f(x) at x = 0

Code
def fact(n):
    return math.factorial(n)

def nth_deriv(f, a, n):
    if isinstance(a, (float, int)):
        a = torch.tensor([a], dtype=torch.float, requires_grad=True)
    else:
        a = a.clone().detach().requires_grad_(True)
    
    y = f(a)
    for i in range(n):
        y = torch.autograd.grad(y, a, create_graph=True)[0]
    return y



# nth degree Taylor polynomial of f(x) around x=a
def taylor(f, x, n):
    result = torch.zeros_like(x)
    for i in range(n+1):
        result += nth_deriv(f, 0, i) / torch.tensor(math.factorial(i), dtype=torch.float32) * (x**i)
    return result
x_vals = torch.linspace(-math.pi, math.pi, 200)
plt.plot(x_vals.numpy(), f(x_vals).numpy(), label='f(x)', lw=5)
plt.plot(x_vals.numpy(), taylor(f, x_vals, 1).detach().numpy(), label='Taylor approximation, n=1')
plt.plot(x_vals.numpy(), taylor(f, x_vals, 3).detach().numpy(), label='Taylor approximation, n=3')
plt.plot(x_vals.numpy(), taylor(f, x_vals, 5).detach().numpy(), label='Taylor approximation, n=5')
plt.legend()
plt.show()

Plot of the function g(x) = x^2 and its Taylor approximations up to degree 1 and degree 2 centered at x=0.

Code
x_vals = torch.linspace(-4, 4, 100)

def g(x):
    return x**2

def taylor(f, x, n):
    x = x.unsqueeze(-1)
    y = f(x)
    for i in range(1, n+1):
        y += (x - x[0])**i / torch.tensor([math.factorial(i)]).float() * f(x[0] + 0.0)
    return y


plt.plot(x_vals.numpy(), g(x_vals).numpy(), label='g(x)', lw=5)
plt.plot(x_vals.numpy(), taylor(g, x_vals, 1).detach().numpy(), label='Taylor approximation, n=1')
plt.plot(x_vals.numpy(), taylor(g, x_vals, 3).detach().numpy(), label='Taylor approximation, n=3')
plt.plot(x_vals.numpy(), taylor(g, x_vals, 5).detach().numpy(), label='Taylor approximation, n=5')
plt.legend()
plt.show()