Project 3: Numerical Derivatives and Taylor Series Approximations

CEGM1000 MUDE: Week 1.5, Friday, Oct 4, 2024.

Overview:

Numerical derivatives are required to solve differential equations in numerical modelling and they can also be applied to data. To understand its intricacies, the ice thickness data of the Nenana River will be applied using the three derivative approximations.

Taylor Series Expansions are necessary to understand the limitations of numerical models based on finite differences/finite elements. In particular, the approximation of derivatives/integrals and its effect on representing physical systems. In the majority of cases, we use Taylor as an approximation of non-linear functions that do not have mathematically simple solutions and instead we use a polynomial to locally approximate. However, it is crucial to understand that the result is only accurate over a small interval. For this project, we will look at a simple function that we can solve analytically and compare our Taylor approximations.

Note on programming implementation: one goal of this assignment is for you to be aware of how the numerical implementation differs from the analytic expressions, as well as how the implementation differs between numerical schemes. Typically this is straightforward for the calculations in general, but special consideration must be made when evaluating equations at the edges of the data set (or function domain). In 1D problems, this means the first and last points.

Note also that we use Pandas in this notebook to make the data handling easier. You are not expected to be able to use it (yet), but now that you are an expert with Python dictionaries, you will probably be able to recognize what it is doing easily.

import numpy as np
import matplotlib.pylab as plt
import pandas as pd

Part 0: Data Import and Exploration

In the following, the ice thickness [cm] data of the Nenana River is plotted over time, and we are interested in the derivative of ice thickness.

Task 0.1: run the cell below to create a plot, then reflect on this question: what does the time derivative represent?

data=pd.read_csv(filepath_or_buffer='justIce.csv',index_col=0)
data.index = pd.to_datetime(data.index, format="%Y-%m-%d")

plt.figure(figsize=(15,4))
plt.scatter(data.index,data, color='green', marker='x')
plt.xlabel('Year')
plt.ylabel('Ice thickness [cm]')
plt.grid()

Task 0.2: Before computing the numerical derivatives, reflect on why it would be a bad idea to do it over the entire data set? In other words, why is it a bad idea to compute the derivatives over the entire nearly 40-year period?

Task 0.3: Run the cell below to visualize what the data looks like for one year. Reflect on whether this is a suitable time period over which to compute derivatives.

data_2021 = data.loc['2021']

plt.figure(figsize=(15,4))
plt.scatter(data_2021.index,data_2021, color='green', marker='x')
plt.xlabel('Date')
plt.ylabel('Ice thickness [cm]')
plt.grid()

Part 1: Numerical Derivatives

For numerics we will use numpy due to its optimization for numerical computations. In the cell below, the conversion to numpy is already there as well as a transformation from date to time_days.

Task 1.1: Use Forward Difference to estimate the growth rate of the ice, using only the variables h_ice and t_days (i.e., no need to write a function).

Remember that one of our goals is to compare the numerical implementation, therefore, in the Tasks below you are asked to keep the function evaluations identical in each case, while only changing the evaluation points (i.e., the indices of the time array).

h_ice = (data_2021.to_numpy()).ravel()
t_days = ((data_2021.index - data_2021.index[0]).days).to_numpy()

# dh_dt_FD = YOUR_CODE_HERE

# SOLUTION:
dh_dt_FD = (h_ice[1:]-h_ice[:-1])/(t_days[1:]-t_days[:-1]) 

Great, you did it! However, computing the derivative is not enough. It is necessary to specify where that derivative was computed.

Task 1.2: Plot the growth rate together with the ice thickness, specifying the proper time where the derivative is computed.

Note that a second axis has been added to enable us to include the original ice thickness measurements for comparison.

fig, ax1 = plt.subplots(figsize=(15,4))

# ax1.scatter(YOUR_CODE_HERE, dh_dt_FE,
#             color='blue', marker='o', label='dh_dt_FE Forward Euler')

# SOLUTION:
ax1.scatter(t_days[:-1], dh_dt_FD,
            color='blue', marker='o', label='dh_dt_FD Forward Difference')

ax1.set_xlabel('Days')
ax1.set_ylabel('Growth Rate [cm/day]', color='blue')
ax1.tick_params(axis='y', labelcolor='blue')
ax1.grid()


ax2 = ax1.twinx()
ax2.scatter(t_days, h_ice,
            color='green', marker='x', label='h_ice Measurements')
ax2.set_ylabel('Ice thickness [cm]', color='green')
ax2.tick_params(axis='y', labelcolor='green')

handles, labels = ax1.get_legend_handles_labels()
handles2, labels2 = ax2.get_legend_handles_labels()
ax1.legend(handles + handles2, labels + labels2, loc='upper left')

plt.show()

Task 1.3: Use backward difference to compute the ice growth rate and plot it, together with the previous estimation and the ice thickness.

# dh_dt_BE = YOUR_CODE_HERE

# SOLUTION:
dh_dt_BD = (h_ice[1:]-h_ice[:-1])/(t_days[1:]-t_days[:-1]) 

fig, ax1 = plt.subplots(figsize=(15,4))


# ax1.scatter(YOUR_CODE_HERE, dh_dt_FD,
#             color='blue', marker='o', label='dh_dt_FE Forward Difference')
# ax1.scatter(YOUR_CODE_HERE, dh_dt_BD,
            # color='red', marker='o', label='dh_dt_BE Backward Difference')

# SOLUTION:
ax1.scatter(t_days[:-1], dh_dt_FD,
            color='blue', marker='o', label='dh_dt_FD Forward Difference')
ax1.scatter(t_days[1:], dh_dt_BD,
            color='red', marker='o', label='dh_dt_BD Backward Difference')

ax1.set_xlabel('Days')
ax1.set_ylabel('Growth Rate [cm/day]', color='blue')
ax1.tick_params(axis='y', labelcolor='blue')
ax1.grid()

ax2 = ax1.twinx()
ax2.scatter(t_days, h_ice,
            color='green', marker='x', label='h_ice Measurements')
ax2.set_ylabel('Ice thickness [cm]', color='green')
ax2.tick_params(axis='y', labelcolor='green')

handles, labels = ax1.get_legend_handles_labels()
handles2, labels2 = ax2.get_legend_handles_labels()
ax1.legend(handles + handles2, labels + labels2, loc='upper left')

plt.show()

Task 1.4: Now apply central differences and plot them.

Beware: if you follow the derivation in the textbook, the distances $x_{i+1}-x_i$ and $x_{i}-x_{i-1}$ are assumed to be the same for all points (and between $x_{i-1}$, $x_i$ and $x_{i+1}$). That is not the case in this data set. Therefore, you must evaluate the derivative at some location between each data point that satisfies this criteria. (Hint: the middle!). Also, remember, we are only changing the time indices.

Note that it is not necessary to define the "new" points that the derivative is evaluated at, you can write them directly in the scatter plot arguments.

# dh_dt_CD = YOUR_CODE_HERE

# SOLUTION:
dh_dt_CD = (h_ice[1:]-h_ice[:-1])/(t_days[1:]-t_days[:-1]) 

fig, ax1 = plt.subplots(figsize=(15,4))

# ax1.scatter(YOUR_CODE_HERE, dh_dt_FE,
#             color='blue', marker='o', label='dh_dt_FE Forward Euler')
# ax1.scatter(YOUR_CODE_HERE, dh_dt_BE,
#             color='red', marker='o', label='dh_dt_BE Backward Euler')
# ax1.scatter(YOUR_CODE_HERE, dh_dt_CD,
#             color='purple', marker='o', label='dh_dt_CD Central Difference')

# SOLUTION:
ax1.scatter(t_days[:-1], dh_dt_FD,
            color='blue', marker='o', label='dh_dt_FE Forward Difference')
ax1.scatter(t_days[1:], dh_dt_BD,
            color='red', marker='o', label='dh_dt_BE Backward Difference')
ax1.scatter((t_days[1:]+t_days[:-1])/2, dh_dt_CD,
            color='purple', marker='o', label='dh_dt_CD Central Difference')
ax1.set_xlabel('Days')
ax1.set_ylabel('Growth Rate [cm/day]', color='blue')
ax1.tick_params(axis='y', labelcolor='blue')
ax1.grid()

ax2 = ax1.twinx()
ax2.scatter(t_days, h_ice, color='green', marker='x', label='h_ice Measurements')
ax2.set_ylabel('Ice thickness [cm]', color='green')
ax2.tick_params(axis='y', labelcolor='green')

handles, labels = ax1.get_legend_handles_labels()
handles2, labels2 = ax2.get_legend_handles_labels()
ax1.legend(handles + handles2, labels + labels2, loc='upper left')

plt.show()

Part 2: Taylor Series

Definition

Recall that the Taylor series for one variable as previously described in the fundamentals section of the MUDE textbook as

$$ f(x)\approx f(x_0)+(x-x_0)\frac{\partial f(x_0)}{\partial x}+\frac{(x-x_0)^2}{2!}\frac{\partial^2 f(x_0)}{\partial x^2}+...+\frac{(x-x_0)^n}{n!}\frac{\partial^n f(x_0)}{\partial x^n} $$

This may also be written as a summation, which may help to visualize the process of writing the terms of the Taylor approximation. The Taylor approximation as a summation becomes

$$ f(x) \approx \sum_{n=0}^{\infty} \frac{(x-x_0)^n}{n!} f^{(n)}(x_0) $$

where $f^{(n)}(x_0)$ indicates the $n$-th derivative of a function $f(x)$ evaluated at $x=x_0$.

Task 2.0: Derive the backward difference method such that it is second order accurate. Refer to the book for an illustration of this approach with first order accuracy. Insert an image of your math below.

(You don't have to use this formula later on in this assignment, but it is useful to understand the approach and preparing for the exam)

You will be asked to include your derivation in the Report, but you can also include it here.

Solution:

For backward euler we need to evaluate at $x_{i-1}$ around $x_{i}$ and $x_{i-2}$ around $x_{i}$. This yields the following Taylor approximation:

$$ f(x_{i-1})\approx f(x_{i})+(x_{i-1}-x_i)\frac{\partial f(x_{i})}{\partial x} +\frac{(x_{i-1}-x_i)^2}{2}\frac{\partial^2 f(x_i)}{\partial x^2}+\mathcal{O}(x_{i-1}-x_i)^3 \newline\newline f(x_{i-2})\approx f(x_{i})+(x_{i-2}-x_i)\frac{\partial f(x_{i})}{\partial x} +\frac{(x_{i-2}-x_i)^2}{2}\frac{\partial^2 f(x_i)}{\partial x^2}+\mathcal{O}(x_{i-2}-x_i)^3 \newline\newline $$

we set $\Delta x = x_i - x_{i-1}$, which also means: $2\Delta x = x_i - x_{i-2}$ $$ f(x_{i-1})\approx f(x_{i})-\Delta x\frac{\partial f(x_{i})}{\partial x} +\frac{\Delta x^2}{2}\frac{\partial^2 f(x_i)}{\partial x^2}+\mathcal{O}(\Delta x)^3 \newline\newline f(x_{i-2})\approx f(x_{i})-2\Delta x\frac{\partial f(x_{i})}{\partial x} +\frac{4\Delta x^2}{2}\frac{\partial^2 f(x_i)}{\partial x^2}+\mathcal{O}(\Delta x)^3 \newline\newline $$ We multiply the first expression by 4 and subtract the second expression:

$$ 4f(x_{i-1})-f(x_{i-2})\approx (4-2)f(x_{i})-(4-2)\Delta x\frac{\partial f(x_{i})}{\partial x} +(4-4)\frac{\Delta x^2}{2}\frac{\partial^2 f(x_i)}{\partial x^2}+\mathcal{O}(\Delta x)^3 \newline\newline 3f(x_i)-4f(x_{i-1})+f(x_{i-2}) \approx 2\Delta x\frac{\partial f(x_{i})}{\partial x} +\mathcal{O}(\Delta x)^3 \newline \frac{3f(x_i)-4f(x_{i-1})+f(x_{i-2})}{2\Delta x} \approx \frac{\partial f(x_{i})}{\partial x}+\mathcal{O}(\Delta x)^2 $$

We don't care if you use AI tools for derivations, but remember: deriving equations in this way will be necessary to complete exam questions, so you really should practice doing it by hand!

Task 2.1: Derive the Taylor series expansion terms:

On paper, obtain the first four derivatives for the expression: $$ f(x)= 2\cos(x)+\sin(x) $$

Once you obtain the expressions, evaluate them around the point $x_0=\pi$.

These terms will be used later to assess the effects of using more or less terms in the approximation as well as the distance from $x_0$ later on in the notebook.

Use the following markdown cell to include your derivation of the Taylor series terms.

You will be asked to include your derivation in the Report, but you can also include it here.

Solution:

Find the first four derivatives of $f(x)$:

• $f'(x) = -2\sin(x) + \cos(x)$
• $f''(x) = -2\cos(x) - \sin(x)$
• $f'''(x) = 2\sin(x) - \cos(x)$
• $f^{(4)}(x) = 2\cos(x) + \sin(x)$

Evaluate the terms at $x_0=\pi$

• $f'(\pi) = -2\sin(\pi) + \cos(\pi) = 0 + (-1) = -1$
• $f''(\pi) = -2\cos(\pi) - \sin(\pi) = -2(-1) - 0 = 2$
• $f'''(\pi) = 2\sin(\pi) - \cos(\pi) = 0 - (-1) = 1$
• $f^{(4)}(\pi) = 2\cos(\pi) + \sin(\pi) = 2(-1) + 0 = -2$

Task 2.2: Plotting reference expression and approximations.

Before continuing with TSE, plot the expression $f(x)=2\cos(x)+\sin(x)$ in the interval $[-3\pi,5\pi]$. This will be used as benchmark to assess your approximations. We will want to produce plots that will include each successive term of the Taylor approximation to see how the approximation improves as we include more terms in the Taylor series.

# x = np.linspace(-3*np.pi, 5*np.pi, 400)

# def f(x):
#     return YOUR_CODE_HERE

# plt.plot(YOUR_CODE_HERE, YOUR_CODE_HERE, color='b', linewidth=2)
# plt.xlabel('x')
# plt.ylabel('f(x)')
# plt.title("Plot of $f(x) = 2cos(x) + sin(x)$");

# SOLUTION:
x = np.linspace(-3*np.pi, 5*np.pi, 400)

def f(x):
    return 2 * np.cos(x) + np.sin(x)

plt.plot(x, f(x), color='b', linewidth=2)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title("Plot of $f(x) = 2cos(x) + sin(x)$");

Task 2.3: Putting pen to python.

Enter the first four derivative terms in the corresponding functions below so that we can use them in our approximation.

Complete the functions in the code cell below, f_# refers to the derivatives of each portion of the function.

# def f_1(x):
#     return YOUR_CODE_HERE

# def f_2(x):
#     return YOUR_CODE_HERE

# def f_3(x):
#     return YOUR_CODE_HERE

# def f_4(x):
#     return YOUR_CODE_HERE

# SOLUTION:
def f_1(x):
    return -2 * np.sin(x) + np.cos(x)

def f_2(x):
    return -2 * np.cos(x) - np.sin(x)

def f_3(x):
    return 2 * np.sin(x) - np.cos(x)

def f_4(x):
    return 2 * np.cos(x) + np.sin(x)

Task 2.4: Define the expansion point and write the Taylor series expansion of first, second, third and fourth order.

Complete the functions in the code cell below, the taylor_# refers to the order of the Taylor approximation and x0 the expansion point, using our previously defined $x_0$.

# x0 = YOUR_CODE_HERE
# taylor_1 = YOUR_CODE_HERE
# taylor_2 = YOUR_CODE_HERE
# taylor_3 = YOUR_CODE_HERE
# taylor_4 = YOUR_CODE_HERE

x0 = np.pi
taylor_1 = f(x0) + f_1(x0)*(x - x0)
taylor_2 = taylor_1 + f_2(x0)*(x - x0)**2/2
taylor_3 = taylor_2 + f_3(x0)*(x - x0)**3/6
taylor_4 = taylor_3 + f_4(x0)*(x - x0)**4/24

Task 2.5: Plot your function along with your Taylor orders to illustrate the local approximations of the inclusions of each extra term. Add your own labels based on the formatting of the first example plot.

plt.figure(figsize=(10, 6))

# plt.plot(YOUR_CODE_HERE, YOUR_CODE_HERE,
#          label='$f(x) = 2cos(x) + sin(x)$', color='b', linewidth=2)
# plt.plot(YOUR_CODE_HERE, YOUR_CODE_HERE,
#          label='First Order', linestyle='--', color='g', linewidth=2)
# plt.plot(YOUR_CODE_HERE, YOUR_CODE_HERE,
#          label='Second Order', linestyle='--', color='r', linewidth=2)
# plt.plot(YOUR_CODE_HERE, YOUR_CODE_HERE,
#          label='Third Order', linestyle='--', color='m', linewidth=2)
# plt.plot(YOUR_CODE_HERE, YOUR_CODE_HERE,
#          label='Fourth Order', linestyle='--', color='y', linewidth=2)

# SOLUTION:
plt.plot(x, f(x),
         label='$f(x) = 2cos(x) + sin(x)$', color='b', linewidth=2)
plt.plot(x, taylor_1,
         label='First Order', linestyle='--', color='g', linewidth=2)
plt.plot(x, taylor_2,
         label='Second Order', linestyle='--', color='r', linewidth=2)
plt.plot(x, taylor_3,
         label='Third Order', linestyle='--', color='m', linewidth=2)
plt.plot(x, taylor_4,
         label='Fourth Order', linestyle='--', color='y', linewidth=2)


plt.scatter([x0], [f(x0)],
            color='k', marker='o', label=f'Expansion Point ($x = {x0:0.3f}$)')

plt.xlabel('x')
plt.ylabel('$f(x)$')
plt.title('Taylor Series Expansion of $f(x) = 2cos(x) + sin(x)$')
plt.legend()
plt.xlim(-1,9)
plt.ylim(-10,10)

plt.grid(True)
plt.show();

Estimating the Truncation error

When we use Taylor series approximations of a function, we truncate the infinite series to a specfic number of terms. Thus, a truncation error is introduced. A simple error metric between the analytical solution and the numerical approximation of our function can be found by computing the absolute value of the difference between the function and the approximation, namely $$ \text{error }=|f(x)-T_n| $$ where $T_n$ refers to the TSE computed using $n$ number of terms (or derivatives) of the Taylor series we use.

Task 2.6: Use your Taylor series approximations and the analytic solution to determine their absolute error. Plot these approximations and vary the x and y limits. Are the larger orders TSE always better?

# error_1 = YOUR_CODE_HERE
# error_2 = YOUR_CODE_HERE
# error_3 = YOUR_CODE_HERE
# error_4 = YOUR_CODE_HERE

# SOLUTION:
error_1 = np.abs(f(x) - taylor_1)
error_2 = np.abs(f(x) - taylor_2)
error_3 = np.abs(f(x) - taylor_3)
error_4 = np.abs(f(x) - taylor_4)

plt.figure(figsize=(10, 4))
plt.plot(x, error_1,
         label='First Order', color='g', linewidth=2)
plt.plot(x, error_2,
         label='Second Order', color='r', linewidth=2)
plt.plot(x, error_3,
         label='Third Order', color='m', linewidth=2)
plt.plot(x, error_4,
         label='Fourth Order', color='y', linewidth=2)

plt.xlabel('x')
plt.ylabel('Absolute Error: $f(x)$-Taylor Order')
plt.title('Absolute Error of Taylor Series Approximations')
plt.xlim(np.pi-1,np.pi+1)
plt.ylim(0,0.01)
plt.legend()

plt.grid(True)
plt.show();

Part 3: Taylor Series Expansion in Two Variables

Let's investigate how the Taylor Series Expansion operates in two dimensions. Expanding a function dependent on two variables, $f(x, y)$, can be expressed as:

$$ f(x, y) \approx f(x_0, y_0) + (x - x_0)\frac{\partial f}{\partial x}\bigg|_{(x_0, y_0)} + (y - y_0) \frac{\partial f}{\partial y}\bigg|_{(x_0, y_0)} \\ + \frac{(x - x_0)^2}{2!}\frac{\partial^2 f}{\partial x^2}\bigg|_{(x_0, y_0)} + \frac{(y - y_0)^2}{2!}\frac{\partial^2 f}{\partial y^2}\bigg|_{(x_0, y_0)} \\ + (x - x_0)(y - y_0)\frac{\partial^2 f}{\partial x \partial y}\bigg|_{(x_0, y_0)} + \ldots $$

The terms in this summation are determined by the partial derivatives of $f$ with respect to $x$ and $y$ evaluated at $(x_0, y_0)$. The expansion includes an infinite series of terms, starting from the first-order terms and continuing to higher-order terms as $n$ increases. The choice of how many terms to include in the expansion depends on the desired level of accuracy.

Task 3.1: Writing out the expansion.

Using the function $f(x,y)=\sin(2x)+\cos(y)$ derive the first and second derivatives of the approximation. Do this on paper, as before. After obtaining derivatives, evaluate them around the expanding point $(x_0=\pi,y_0=\pi)$. Also, add the evaluation of $f(x_0,y_0)$.

Use the following cell to include your workings for the solution.

You will be asked to include your derivation in the Report, but you can also include it here.

Solution:

Compute the derivatives:

• $\frac{\partial f}{\partial x}=2 \cos (2x)$
• $\frac{\partial f}{\partial y}=- \sin (y)$
• $\frac{\partial ^2 f}{\partial x^2}=-4 \sin (2x)$
• $\frac{\partial ^2 f}{\partial y^2}=-\cos(y)$
• $\frac{\partial ^2 f}{\partial x \partial y}=0$

Evaluate Derivatives at point $(x_0,y_0)=(\pi, \pi)$

Evaluate $(f(\pi, \pi))$ and its first and second partial derivatives at $(\pi, \pi)$:

• $f(\pi, \pi) = \sin(2\pi) + \cos(\pi) = 0 - 1 = -1$
• $\frac{\partial f}{\partial x}(\pi, \pi) = 2\cos(2\pi) = 2$
• $\frac{\partial f}{\partial y}(\pi, \pi) = -\sin(\pi) = 0$
• $\frac{\partial^2 f}{\partial x^2}(\pi, \pi) = -4\sin(2\pi)= 0$
• $\frac{\partial^2 f}{\partial y^2}(\pi, \pi) = -\cos(\pi) = 1$
• $\frac{\partial^2 f}{\partial x \partial y}(\pi, \pi) = 0$

Write out the second-order Taylor series expansion

The second-order Taylor series expansion for $f(x, y)$ about the point $(\pi, \pi)$ is given by:

$$ f(x, y) \approx f(\pi, \pi) + \frac{\partial f}{\partial x}(\pi, \pi)(x - \pi) + \frac{\partial f}{\partial y}(\pi, \pi)(y - \pi) + \frac{1}{2}\frac{\partial^2 f}{\partial x^2}(\pi, \pi)(x - \pi)^2 + \frac{1}{2}\frac{\partial^2 f}{\partial y^2}(\pi, \pi)(y - \pi)^2 $$$$ f(x, y) \approx -1 + 2(x - \pi) +\frac{(y - \pi)^2}{2} $$

Now, we can substitute the evaluated values from Step 2 into this formula.

Task 3.2: Transfer the Taylor Approximation

Enter the expansion point and include the Taylor approximation in the function definition of the code cell below.

Remember, to break a line in the code and keep it readable you can surround the expression with an extra set of parentheses, ( ); this will allow arbitrarily-placed line breaks within the code.

Note that here we define x0 and y0 as global variables (outside the function). This is a bad practice in general, but in this case it is OK because we are limiting our use of the code to a few small functions and tasks in a single notebook.

# def f2D(x, y):
#     return YOUR_CODE_HERE

# x0 = YOUR_CODE_HERE
# y0 = YOUR_CODE_HERE

# def taylor2D(x, y):
#     return (YOUR_CODE_HERE)


# SOLUTION:
def f2D(x, y):
    return np.sin(2*x) + np.cos(y)

x0 = np.pi
y0 = np.pi

def taylor2D(x, y):
    return ((np.sin(2*x0)
            + np.cos(y0))
            + (2*np.cos(2*x0)*(x - x0))
            - (np.sin(y0)*(y-y0))
            - (1/2)*(4*np.sin(2*x0))*(x - x0)**2
            - (1/2)*(np.cos(y0))*(y - y0)**2)

# this is a second way to implement the function above
def taylor2D(x, y):
    return -1 + 2*(x - x0) + (y - x0)**2/2

Visualize Results

Let's plot a 3D surface plot of the results and try to visualize the intersection of the function and its Taylor approximation. In order to do this, we need to create a meshgrid over the region of interest (this is done with np.meshgrid). Using our def functions above we can determine the surface.

Task 3.3: Scrutinize the code to understand what is being plotted, then run it to visualize the results.

# Create a meshgrid of x and y values
x = np.linspace(-2+x0, 2+x0, 100)
y = np.linspace(-2+y0, 2+y0, 100)
X, Y = np.meshgrid(x, y)

# Calculate the original function values and the approximation
Z_orig = f2D(X, Y)
Z_approx = taylor2D(X, Y)

# Create a 3D plot
fig = plt.figure(figsize=(9,9))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z_orig, rstride=1, cstride=1, cmap='Reds',
                label='Original Function')
ax.plot_surface(X, Y, Z_approx, rstride=1, cstride=1, cmap='Blues',
                alpha=0.7, label='Taylor Approximation')

# Set labels and legend
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('f(x, y)')
plt.legend()

# Show the plot
plt.title('Original Function vs. Taylor Approximation')
plt.show()

c:\Users\beren\anaconda3\envs\TAMude\Lib\site-packages\IPython\core\pylabtools.py:170: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  fig.canvas.print_figure(bytes_io, **kw)

Error Analysis

How good is this approximation?

Task 3.4: Use the following cell to calculate the absolute error between the analytical and Taylor approximation as previously and plot the results.

# error_2d = YOUR_CODE_HERE

# SOLUTION:
error_2d = np.abs(Z_orig - Z_approx)

fig = plt.figure(figsize=(9,9))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, error_2d, cmap='Reds', label='Absolute Error')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Absolute Error')
plt.title('Absolute Error between $f(x, y)$ and Taylor Approximation')
plt.show()

End of notebook.

By MUDE Team © 2024 TU Delft. CC BY 4.0. doi: 10.5281/zenodo.16782515.