#hw13 #Bo Zhen #em14 #3/9/2014


################################content of C13.txt and C13a.txt###################################
#C13.txt; March 6, 2014; 
#Implementing Jesus' Guillera PROVED formula for Pi

Help0:=proc(): print(` JesusPi(n), JesusPi1(n), JesusPi1e(n) `):  end:

#JesusPi1(n): computes the n-th partial sum of
#https://sites.google.com/site/guilleramath/
#due to Jesus Guillera. It is PROVED, and
#much more importantly using the so-called WZ
#method (that stands for Wilf and Zeilberger)
JesusPi1:=proc(n) local c,s,i,c1:
c:=1:
s:=13:
for i from 1 to n do
c:=c*(-1)/2^10*((2*i-1)/(2*i))^5:
c1:=c*((20*i)*(41*i+9)+13):
s:=s+c1:
od:

s:

end:

#JesusPi(n): computes Pi from 128/Pi^2
JesusPi:=proc(n): evalf(sqrt(128/JesusPi1(n))):end:


#floating point version
JesusPi1e:=proc(n) local c,s,i,c1:
c:=1:
s:=13:
for i from 1 to n do
c:=evalf(c*(-1)/2^10*((2*i-1)/(2*i))^5):
c1:=evalf(c*((20*i)*(41*i+9)+13)):
s:=evalf(s+c1):
od:

s:

end:



#March 6, 2014; Starting Numerical solutions of ODEs

Help:=proc(): print(`Euler1(f,x,y,x0,y0,h),  `): 
print(`Euler(f,x,y,x0,y0,x1,h)`):
print(`impEuler1(f,x,y,x0,y0,h) `):
print(`impEuler(f,x,y,x0,y0,x1,h) `):
end:

#Euler1(f,x,y,x0,y0,h);
Euler1:=proc(f,x,y,x0,y0,h) 
y0+h*subs({x=x0,y=y0},f):
end:

#Euler(f,x,y,x0,y0,x1,h): Euler's method approximating
#y'(x)=f(x,y), y(x0)=y0, all the way to x1
#with increments of h
Euler:=proc(f,x,y,x0,y0, x1, h) local L,xc,yc:
L:=[y0]:

xc:=x0:
yc:=y0:

while xc<x1 do

yc:=Euler1(f,x,y,xc,yc,h):
xc:=xc+h:
L:=[op(L),yc]:

od:
L:

end:

#impEuler1(f,x,y,x0,y0,h);
impEuler1:=proc(f,x,y,x0,y0,h) local ystar:

ystar:=y0+h*subs({x=x0,y=y0},f):
y0+h*(subs({x=x0,y=y0},f)+subs({x=x0+h,y=ystar},f))/2:


end:

#impEuler(f,x,y,x0,y0,x1,h): Euler's method approximating
#y'(x)=f(x,y), y(x0)=y0, all the way to x1
#with increments of h
impEuler:=proc(f,x,y,x0,y0, x1, h) local L,xc,yc:
L:=[y0]:

xc:=x0:
yc:=y0:

while xc<x1 do

yc:=impEuler1(f,x,y,xc,yc,h):
xc:=xc+h:
L:=[op(L),yc]:

od:
L:

end:

#################################################end old content########################################


#1

OtherJesusPi1:=proc(n) local c,s,i,c1:

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

s:

end:

OtherJesusPi:=proc(n): 
evalf(sqrt(128*sqrt(5)/OtherJesusPi1(n))):
end:


#2


#RK1(f,x,y,x0,y0,h);
RK1:=proc(f,x,y,x0,y0,h) local k1,k2,k3,k4:

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

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

end:

#Euler(f,x,y,x0,y0,x1,h): Euler's method approximating
#y'(x)=f(x,y), y(x0)=y0, all the way to x1
#with increments of 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:=RK1(f,x,y,xc,yc,h):
xc:=xc+h:
L:=[op(L),yc]:

od:
L:

end:


#3
#the codes here run for about Di=5 and takes incredibly long for anything greater than that.
#not sure how to resolve it.  


EulerH:=proc(f,x,y,x0,y0,x1,Di) local eps,aList,a,b,h,Y:

Digits:=100:

h:=.1:
eps:=(10^(-Di)):
aList:=Euler(f,x,y,x0,y0,x1,h):
a:=aList[nops(aList)]:
b:=subs(x=x1,rhs(dsolve({diff(y(x),x)=subs(y=y(x),f),y(x0)=y0}))):

while evalf(abs(a-b))>eps do
	h:=.5*h:
	print(h):
	aList:=Euler(f,x,y,x0,y0,x1,h):
	a:=aList[nops(aList)]:
	b:=subs(x=x1,rhs(dsolve({diff(y(x),x)=subs(y=y(x),f),y(x0)=y0}))):
	print(evalf(a-b)):
od:  

h:
	

end:

impEulerH:=proc(f,x,y,x0,y0,x1,Di) local eps,aList,a,b,h,Y:

Digits:=100:

h:=.1:
eps:=(10^(-Di)):
aList:=impEuler(f,x,y,x0,y0,x1,h):
a:=aList[nops(aList)]:
b:=subs(x=x1,rhs(dsolve({diff(y(x),x)=subs(y=y(x),f),y(x0)=y0}))):

while evalf(abs(a-b))>eps do
	h:=.5*h:
	print(h):
	aList:=impEuler(f,x,y,x0,y0,x1,h):
	a:=aList[nops(aList)]:
	b:=subs(x=x1,rhs(dsolve({diff(y(x),x)=subs(y=y(x),f),y(x0)=y0}))):
	print(evalf(a-b)):
od:  

h:
	
end:

RK4H:=proc(f,x,y,x0,y0,x1,Di) local eps,aList,a,b,h,Y:

Digits:=100:

h:=.1:
eps:=(10^(-Di)):
aList:=RK4H(f,x,y,x0,y0,x1,h):
a:=aList[nops(aList)]:
b:=subs(x=x1,rhs(dsolve({diff(y(x),x)=subs(y=y(x),f),y(x0)=y0}))):

while evalf(abs(a-b))>eps do
	h:=.5*h:
	print(h):
	aList:=RK4H(f,x,y,x0,y0,x1,h):
	a:=aList[nops(aList)]:
	b:=subs(x=x1,rhs(dsolve({diff(y(x),x)=subs(y=y(x),f),y(x0)=y0}))):
	print(evalf(a-b)):
od:  

h:
	

end:

impEulerH:=proc(f,x,y,x0,y0,x1,Di) local eps:

end:

RK4H:=proc(f,x,y,x0,y0,x1,Di) local eps:

end: