#restart; read `W:/hw4.txt`;

Help := proc();

print(`RichExtra, RichTable`);
print(`Trapezoidal, Simpson, Romberg`);
print(`NewtonCotes1`);
end:

# Section 4.2 
# Richard's Extrapolation

# Input: function f, the initial point start,
# step size h, number of iteration N0;
#
# Output: the approximation of f'(start)
# (using three-point formula): N_N0(start)

# Example: problem 1a)
# Check the answer in the back to the book
# Try: 
# RichExtra(x->ln(x),1,0.4,3);

# Another example: problem 1b)
# Try:
# RichExtra(x->x+exp(x),0,0.4,3);

RichExtra := proc(f, start, h, N0) local A;

A := [seq((f(start+h/(2^i))-f(start-h/(2^i)))/
     (2*h/(2^i)), i = 0..N0-1)];

RichTable(A)[N0][N0];
end:

# Section 4.2 
# Richard's Extrapolation

# Input: the list of value 
# N(h), N(h/2), N(h/4), ...
# Output: the table of N_i(h)

# Example: Table 4.6 page 183 in the book
# Try: 
# RichTable([22.414160,22.228786,22.182564]); 

# Another example: problem 5)
# Check the answer in the back to the book
# Try: 
# RichTable([1.57096, 1.896119, 1.974232, 1.993570]);

#############################
#
# Note: the formula is form 2 
# under assumption
# M = N(h) + K_1h^2 + K_2h^4  
#     + K_3h^6  + ... 
#
#############################

RichTable := proc(A) local i,j,T,n;

n := nops(A);

T := [seq([A[i],0$(n-1)],i=1..n)];

for i from 2 to n do
	for j from 2 to i do
		T[i][j] := T[i][j-1]+
		(T[i][j-1]-T[i-1][j-1])/(4^(j-1)-1);
	od:
od:

return(T);
end:

# Section 4.3, 4.4 
# Trapezoidal Rule

# Input: function f, the initial point a,
# end point b, number of interval k;
#
# Output: approximation of intergration
# integrate f from a to b with k intervals
# using Trapezoidal rule

# Example: section 4.3 problem 1a) 
# Check the answer in the back to the book
# Try: 
# Trapezoidal(x->x^4 , 1/2, 1, 1);

# Another example: section 4.4 problem 1a) 
# Check the answer in the back to the book
# Try: 
# Trapezoidal(x->x*ln(x), 1, 2, 4);

Trapezoidal := proc(f,a,b,k) local h;

h := (b-a)/k;

evalf(sum( (f(a+h*(i-1))+f(a+h*i))
*h/2, i = 1..k));

end:

# Section 4.3, 4.4 
# Simpson Rule

# Input: function f, the initial point a,
# end point b, number of interval k;
#
# Output: approximation of intergration
# integrate f from a to b with k intervals
# using Simpson rule

# Example: section 4.3 problem 5a) 
# Check the answer in the back to the book
# Try: 
# Simpson(x->x^4 ,1/2, 1, 2);

# Another example: section 4.4 problem 3a) 
# Check the answer in the back to the book
# Try: 
# Simpson(x -> x*ln(x), 1, 2, 4);

Simpson := proc(f,a,b,k) local h;

h := (b-a)/k;

evalf(sum((f(a+h*(2*i-2))+
4*f(a+h*(2*i-1))+f(a+h*(2*i)))
*h/3, i = 1..k/2));

end:

# Section 4.5 
# Romberg Integration

# Input: function f, the initial point a,
# end point b, number of iteration N0;
#
# Output: approximation of intergration
# R_N0_N0

# Example: section 4.5 problem 1a)
# Check the answer in the back to the book
# Try: 
# Romberg(x->x^2*ln(x),1,1.5,3);

# Another example: problem 1b)
# Try:
# Romberg(x->x^2*exp(-x),0,1,3);

Romberg := proc(f,a,b,N0) local i,A;

A := [seq(Trapezoidal(f,a,b,2^(i-1)),i=1..N0)];

RichTable(A)[N0][N0];

end: