% 4-29-08
% Differentiaion Demo

%% Basic Derivative examples
clear all; close all; figure(1);

dx = 0.1;
x = -10:dx:10;

y = sech(x);
yx_exact = -sech(x).*tanh(x);

plot(x,y,'b-',x,yx_exact,'r-');
hold on;

n = length(x);

% Second Order
yx(1) = (-3*y(1)+4*y(2)-y(3))/(2*dx);
for i=2:n-1
    yx(i) = (y(i+1)-y(i-1))/(2*dx);
end
yx(n) = (3*y(n) - 4*y(n-1) + y(n-2))/(2*dx);
plot(x,yx,'ko');

% Fourth Order
yx(1) = (-3*y(1)+4*y(2)-y(3))/(2*dx);
yx(2) = (y(3)-y(1))/(2*dx);
for i=3:n-2
    yx(i) = (-y(i+2)+8*y(i+1)-8*y(i-1)+y(i-2))/(12*dx);
end
yx(n-1) = (y(n)-y(n-2))/(2*dx);
yx(n) = (3*y(n) - 4*y(n-1) + y(n-2))/(2*dx);
plot(x,yx,'go');
hold off;

%% Date differentiation example
clear all; figure(2);
load spring.dat
x=spring(:,1);
u=spring(:,2);
dx=mean(diff(x));
ux_exact = -2*x.*exp(-x.^2).*sin(10*x).^2+20*exp(-x.^2).*sin(10*x).*cos(10*x);
plot(x,u,'ko',x,ux_exact,'go');
hold on;

n=length(x);

ux(1) = (-3*u(1) + 4*u(2) -u(3) ) / (2*dx);
for i=2:n-1
    ux(i) = (u(i+1) - u(i-1))/(2*dx);
end
ux(n) = (3*u(i) - 4*u(i-1) + u(i-2))/(2*dx);

plot(x,ux,'r')

dx=0.02;
xx = 0:dx:1;

u_interp = spline(x,u,xx);

n=length(xx);

ux(1) = (-3*u_interp(1) + 4*u_interp(2) -u_interp(3) ) / (2*dx);
for i=2:n-1
    ux(i) = (u_interp(i+1) - u_interp(i-1))/(2*dx);
end
ux(n) = (3*u_interp(i) - 4*u_interp(i-1) + u_interp(i-2))/(2*dx);

plot(xx,ux,'k')
plot(xx,u_interp,'c')
hold off;

%% Discontinuous Function Example
clear all; figure(3);

dx = 0.5;
x=-10:dx:10;

u=(x<0).*(-1/10*x+1) + (x>=0).*(1/10*x+1);
ux_exact = (x<0).*-1/10 + (x>=0).*1/10;
plot(x,u,'b',x,ux_exact,'g')
hold on;

n=length(x);

ux(1) = (-3*u(1) + 4*u(2) -u(3) ) / (2*dx);
for i=2:n-1
    ux(i) = (u(i+1) - u(i-1))/(2*dx);
end
ux(n) = (3*u(i) - 4*u(i-1) + u(i-2))/(2*dx);

plot(x,ux,'or')

ux(1) = (-3*u(1) + 4*u(2) -u(3) ) / (2*dx);
ux(2) = (u(3) - u(1))/(2*dx);
for i=3:n-2
    ux(i) = (-u(i+2) + 8*u(i+1) - 8*u(i-1) + u(i-2))/(12*dx);
end
ux(n-1) = (u(n) - u(n-2))/(2*dx);
ux(n) = (3*u(i) - 4*u(i-1) + u(i-2))/(2*dx);

plot(x,ux,'ok')
hold off;