Script 2 For Week 9

% Below are ten integration and differentiation functions which you may use for any problem. Lower ones sometimes call higher ones, so you should probably just create M-files for all 10 and put them all in a working directory.

% Here is the data from problem 21.15 to test the derivative M-files on.

t=[1 2 3.25 4.5 6 7 8 8.5 9.3 10]
v=[10 12 11 14 17 16 12 14 14 10]

plot(t,v)

% Using the functions below
%  we find the power k_1 so we can use Richardson 
% extrapolation on the trapezoid rule in richtrap.

ErrorPowerFinder(@(t) exp(t),0,8.1)

richtrap(@(t) exp(t),0,8.1)

% Here we first find the power k_1 so we can use Richardson 
% extrapolation on forward 1-step differentiation at a point.

ErrorPowerFinder2(@(t) exp(t),8.1)

RichDeriv(@(t) exp(t),8.1)

% Last we find the power k_1 so we can use Richardson 
% extrapolation on forward 1-step centered differentiation at a point.

ErrorPowerFinder3(@(t) exp(t),8.1)

RichDeriv2(@(t) exp(t),8.1)

% Immediately below is a bunch of my M-files, which could probably stand some improvement (error checking, etc.)
%*********************************************

%**********************************************************
function [E]=ErrorPowerFinder(fun,a,b)

% This script allows you to estimate or guess the order of the error term 
% in a numerical approximation method. You need this exponent to apply the 
% Richardson technique to try to improve order of the error of your 
% estimation method. 
% This one is for the trapezoid rule.


for k=0:10
t=linspace(a,b,6^k+1);
yt=fun(t);
A(k+1)=trapuneq(t,yt);
if k==1 break;
else
E(k+1)=log(abs((A(k+1)-A(k))/A(k)))/log(6);
end

end
%**********************************************************
%**********************************************************
function [E]=ErrorPowerFinder2(fun,a)

% This script allows you to estimate or guess the order of the error term 
% in a numerical approximation method. You need this exponent to apply the 
% Richardson technique to try to improve order of the error of your 
% estimation method. 
% This one is for the forward 1-step differentiation formula.


for k=1:6
h=1/2^k;
A(k)=(fun(a+h)-fun(a))/h;
end

for k=2:6
E(k-1)=log(abs((A(k)-A(k-1))/A(k)))/log(2);
end
%**********************************************************
%**********************************************************
function [E]=ErrorPowerFinder3(fun,a)

% This script allows you to estimate or guess the order of the error term 
% in a numerical approximation method. You need this exponent to apply the 
% Richardson technique to try to improve order of the error of your 
% estimation method. 
% This one is for the centered 1-step differentiation formula.


for k=1:6
h=1/2^k;
A(k)=(fun(a+h)-fun(a-h))/(2*h);
end

for k=2:6
E(k-1)=log(abs((A(k)-A(k-1))/A(k)))/log(2);
end

%**********************************************************
%**********************************************************
function [E R]=RichDeriv(fun,a)

% Looks at the forward difference estimate for the derivative

m=1/2;
for k=1:6

h=1/6^k;
A2(k)=(fun(a+h)-fun(a))/h;
A1(k)=(fun(a+h/2)-fun(a))/(h/2);
R(k)=(A1(k)-m*A2(k))/(1-m);
end

for k=2:6
E(1,k-1)=log(abs((A2(k)-A2(k-1))/A2(k)))/log(6);
E(2,k-1)=log(abs((A1(k)-A1(k-1))/A1(k)))/log(6);
E(3,k-1)=log(abs((R(k)-R(k-1))/R(k)))/log(6);
end



end
%**********************************************************
function [E R]=RichDeriv2(fun,a)

% Looks at the centered difference estimate for the derivative


m=1/2;
for k=1:6

h=1/6^k;
A2(k)=(fun(a+h)-fun(a-h))/(2*h);
A1(k)=(fun(a+h/2)-fun(a-h/2))/h;
R(k)=(A1(k)-m^2*A2(k))/(1-m^2);
end

for k=2:6
E(1,k-1)=log(abs((A2(k)-A2(k-1))/A2(k)))/log(6);
E(2,k-1)=log(abs((A1(k)-A1(k-1))/A1(k)))/log(6);
E(3,k-1)=log(abs((R(k)-R(k-1))/R(k)))/log(6);
end



end
%**********************************************************
%**********************************************************
function deriv=diff3uneq(x,y,t)

% t is a number
% x and y must be rows and have length 3
% the x entries must be in real-number order and not repeat
% This script creates the derivative of the Lagrange polynomial through three
% points and evaluates it at t. If t is one of the x values this corresponds
% to a backward, forward or central difference formula for the derivative.

if length(x)~=length(y)
disp('the two data vectors must have equal length') 

elseif length(x)~=3
disp('the data vectors must have length 3') 

elseif (sort(x)~=x) | (x(1)==x(2)) | (x(2)==x(3))
disp('the domain values must be in order and not repeat') 

elseif (length(t)~=1) | (sum(size(t))~=2)
disp('the third entry must be a number') 

else

deriv=(2*t-x(2)-x(3))/(x(1)-x(2))/(x(1)-x(3))*y(1) + (2*t-x(1)-x(3))/(x(2)-x(1))/(x(2)-x(3))*y(2) +(2*t-x(1)-x(2))/(x(3)-x(1))/(x(3)-x(2))*y(3);

end
%*********************************************
%*********************************************
function yd=FirstDeriv(x,y)

% x and y must be row vectors
% the x entries must be in real-number order and not repeat
% endpoint derivatives are calculated by a forward or 
% backward difference formula involving three points,
% while the middle values are calculated by a central difference formula
% which reduces to the 2 point central difference formula if 
% the x values are equally spaced.
% yd is the approximate derivative at the various values in x.
% The output yd could be used in combination with x and an 
% interpolation method to estimate a derivative function.

if (length(x)~=length(y)) | (sum(size(x))~=length(x)+1) | (sum(size(y))~=length(y)+1) | (length(x)<3)
disp('the two data vectors must be row vectors of equal length, and that length must be at least 3') 

elseif (unique(sort(x))~=x) 
disp('the domain values must be in order and not repeat') 

else

yd(1)=diff3uneq(x(1:3),y(1:3),x(1));
yd(length(x))=diff3uneq(x((length(x)-2):length(x)),y((length(x)-2):length(x)),x(length(x)));

for i=2:(length(x)-1)
yd(i)=diff3uneq(x((i-1):(i+1)),y((i-1):(i+1)),x(i));
end

end
%*********************************************
%*********************************************
function deriv=diff6eq(x,y,t)


% x and y must be rows and have length 6
% t is a number between x(3) and x(4)
% the x entries must be equally spaced
% This script creates the derivative of the data at t using 
% a weighted average producing a fourth-order accurate central difference 
% formula evaluated at x(3) and x(4).

if length(x)~=length(y)
disp('the two data vectors must have equal length') 

% elseif length(x)~=6
% disp('the data vectors must have length 6') 

elseif (sort(x)~=x) | (x(1)==x(2)) | (x(2)==x(3)) | (x(3)==x(4)) | (x(4)==x(5)) | (x(5)==x(6))
disp('the domain values must be in order and not repeat') 

elseif ( t < x(3) ) | ( t>x(4) )
disp('the location of the derivative estimate must be in the interval of the third and fourth data domain values') 

elseif diff(x)-x(2)+x(1)~=zeros(1,length(x)-1)
disp('data vectors must be rows and each data domain increment must be the same length')

elseif (length(t)~=1) | (sum(size(t))~=2)
disp('the third entry must be a number') 

else

h=x(4)-x(3);

deriv=(t-x(3))/h*( -y(6)+8*y(5)-8*y(3)+y(2) )/(12*h)+(x(4)-t)/h*( -y(5)+8*y(4)-8*y(2)+y(1) )/(12*h);



end
%*********************************************
%*********************************************
function [deriv7,deriv5, reldifpercent]=diff7eq(x,y,t)


% x and y must be rows and have length at least 7
% t is a number between x(3) and x(4)
% the x entries must be equally spaced
% This script creates the derivative of the data at t using 
% a weighted average producing a fourth-order accurate central difference 
% formula evaluated at x(3) and x(4).
% deriv5 should also be 4th order accurate and only uses 
% the first 5 points, and can be computed twice as fast. 
% But I am still thinking about it.

if length(x)~=length(y)
disp('the two data vectors must have equal length') 

 elseif length(x)<7
 disp('the data vectors must have length at least 7') 

elseif (sort(x)~=x) | (x(1)==x(2)) | (x(2)==x(3)) | (x(3)==x(4)) | (x(4)==x(5)) | (x(5)==x(6))
disp('the domain values must be in order and not repeat') 

elseif ( t < x(3) ) | ( t>x(4) )
disp('the location of the derivative estimate must be in the interval of the third and fourth data domain values') 

elseif diff(x)-x(2)+x(1)~=zeros(1,length(x)-1)
disp('data vectors must be rows and each data domain increment must be the same length')

elseif (length(t)~=1) | (sum(size(t))~=2)
disp('the third entry must be a number') 

else

h=x(4)-x(3);

deriv7= 4/3*(  0.5*diff3uneq([x(2) x(3) x(4)],[y(2) y(3) y(4)],t)+ 0.5*diff3uneq([x(3) x(4) x(5)],[y(3) y(4) y(5)],t) )-1/3*(  0.5*diff3uneq([x(1) x(3) x(5)],[y(1) y(3) y(5)],t)+ 0.5*diff3uneq([x(3) x(5) x(7)],[y(3) y(5) y(7)],t) );

deriv5=4/3*diff3uneq([x(2) x(3) x(4)],[y(2) y(3) y(4)],t)-1/3*diff3uneq([x(1) x(3) x(5)],[y(1) y(3) y(5)],t);

reldifpercent=100*(deriv5-deriv7)/deriv7;

end
%*********************************************