Taylor Series Approximation

In some cases, in engineering or real world technical problems, we are not interested to find the exact solution of a problem. If we have a good enough approximation, we can consider that we’ve found the solution of the problem.

For example, if we want to compute the trigonometric function f(x)=sin(x) with a hand held calculator, we have two options:

• use the actual trigonometric function sin(x), if the calculator has the function embedded (available if it’s a scientific calculator)
• use a polynomial as an approximation of the sin(x) function and compute the result with any calculator, or even by hand

In general, any mathematical function f(x), with some constraints, can be approximated by a polynomial P(x):

\[P(x)=a_{0}+a_{1} \cdot x+a_{2} \cdot x^{2} + … + a_{n} \cdot x^{n}\]

Weierstrass approximation theorem

First, let’s put down what the theorem sounds like. After, we’ll try to explain it a bit.

Theorem: For a given function f(x), which is defined and continuous on the interval [a, b], there is always a polynomial P(x), also defined on the interval [a, b], with the property:

\[\left | f(x)-P(x) \right | < \varepsilon\]

for any x ∈ [a, b]. and a given ε > 0.

The theorem says that for any function f(x), which is continuous and defined between the points a and b, there is always a polynomial P(x), which can approximate the function f(x) with a small error ε, in the same interval [a, b].

Taylor’s polynomials

Example: Let’s approximate the function f(x)=sin(x) with a polynomial of order 3, around the point x₀ = 0. Using the determined polynomial, approximate the value of sin(0.1).

Explanation: “around the point x₀” means that f⁽ⁿ⁾(x₀) = P⁽ⁿ⁾(x₀), which means that the evaluation of the function and its derivatives in the point x₀ is equal to the evaluation of the polynomial and its derivatives.

Step 1. Write the polynomial of order 3.

\[P(x)=a_{0}+a_{1} \cdot x+a_{2} \cdot x^{2} +a_{3} \cdot x^{3}\]

Step 2. Calculate the 3^rd order derivatives of P(x). We need them in order to find out the values of the coefficients a₀, a₁, a₂ and a₃.

\[ \begin{split}
P(x)&=a_{0}+a_{1}x+a_{2}x^{2}+a_{3}x^{3}\\
P^{\prime}(x)&=a_{1}+2a_{2}x+3a_{3}x^{2}\\
P^{\prime\prime}(x)&=2a_{2}+6a_{3}x\\
P^{\prime\prime\prime}(x)&=6a_{3}
\end{split} \]

Step 3. Calculate P⁽ⁿ⁾(x₀).

\[ \begin{split}
P(0)&=a_{0}\\
P^{\prime}(0)&=a_{1}\\
P^{\prime\prime}(0)&=2a_{2}\\
P^{\prime\prime\prime}(0)&=6a_{3}
\end{split} \]

Step 4. Calculate the 3^rd order derivatives of f(x).

\[ \begin{split}
f(x)&=\sin(x)\\
f^{\prime}(x)&=\cos(x)\\
f^{\prime\prime}(x)&=-\sin(x)\\
f^{\prime\prime\prime}(x)&=-\cos(x)
\end{split} \]

Step 5. Calculate f⁽ⁿ⁾(x₀).

\[ \begin{split}
f(0)&=\sin(0)&=0\\
f^{\prime}(0)&=\cos(0)&=1\\
f^{\prime\prime}(0)&=-\sin(0)&=0\\
f^{\prime\prime\prime}(0)&=-\cos(0)&=-1
\end{split} \]

Step 6. Calculate the coefficients a₀, a₁, a₂ and a₃.

\[ \begin{split}
P(0)&=f(0) &\Rightarrow a_{0}&=0\\
P^{\prime}(0)&=f^{\prime}(0) &\Rightarrow a_{1}&=1\\
P^{\prime\prime}(0)&=f^{\prime\prime}(0) &\Rightarrow a_{2}&=0\\
P^{\prime\prime\prime}(0)&=f^{\prime\prime\prime}(0) &\Rightarrow a_{3}&=-\frac{1}{6}
\end{split} \]

Step 7. Write the final form of the approximation polynomial

\[P(x)=x-\frac{1}{6}x^{3}\]

Step 8. Approximate the value of sin(0.1) using the polynomial.

\[f(0.1)=\sin(0.1)\cong P(0.1)=0.1-\frac{1}{6} \cdot 0.1^{3}\approx 0.0998333\]

Using Scilab we can compute sin(0.1) just to compare with the approximation result:

--> sin(0.1)

 ans  =
   0.0998334

As you can see, the approximation with the polynomial P(x) is quite accurate, the result being equal up to the 7^th decimal.

Accuracy of the approximation

Since the approximation is done around the point x₀, evaluating the polynomial in x₀ gives pretty good results. But what happens if we move away from x₀, how accurate the approximation will still be?

To evaluate the polynomial for several points around x0, we are going to use a Scilab script.

deff('y=P(x)','y=x-(1/6)*x^3')
x=-%pi/2:%pi/10:%pi/2;
e_x=abs(1-P(x)./sin(x))*100;
subplot(2,1,1)
plot(x,P(x),'sb-')
plot(x,sin(x),'sr-')
legend('$\Large{P(x)}$','$\Large{\sin(x)}$',2)
ylabel('P(x), sin(x)','FontSize',2)
title('x-engineer.org','FontSize',2)
xgrid()
subplot(2,1,2)
plot(x,e_x,'sb-')
xgrid()
xlabel('x [rad]','FontSize',2)
ylabel('Relative error [%]','FontSize',2)

What the Scilab script does:

• defines the polynomial P(x) as a function
• defines the range of x points from -π/2 up to π/2, using π/10 increments
• calculates the relative error e(x) between sin(x) and P(x)
• plots the values of sin(x) and P(x)
• plots the relative error for each x point

The relative error e(x), in percentage, is calculated as:

\[e(x) = \left | 1-\frac{P(x)}{\sin(x)} \right | \cdot 100\]

Running the Scilab script will plot the following:

As you can see, the polynomial approximation of sin(x) is quite accurate around x₀ but starts to diverge as we move toward the margins of the interval. For the extreme points, -π/2 and π/2, the relative error is around 7.5%.

In order to increase the accuracy of the polynomial approximation we need to increase the order of the polynomial. The higher the degree of the polynomial the better the approximation for a wider interval.

Taylor’s theorem

The polynomial P(x) used in the example above is a specific case of a Taylor series for function approximation.

Theorem: Any function f(x) can be written as:

\[f(x)=P_{n}(x)+R_{n}(x)\]

with P(x) being Taylor’s polynomial and R(x) being Taylor’s remainder:

\[ \begin{split}
P_{n}(x)&=f(x_{0})+f^{\prime}(x_{0})(x-x_{0})+\frac{f^{\prime\prime}(x_{0})}{2!}(x-x_{0})^{2}+…+\frac{f^{(n)}(x_{0})}{n!}(x-x_{0})^{n}\\
R_{n}(x)&=\left | \frac{f^{(n+1)}(\xi(x)) }{(n+1)!}(x-x_{0})^{n+1}\right |
\end{split} \]

P(x) can be written is a short form, as:

\[P_{n}(x)=\sum_{k=0}^{n}\frac{f^{(k)}(x_{0})}{k!}(x-x_{0})^{k}\]

where:

• f ∈ Cⁿ([a, b]) – which means that f(x) is continuous and derivable on an interval [a, b]
• f⁽ⁿ⁺¹⁾ ∈ [a, b] – which means that all n+1 derivatives of f(x) exist in [a, b]
• x₀ ∈ [a, b] – which means that x₀ is a point between a and b
• ξ(x) ∈ (a, b) – which means that there is a point ξ between a and b which is greater than a and smaller than b

Example. For the function f(x), defined as:

\[f(x)=\sqrt{1+x}=(1+x)^{\frac{1}{2}}\]

1. find the Taylor polynomial of order 3 around x₀ = 0.
2. approximate the value of √1.1
3. find the maximum error for the approximation

Step 1. Calculate the 3^rd order derivatives of f(x).

\[ \begin{split}
f(x)&=(1+x)^{\frac{1}{2}}\\
f^{\prime}(x)&=\frac{1}{2}(1+x)^{-\frac{1}{2}}\\
f^{\prime\prime}(x)&=-\frac{1}{4}(1+x)^{-\frac{3}{2}}\\
f^{\prime\prime\prime}(x)&=\frac{3}{8}(1+x)^{-\frac{5}{2}}
\end{split} \]

Step 2. Calculate f⁽ⁿ⁾(x₀).

\[ \begin{split}
f(0)&=1\\
f^{\prime}(0)&=\frac{1}{2}\\
f^{\prime\prime}(0)&=-\frac{1}{4}\\
f^{\prime\prime\prime}(0)&=\frac{3}{8}
\end{split} \]

Step 3 (solves a). Calculate Taylor’s polynomial P₃(x)

\[ \begin{split}
P_{3}(x)&=f(0)+f^{\prime}(0)x+\frac{f^{\prime\prime}(0)}{2!}x^{2}+\frac{f^{\prime\prime\prime}(0)}{3!}x^{3}\\
&=1+\frac{1}{2}x-\frac{1}{4}\cdot\frac{1}{2!}x^{2}+\frac{3}{8}\cdot\frac{1}{3!}x^{3}\\
&=1+\frac{1}{2}x-\frac{1}{8}x^{2}+\frac{1}{16}x^{3}
\end{split} \]

Step 4 (solves b). Approximate the value of √1.1

1.1 = 1 + 0.1, which means that x=0.1.

\[P_{3}(0.1)=1+\frac{1}{2}(0.1)-\frac{1}{8}(0.1)^{2}+\frac{1}{16}(0.1)^{3} \approx 1.0488125\]

To verify our result let’s use Scilab function sqrt():

--> sqrt(1+0.1)

 ans  =
   1.0488088

Step 5. Calculate the 4^th derivative of f(x).

\[f^{\prime\prime\prime\prime}(x)=-\frac{15}{16}(1+x)^{-\frac{7}{2}}\]

Step 6. Calculate the 4^th derivative of f(ξ).

\[f^{\prime\prime\prime\prime}(\xi)=-\frac{15}{16}(1+\xi)^{-\frac{7}{2}}\]

where ξ ∈ [0, x].

Step 7 (solves c). Calculate R₃(x) for a given ξ.

ξ must be chosen in such a way that the 4^th derivative of f(ξ) has maximum value. In our case this makes ξ = 0.

\[ \begin{split}
R_{3}(0.1)&=\left | \frac{-\frac{15}{16}(1+\xi)^{-\frac{7}{2}}}{4!}(0.1)^{4} \right | \\
&= \frac{15}{16 \cdot 24}(0.1)^{4} \cdot \underset{\xi\in [0, 0.1]}{\max} (1+\xi)^{-\frac{7}{2}}\\
&= \frac{15}{16 \cdot 24}(0.1)^{4} \cdot 1\\
&= 0.00000391
\end{split} \]

This means that the error of the polynomial approximation will not exceed 3.91·10^-6.

Examples with Scilab programming

Let’s use Scilab to calculate the Taylor series approximations for a couple of functions. To visualise the impact of the order of the approximation polynomial, we’ll use Scilab plot() function. For the functions f(x) and P(x) given below, we’ll plot the exact solution and Taylor approximation using a Scilab script.

Example 1. Exponential function f(x)=e^x.

\[ \begin{split}
f(x)&=e^{x}\\
P(x)&=1+x+\frac{x^{2}}{2!}+\frac{x^{3}}{3!}+…\\
&=\sum_{n=0}^{\infty}\frac{x^{n}}{n!}
\end{split} \]

The Scilab script will define a custom Scilab function for P(x).

// Define approximation polynomial
function [y] = P(x,n)
    y=zeros(1,length(x));
    for i=1:n+1
        y=y+x.^(i-1)./factorial(i-1);
    end
endfunction

// Initialize figure
figure(0)
clf()
ha=gca();
hf=gcf();
hf.background = color('white');
color_list=['red' 'orange' 'magenta' 'darkmagenta' 'purple' 'green' 'darkcyan' 'blue'];

// Calculate and plot approximations P(x)
x=-3:.01:2;
n=8;
for i= 1:n
    plot(x,P(x,i-1));
end
for i=1:n
    ha.children(n-i+1).children.foreground = color(color_list(i));
end

// Plot f(x)
plot(x,exp(x),'k','LineWidth',2)
xgrid()
xlabel('x','FontSize',2)
ylabel('$\Large{P_{n}(x),\quad e^{x}}$')
legend(['n=0' 'n=1' 'n=2' 'n=3' 'n=4' 'n=5' 'n=6' 'n=7' '$\Large{e^{x}}$'],2)
title('x-engineer.org','FontSize',2)

The output of the Scilab script is plotted below, together with an animation of the same approximation.

Example 2. Natural logarithm function f(x)=ln(1+x).

\[ \begin{split}
f(x)&=\ln(1+x)\\
P(x)&=x-\frac{x^{2}}{2}+\frac{x^{3}}{3}-…\\
&=\sum_{n=1}^{\infty}(-1)^{n+1}\frac{x^{n}}{n}
\end{split} \]

The Scilab script will define a custom Scilab function for P(x).

// Define approximation polynomial
function [y] = P(x,n)
    y=zeros(1,length(x));
    for i=1:n
        y=y+(-1)^(i+1).*(x.^i)/(i);
    end
endfunction

clc()
// Initialize figure
figure(0)
clf()
ha=gca();
hf=gcf();
hf.background = color('white');
color_list=['red' 'orange' 'magenta' 'darkmagenta' 'purple' 'green' 'darkcyan' 'blue'];

// Calculate and plot approximations P(x)
x=-1:.01:1;
n=8;
for i= 1:n
    plot(x,P(x,i-1));
end
for i=1:n
    ha.children(n-i+1).children.foreground = color(color_list(i));
end

// Plot f(x)
plot(x,log(1+x),'k','LineWidth',2)
xgrid()
xlabel('x','FontSize',2)
ylabel('$\Large{P_{n}(x),\quad \ln(1+x)}$')
legend(['n=0' 'n=1' 'n=2' 'n=3' 'n=4' 'n=5' 'n=6' 'n=7' '$\Large{\ln(1+x)}$'],4)
title('x-engineer.org','FontSize',2)

The output of the Scilab script is plotted below.

In order to try out other functions and their Taylor series approximation, redefine the P(x,n) function from the Scilab script.