#hw13.txt. March 9, 2014. Katie McKeon.

read(`ExpMath/C13a.txt`):
read(`ExpMath/C13.txt`):

Help13:=proc(): 
print(`OtherJesusPi(n), OtherJesusPi1(n)`):  
print(`RK41(f,x,y,x0,y0,h), RK4(f,x,y,x0,y0,x1,h)`):
print(`EulerH(f,x,y,x0,y0,x1,Di), impEulerH(f,x,y,x0,y0,x1,Di), RK4H(f,x,y,x0,y0,x1,Di)`):
end:


#OtherJesusPi1(n) computes the nth partial sum in Jesus Guillera's approximation
#of 128*sqrt(5)/pi^2
OtherJesusPi1:=proc(n) local c,d,s,i:

c:=1:
d:=1:
s:=29:
for i from 1 to n do
	c:=c*(6*i)*(6*i-1)*(6*i-2)*(6*i-3)*(6*i-4)*(6*i-5)/(i^6):
	d:=(-1)*d/(2880^3):
	s:=s+c*d*(5418*i^2+693*i+29):	
od:

s:

end:

#OtherJesusPi(n) computes the approximation of pi using the partial sum 
#OtherJesusPi1(n) - accurate to ~ 5n digits
OtherJesusPi:=proc(n) 
evalf(sqrt(128*sqrt(5)/OtherJesusPi1(n))):
end:


#RK41(f,x,y,x0,y0,h) performs one iteration of the
#Runge-Kutta method approximating y(x0) for the differential equation
#y'(x) = f(x,y(x)) with y(0)=y0 and interval length h
RK41:=proc(f,x,y,x0,y0,h) local k1, k2, k3, k4:

k1:=subs({x=x0,y=y0},f):
k2:=subs({x=x0+h/2,y=y0+h/2*k1},f):
k3:=subs({x=x0+h/2,y=y0+h/2*k2},f):
k4:=subs({x=x0+h,y=y0+h*k3},f):

y0+h/6*(k1+2*k2+2*k3+k4):

end:


#RK4(f,x,y,x0,y0,x1,h) performs the Runge-Kutta method approximating 
#y(x1) for the differential equation y'(x) = f(x,y(x)) with y(0)=y0 
#and interval length h
RK4:=proc(f,x,y,x0,y0,x1,h) local L,xc,yc:
L:=[y0]:
xc:=x0:
yc:=y0:

while xc<x1 do
	yc:=RK41(f,x,y,xc,yc,h):
	xc:=xc+h:
	L:=[op(L),yc]:
od:

L:

end:



#EulerH(f,x,y,x0,y0,x1,Di) finds the interval length h needed for Euler(f,x,y,x0,y0,x1,h)
#to approximate y(x1) accurate to Di digits
EulerH:=proc(f,x,y,x0,y0,x1,Di) local s,h,d,i,j,L:

s:=dsolve({diff(y(x),x)=subs(y=y(x),f), y(x0)=y0},y(x)):

if s=NULL then
	RETURN(FAIL):
else
	s:=op(2,s):
	d:=x1-x0:
	i:=1:
	do
		h:=d/i:
		L:=Euler(f,x,y,x0,y0,x1,h):
		if evalf(abs(subs(x=x1,s)-L[i+1]))< evalf(10^(-Di)) then
			RETURN(h):
		fi:
		i:=i+1:
	od:
fi:

end:



#impEulerH(f,x,y,x0,y0,x1,Di) finds the interval length h needed for impEuler(f,x,y,x0,y0,x1,h)
#to approximate y(x1) accurate to Di digits
impEulerH:=proc(f,x,y,x0,y0,x1,Di) local s,h,d,i,j,L:

s:=dsolve({diff(y(x),x)=subs(y=y(x),f), y(x0)=y0},y(x)):

if s=NULL then
	RETURN(FAIL):
else
	s:=op(2,s):
	d:=x1-x0:
	i:=1:
	do
		h:=d/i:
		L:=impEuler(f,x,y,x0,y0,x1,h):
		if evalf(abs(subs(x=x1,s)-L[i+1]))< evalf(10^(-Di)) then
			RETURN(h):
		fi:
		i:=i+1:
	od:
fi:

end:


#RK4H(f,x,y,x0,y0,x1,Di) finds the interval length h needed for RK4(f,x,y,x0,y0,x1,h)
#to approximate y(x1) accurate to Di digits
RK4H:=proc(f,x,y,x0,y0,x1,Di) local s,h,d,i,j,L:

s:=dsolve({diff(y(x),x)=subs(y=y(x),f), y(x0)=y0},y(x)):

if s=NULL then
	RETURN(FAIL):
else
	s:=op(2,s):
	d:=x1-x0:
	i:=1:
	do
		h:=d/i:
		L:=RK4(f,x,y,x0,y0,x1,h):
		if evalf(abs(subs(x=x1,s)-L[i+1]))< evalf(10^(-Di)) then
			RETURN(h):
		fi:
		i:=i+1:
	od:
fi:

end:


#DATA - In all cases x0=0, y0=1, x1=1, Di=10
#	#  EulerH	#  impEulerH	#  RK4H
###################################################
# y	#	?	#	?	# 1/123
# 1+x+y	#	?	#	?	# 1/162
# 1+x*y	#	?	#	?	# 1/96
#
#Computation for EulerH and impEulerH had to be interrupted in all cases
#Here are the results when Di=5
#	#  EulerH	#  impEulerH	#  RK4H
###################################################
# y	#	?	#  1/213	# 1/7
# 1+x+y	#	?	#  1/369	# 1/9
# 1+x*y	#	?	#  1/67		# 1/6
#for all three EulerH cases, the interval length is smaller than 1/1000