Derivative Matrix for a Function Vector

Author

John Hawkins

Problem

The following problem statement is that proposed by Prof. Frohne on the main class page:

1. Explore how a linear operator, like for example $d \over dt$ can be represented as some kind of a matrix multiply (with perhaps an infinite number of dimensions).

Solution

The solution presented here is limited to presenting a matrix representation of the operator $d/dx$.

First, let us consider the finite case. Let $x$ be the independent variable with $n$ components and let the components of $x$ be defined such that $x_{m+1}=x_{m}+\Delta x$ for all $m<n$ and for some constant $\Delta x$. Also let $y=f(x)$. We wish to find the slope of the function at $y(x_{m})$ for any $m$. To do this for the finite case, we will approximate the slope at a point $y(x_{m})$ as being the slope of the line between the two points $y(x_{m-1})$ and $y(x_{m+1})$ given that $1<m<n$. For the case when $m=1$ or $m=n$, we approximate $dy/dx$ as the slope the lines between $y(x_{1})$ and $y(x_{n})$ respectively. These approximations are accomplished by the following matrix:

${\frac {d}{dx}}={\dfrac {1}{2(\Delta x)}}{\begin{bmatrix}-2&2&0&0&&&&&\\-1&0&1&0&&&\cdots &&\\0&-1&0&1&&&&&\\\\&&&&\ddots &&&&\\\\&\vdots &&&&-1&0&1&0\\&&&&&0&-1&0&1\\&&&&&0&0&-2&2\end{bmatrix}}$

Note that we need to divide by $1/2(\Delta x)$ because the distance between $x_{m-1}$ and $x_{m+1}$ is $2\Delta x$. Similarly, the nonzero entries in the first and last rows are 2 to cancel out the 2 in the denominator. Using this matrix, we find that

$\left({\frac {d}{dx}}\right)y={\dfrac {1}{2(\Delta x)}}{\begin{bmatrix}-2&2&0&0&&&&&\\-1&0&1&0&&&\cdots &&\\0&-1&0&1&&&&&\\\\&&&&\ddots &&&&\\\\&\vdots &&&&-1&0&1&0\\&&&&&0&-1&0&1\\&&&&&0&0&-2&2\end{bmatrix}}{\begin{bmatrix}y(x_{1})\\y(x_{2})\\y(x_{3})\\\\\vdots \\\\y(x_{n-2})\\y(x_{n-1})\\y(x_{n})\end{bmatrix}}={\begin{bmatrix}\left[y(x_{2})-y(x_{1})\right]/\Delta x\\\left[y(x_{3})-y(x_{1})\right]/2(\Delta x)\\\left[y(x_{4})-y(x_{2})\right]/2(\Delta x)\\\\\vdots \\\\\left[y(x_{n-1})-y(x_{n-3})\right]/2(\Delta x)\\\left[y(x_{n})-y(x_{n-2})\right]/2(\Delta x)\\\left[y(x_{n})-y(x_{n-1})\right]/\Delta x\end{bmatrix}}$

which corresponds with the slopes we wished to obtain.

Solution Picture

The figure below was created to demonstrate how the above matrix works. First, 10 points of the function are displayed. Then for each point the line between the two adjacent points is drawn and the slope of that line is shown to be $dy/dx$ for the original point, calculated using the matrix shown above. After this has been done in the case of 10 points, the function and derivative vectors are shown for 15, 20, 30, 50, 100, 200, and 500 points.

Code

The figure above was generated with the following code:

<source lang="matlab">

%

% _-_,, ,, _-_- ,,

% ( // || /, _  ; || '

% _|| /'\\ ||/\\ \\/\\ || __ < \, \\/\/\ ||/\ \\ \\/\\ _-_,

% _|| || || || || || || ~||- - /-|| || | | ||_< || || || ||_.

% || || || || || || || ||===|| (( || || | | || | || || || ~ ||

% -__-, \\,/ \\ |/ \\ \\ ( \_, | \/\\ \\/\\/ \\,\ \\ \\ \\ ,-_-

% _/ `

%

% John Hawkins

% LNA

% Derivative Matrix Example

% 2 Feb 2010

function DerivativeExample

x=linspace(0,4*pi,500)';

y=sin(x);

dydx=cos(x);

h.fig=figure('color','w','units','normalized','position',[.2 .2 .6 .6]);

plot(x,y,'color',[.8 .8 .8]);

hold on;

plot(x,dydx,'color',[.8 .8 .8]);

numStartPoints=10;

x=linspace(0,4*pi,numStartPoints)';

y=sin(x);

dydx=DERIV(x,y);

h.y=plot(x,y,'b.');

h.construct=plot(0,0,'r');

h.dydx=plot(0,0,'r.');

legend([h.y,h.construct,h.dydx],'y','Construction Line','dy/dx');

xlabel('x'); ylabel('y');

axis equal;

title('Taking the Derivative of a Finite Function Vector for y=sin x',...

   'fontname','times','fontsize',20);

set([h.construct,h.dydx],'visible','off');

fig = getframe(h.fig);

[im,map] = rgb2ind(fig.cdata,256,'nodither');

set(h.construct,'visible','on','xdata',x([1,2]),'ydata',y([1,2]));

im=CAPTUREFIG(im,map,h);

set(h.dydx,'visible','on','xdata',x(1),'ydata',dydx(1));

im=CAPTUREFIG(im,map,h);

for n=2:numStartPoints-1

   set(h.construct,'xdata',x([n-1,n+1]),'ydata',y([n-1,n+1]));

   im=CAPTUREFIG(im,map,h);

   set(h.dydx,'xdata',x(1:n),'ydata',dydx(1:n));

   im=CAPTUREFIG(im,map,h);

end

set(h.construct,'xdata',x([numStartPoints-1,numStartPoints]),...

   'ydata',y([numStartPoints-1,numStartPoints]));

im=CAPTUREFIG(im,map,h);

set(h.dydx,'xdata',x,'ydata',dydx);

im=CAPTUREFIG(im,map,h);

set(h.construct,'visible','off');

im=CAPTUREFIG(im,map,h);

for k=1:5

   im=CAPTUREFIG(im,map,h);

end

for n=[15 20 30 50 100 200 500]

   x=linspace(0,4*pi,n)';

   y=sin(x);

   dydx=DERIV(x,y);

   set(h.y,'xdata',x,'ydata',y);

   set(h.dydx,'xdata',x,'ydata',dydx);

   im=CAPTUREFIG(im,map,h);

end

for k=1:5

   im=CAPTUREFIG(im,map,h);

end

imwrite(im,map,'Derivative.gif','DelayTime',.5,'LoopCount',inf)

function dydx=DERIV(x,y)

num=numel(y);

deltax=x(2)-x(1);

operator=zeros(num);

operator(1,1)=-2;

operator(1,2)=2;

for k=2:num-1

   operator(k,k+1)=1;

   operator(k,k-1)=-1;

end

operator(num,num-1)=-2;

operator(num,num)=2;

operator=operator/(2*deltax);

dydx=operator*y;

function im=CAPTUREFIG(im,map,h)

fig = getframe(h.fig);

im(:,:,1,size(im,4)+1) = rgb2ind(fig.cdata,map,'nodither');

</source>