# Differentiation and Integration

This script is a tutorial on how to differentiate and integrate numerically. No fancy techniques are used and it is really an exercise in using the data structures learned in the first lesson.

## Define parameters and set up the grid

```numpts=1000;
L=10;
```

Set up the 1d grid.

```xg=10*linspace(0,1,numpts+1); dx=xg(2)-xg(1);
xgm=0.5*(xg(2:numpts+1)+xg(1:numpts));
dxg=(xg(2:numpts+1)-xg(1:numpts));
```

## Define the function

Define the function of interest, its derivative and integral (I used Maple).

```myf=@(x) x.*(10-x).*sin(4*pi*x/L);
myfp=@(x) (10-x).*sin((2/5)*pi.*x)-x.*sin((2/5)*pi.*x)+(2/5)*x.*(10-x).*cos((2/5)*pi.*x)*pi;
myintf=@(x) (25/4)*(10*sin((2/5).*pi.*x)-4*x.*cos((2/5)*pi*x)*pi-(5/2)*(-(4/25)*pi^2*x.^2.*cos((2/5)*pi*x)+2*cos((2/5)*pi*x)+(4/5)*pi*x.*sin((2/5)*pi*x))/pi)/pi^2;
myg=@(x) 10*x-0.1*3*x.^2;
myintg=@(x) 5*x.^2-0.1*x.^3;
```

Get the function and the derivative on the grid.

```fg=myf(xg);
```

For the derivative use the grid at the mid points.

```fgp=myfp(xgm);
```

## Finite Difference approximation

Now use a finite difference approximation that uses the info at the closest point on the right and the closest point on the left.

```fpmethod1=(fg(2:end)-fg(1:end-1))./dxg;
figure(1)
clf
subplot(2,1,1)
plot(xgm,fgp,'b',xgm,fpmethod1,'ro')
xlabel('x')
ylabel('derivative')
subplot(2,1,2)
semilogy(xgm,abs(fpmethod1-fgp))
xlabel('x')
ylabel('error')
``` Try changing numpts to see how the error changes. The finite difference method in the above is of second order (you can check this from the Taylor series) but uses a different grid. We will talk about generating differentiation matrices in a different M-file but for now let's try to use interpolation onto the original grid to see how that effects the error.

```fporig1=interp1(xgm,fpmethod1,xg,'linear');
fporig2=interp1(xgm,fpmethod1,xg,'pchip');
fporig3=interp1(xgm,fpmethod1,xg,'spline');
```

The moral of the story is that the built in spline interpolation is both quick and accurate.

```fgp0=myfp(xg);
mxerrlin=max(abs(fporig1-fgp0))
mxerrcubic=max(abs(fporig2-fgp0))
mxerrspline=max(abs(fporig3-fgp0))
```
```mxerrlin =

9.5249e-04

mxerrcubic =

9.5249e-04

mxerrspline =

2.3813e-04

```

## Integration

```% Integration
gg=myg(xg);
% First order estimate of integral from the left and right Riemann sums.
myint1l=sum(gg(1:end-1))*dx;
myint1r=sum(gg(2:end))*dx;
```

Second order estimate of integral from the trapezoid rule. Notice that for not a lot of extra work we get a real improvement.

```myint2=(0.5*gg(1)+sum(gg(2:end-1))+0.5*gg(end))*dx;
myinta=myintg(L)-myintg(0);
relerr1l=(myint1l-myinta)/myinta
relerr1r=(myint1r-myinta)/myinta
relerr2=(myint2-myinta)/myinta
```
```relerr1l =

-8.7513e-04

relerr1r =

8.7487e-04

relerr2 =

-1.2500e-07

```