18.2: Numerical Differentiation
- Page ID
- 86415
\( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)
\( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)
\( \newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\)
( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\)
\( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\)
\( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\)
\( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\)
\( \newcommand{\Span}{\mathrm{span}}\)
\( \newcommand{\id}{\mathrm{id}}\)
\( \newcommand{\Span}{\mathrm{span}}\)
\( \newcommand{\kernel}{\mathrm{null}\,}\)
\( \newcommand{\range}{\mathrm{range}\,}\)
\( \newcommand{\RealPart}{\mathrm{Re}}\)
\( \newcommand{\ImaginaryPart}{\mathrm{Im}}\)
\( \newcommand{\Argument}{\mathrm{Arg}}\)
\( \newcommand{\norm}[1]{\| #1 \|}\)
\( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\)
\( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\AA}{\unicode[.8,0]{x212B}}\)
\( \newcommand{\vectorA}[1]{\vec{#1}} % arrow\)
\( \newcommand{\vectorAt}[1]{\vec{\text{#1}}} % arrow\)
\( \newcommand{\vectorB}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)
\( \newcommand{\vectorC}[1]{\textbf{#1}} \)
\( \newcommand{\vectorD}[1]{\overrightarrow{#1}} \)
\( \newcommand{\vectorDt}[1]{\overrightarrow{\text{#1}}} \)
\( \newcommand{\vectE}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash{\mathbf {#1}}}} \)
\( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)
\( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)
\(\newcommand{\avec}{\mathbf a}\) \(\newcommand{\bvec}{\mathbf b}\) \(\newcommand{\cvec}{\mathbf c}\) \(\newcommand{\dvec}{\mathbf d}\) \(\newcommand{\dtil}{\widetilde{\mathbf d}}\) \(\newcommand{\evec}{\mathbf e}\) \(\newcommand{\fvec}{\mathbf f}\) \(\newcommand{\nvec}{\mathbf n}\) \(\newcommand{\pvec}{\mathbf p}\) \(\newcommand{\qvec}{\mathbf q}\) \(\newcommand{\svec}{\mathbf s}\) \(\newcommand{\tvec}{\mathbf t}\) \(\newcommand{\uvec}{\mathbf u}\) \(\newcommand{\vvec}{\mathbf v}\) \(\newcommand{\wvec}{\mathbf w}\) \(\newcommand{\xvec}{\mathbf x}\) \(\newcommand{\yvec}{\mathbf y}\) \(\newcommand{\zvec}{\mathbf z}\) \(\newcommand{\rvec}{\mathbf r}\) \(\newcommand{\mvec}{\mathbf m}\) \(\newcommand{\zerovec}{\mathbf 0}\) \(\newcommand{\onevec}{\mathbf 1}\) \(\newcommand{\real}{\mathbb R}\) \(\newcommand{\twovec}[2]{\left[\begin{array}{r}#1 \\ #2 \end{array}\right]}\) \(\newcommand{\ctwovec}[2]{\left[\begin{array}{c}#1 \\ #2 \end{array}\right]}\) \(\newcommand{\threevec}[3]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \end{array}\right]}\) \(\newcommand{\cthreevec}[3]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \end{array}\right]}\) \(\newcommand{\fourvec}[4]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \\ #4 \end{array}\right]}\) \(\newcommand{\cfourvec}[4]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \\ #4 \end{array}\right]}\) \(\newcommand{\fivevec}[5]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \\ #4 \\ #5 \\ \end{array}\right]}\) \(\newcommand{\cfivevec}[5]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \\ #4 \\ #5 \\ \end{array}\right]}\) \(\newcommand{\mattwo}[4]{\left[\begin{array}{rr}#1 \amp #2 \\ #3 \amp #4 \\ \end{array}\right]}\) \(\newcommand{\laspan}[1]{\text{Span}\{#1\}}\) \(\newcommand{\bcal}{\cal B}\) \(\newcommand{\ccal}{\cal C}\) \(\newcommand{\scal}{\cal S}\) \(\newcommand{\wcal}{\cal W}\) \(\newcommand{\ecal}{\cal E}\) \(\newcommand{\coords}[2]{\left\{#1\right\}_{#2}}\) \(\newcommand{\gray}[1]{\color{gray}{#1}}\) \(\newcommand{\lgray}[1]{\color{lightgray}{#1}}\) \(\newcommand{\rank}{\operatorname{rank}}\) \(\newcommand{\row}{\text{Row}}\) \(\newcommand{\col}{\text{Col}}\) \(\renewcommand{\row}{\text{Row}}\) \(\newcommand{\nul}{\text{Nul}}\) \(\newcommand{\var}{\text{Var}}\) \(\newcommand{\corr}{\text{corr}}\) \(\newcommand{\len}[1]{\left|#1\right|}\) \(\newcommand{\bbar}{\overline{\bvec}}\) \(\newcommand{\bhat}{\widehat{\bvec}}\) \(\newcommand{\bperp}{\bvec^\perp}\) \(\newcommand{\xhat}{\widehat{\xvec}}\) \(\newcommand{\vhat}{\widehat{\vvec}}\) \(\newcommand{\uhat}{\widehat{\uvec}}\) \(\newcommand{\what}{\widehat{\wvec}}\) \(\newcommand{\Sighat}{\widehat{\Sigma}}\) \(\newcommand{\lt}{<}\) \(\newcommand{\gt}{>}\) \(\newcommand{\amp}{&}\) \(\definecolor{fillinmathshade}{gray}{0.9}\)1-Dimensional Derivatives
In calculus, you learned that the derivative of a function is: limit of Δy/Δx as x goes to 0.
Sometimes, we only have measured data, not an analytic function, so we cannot get the derivative analytically. When the analytic derivative is unknown, we can approximate it from the data using Δy/Δx, without taking a limit. The numerical derivative at point (xn, yn) is:
(xn, yn) ≈ (yn+1 – yn-1) / (xn+1 - xn-1,)
Matlab’s function for this is named gradient(). There are 2 methods of calling gradient(). Both methods take 2 arguments, and the 1st argument is the vector of y-values. The 2nd argument depends on the method:
Method A: dydx_gradient = gradient(y, x) where the 2nd argument is the vector of the x values corresponding to y.
Note: Both y and x need to be inputs.
When y is the only input, gradient(y) assumes that the x spacing = 1, which gives the wrong answer by a factor of 10!
Method B. dydx_gradient = gradient(y, dx), where the 2nd argument is the spacing, dx, between the x values (dx = a single number).
Both methods give the same results when the x values are equally spaced.
dx = 0.05*pi
x = 0: dx: 0.5*pi;
y = sin(x);
%% First, plot y vs. x
figure;
plot(x,y,'-o', 'LineWidth',2);
grid on;
title('y = sin(x)')
xlabel('x (rads)')
%% Label the (x, y) values of 3 consecutive points
text(x(5), (y(5)-0.02),['y= ',num2str(y(5),3)])
text(x(6), (y(6)-0.02),['y= ',num2str(y(6),3)])
text(x(7), (y(7)-0.02),['y= ',num2str(y(7),3)])
%% Second, plot the analytic derivative
dy_anal = cos(x); % This is the analytic derivative
figure;
%subplot(2,1,1);
plot(x,dy_anal,'s','MarkerSize',10);
grid on;
xlabel('x (rads)')
%% Approximate the derivative with the gradient function: grad()
dydx_gradientA = gradient(y, x) % Method A
dydx_gradientB = gradient(y, dx) % Method B
dydx_gradientA =
Columns 1 through 4
0.9959 0.9836 0.9472 0.8873
Columns 5 through 8
0.8057 0.7042 0.5854 0.4521
Columns 9 through 11
0.3077 0.1558 0.0784
dydx_gradientB =
Columns 1 through 4
0.9959 0.9836 0.9472 0.8873
Columns 5 through 8
0.8057 0.7042 0.5854 0.4521
Columns 9 through 11
0.3077 0.1558 0.0784
hold on;
plot(x, dydx_gradientA, '*', 'MarkerSize',10);
legend('dy\_anal\_cos','dy\_gradient')
title({'Numerical Gradient Approximates',
'the Analytic Derivative of sin(x)'})
Solution
.
Figure \(\PageIndex{1}\): y = sin(x) |
Figure \(\PageIndex{1}\): Derivatives of y = sin(x) The numerical derivative is very close to the analytic values, except for the last point, because there is no point beyond the last point to improve the estimated value. |
Add example text here.
.
Homework:
Download the attached norm.m file which computes the normal distribution:
% Start a script .m file with this line:
clear all; close all; clc; format compact
% Define this vector
x = -3: 0.1 : 3;
% Compute y:
y = norm_dist(x);
%% Open a figure and plot y as a function of x:
% Turn on the grid:
%% Use gradient() to compute the derivative.
%% Then plot the gradient result on the same figure;
% If computed correctly:
% the derivative will be > 0 for x < 0
% the derivative will be = 0 for x = 0
% the derivative will be < 0 for x > 0
% Use the legend function to identify the 2 curves:
legend('norm\_dist', 'derivative')
% Add this title:
title('norm_dist and its derivative')
- Answer
-
The answer is not given here.
.
2 Dimensional Derivatives
Matlab’s gradient function can also be used to compute the gradients of a function defined on x and y grids created by meshgrid. This example uses Matlab’s built-in peaks() 3D function.
%% Numerical Derivatives of the Peaks 2D function
clear all; close all; clc; format compact
dx = 0.5;
x1D = (-3: dx : 3);
y1D = x1D;
[x2Da, y2Da] = meshgrid(x1D,y1D);
[x2D, y2D, z2D] = peaks(x2Da,y2Da);
figure;
mesh(x2D, y2D, z2D)
title('Peaks')
colorbar
%%
[dzdx, dzdy] = gradient(z2D, x1D, y1D);
figure;
contour(x2D, y2D, z2D);
axis equal
hold on;
% Draw the x & y cradients as arrows with the quiver plot function
quiver(x2D, y2D, dzdx, dzdy);
title('peaks countour and x and y derivatives (Fig. 13.27b)')
colorbar
Solution
The resulting plots are:
![]() |
![]() |
.