#hw13.txt Cole Franks 08/03/2014

Help13:=proc(): print(): end:

read `C13a.txt`:

###### Problem 1 #####

OtherJesusPi1:=proc(n) local c,s,i,c1,j:
c:=1:
s:=29:
for i from 1 to n do
c:=c*(-1)*product(6*i - j, j=0..5)/((i^6)*(2880^3)):
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:

##### Problem 2 ######

#computes the next RK estimate

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*k1/2}, f):

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

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

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


end:

RK4:=proc(f,x,y,x0,y0,x1,h) local xcur, ycur,L:

xcur:=x0:

ycur:=y0:

while xcur < x1 do 



ycur:=RK41(f, x, y, xcur, ycur, h):

xcur:=xcur + h:


od:

ycur:


end:

####### Problem 3 #####

#I altered euler and impeuler not to make a giant list, 
#because my computer was running out of memory.
#the procedures are EulerN and impEulerN.

EulerN:=proc(f,x,y,x0,y0, x1, h) local xc,yc:


xc:=x0:
yc:=y0:

while xc<x1 do

yc:=Euler1(f,x,y,xc,yc,h):
xc:=xc+h:


od:
yc:

end:


impEulerN:=proc(f,x,y,x0,y0, x1, h) local xc,yc:


xc:=x0:
yc:=y0:

while xc<x1 do

yc:=impEuler1(f,x,y,xc,yc,h):
xc:=xc+h:


od:
yc:

end:

####



EulerH:=proc(f, x, y, x0,y0,x1,Di) local n, g:

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

if g = NULL then 
	return(FAIL):
else g:=rhs(g):
fi:

n:=0:

while evalf(abs(EulerN(f,x,y,x0,y0,x1,(x1 - x0)/2^n) - subs(x=x1,g))) > 10^(-Di) do
n:=n+1:
#print(evalf(abs(EulerN(f,x,y,x0,y0,x1,(x1 - x0)/2^n) - subs(x=x1,g)))):
od:

evalf((x1 - x0)/2^n,10):

end:

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

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

if g = NULL then 
	return(FAIL):
else g:=rhs(g):
fi:

n:=0:

while evalf(abs(impEulerN(f,x,y,x0,y0,x1,(x1 - x0)/2^n) - subs(x=x1,g))) > 10^(-Di) do
n:=n+1:
#print(evalf(abs(impEulerN(f,x,y,x0,y0,x1,(x1 - x0)/2^n) - subs(x=x1,g)))):
od:

evalf((x1 - x0)/2^n,10):

end:

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

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

if g = NULL then 
	return(FAIL):
else g:=rhs(g):
fi:

n:=0:

while evalf(abs(RK4(f,x,y,x0,y0,x1,(x1 - x0)/2^n) - subs(x=x1,g))) > 10^(-Di) do
n:=n+1:
#print(evalf(abs(RK4(f,x,y,x0,y0,x1,(x1 - x0)/2^n) - subs(x=x1,g)))):
od:

evalf((x1 - x0)/2^n,10):

end:
 
#I this code seems perfectly fine, but for some reason Euler1 and EulerN were
#still taking up too much memory and causing the computation to stop.
# I could not figure out why, because the programs have nothing that seems
# to take up lots of memory.