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

Help := proc();

print(` AdaptiveQuad, Euler`);
print(` Taylor2 `);
end:

# Section 4.6 
# Adaptive Quadrature method
# The code here follow the same method
# we discussed during the class with
# the safety factor 10 insted of the usual 15.

# Input: function f, the initial point a,
# end point b, tolerance TOL;
#
# Output: the approximation of int(f) 
# from a to b
# with error bound less than TOL.

# Example: problem 3a)
# Check the answer in the back to the book
# Try: 
# AdaptiveQuad(x-> exp(2*x)*sin(3*x),1,3,10^(-5));

# Another example: problem 3b)
# Try:
# AdaptiveQuad(x-> exp(3*x)*sin(2*x),1,3,10^(-5));


AdaptiveQuad := proc(f,a,b,TOL)  local h1, S1, h2, S2; 

# Simpson's rule approximation from a to b
h1 := (b-a)/2;
S1 := h1/3*(f(a)+4*f(a+h1)+f(a+2*h1));

# Simpson's rule approximation 
# from a to (a+b)/2 and (a+b)/2 to b
h2 := (b-a)/4;
S2 := h2/3*(f(a)+4*f(a+h2)+f(a+2*h2))
	 + h2/3*(f(a+2*h2)+4*f(a+3*h2)+f(a+4*h2));

if evalf(abs(S1-S2)/10) < TOL then
	return(evalf(S2));
else
	return (AdaptiveQuad(f, a,(a+b)/2, TOL/2) 
		  + AdaptiveQuad(f,(a+b)/2, b ,TOL/2));
fi:

end:


# Section 5.2 
# Euler's method
# Follow algorithm in the book page 257

# Input: function f, the initial point a,
# end point b, step size h , initial value alpha;
#
# Output: print the value of an approximation at 
# each t_i, t_i = a+ih

# Example: problem 1a)
# Check the answer in the back to the book
# Try: 
# Euler((t,y)-> t*exp(3*t)-2*y,0,1,0.5,0);

# Another example: problem 5a)
# Try:
# Euler((t,y)-> y/t-(y/t)^2,1,2,0.1,1);


Euler := proc(f, a, b, h, alpha) 

local i, t, w;

t := a;
w := alpha;
print(t,w);

for i from 1 to (b-a)/h do

	w := evalf(w+h*f(t,w));
	t := evalf(a+i*h);
	print(t,w);

od:

return();

end:

# Section 5.3 
# Taylor's method of order 2

# Note: you have to calculate fprime by hand.
# Then use it as part of the input.

# Input: function f, fprime, the initial point a,
# end point b, step size h , initial value alpha;
#
# Output: print the value of an approximation at 
# each t_i, t_i = a+ih

# Example: problem 1a)
# Check the answer in the back to the book
# Try: 
# Taylor2((t,y)-> t*exp(3*t)-2*y,
# (t,y)-> (1+t)*exp(3*t)+4*y,  
#  0, 1, 0.5, 0);

Taylor2 := proc(f, fprime, a, b, h, alpha) local i, t, w;

t := a;
w := alpha;
print(t,w);

for i from 1 to (b-a)/h do

	w := evalf(w+h*f(t,w)+(h^2)/2*fprime(t,w));
	t := evalf(a+i*h);
	print(t,w);

od:

return();
end: