#HW13.txt  Ross Berkowitz
#read("E:\\Documents\\Maple\\em14\\hw13.txt");

Help:=proc():
print(``):
end:
################### 1   ##################

OtherJesusPi:=proc(N) local i,pi,n,term:
Digits:=max(10,5*N+floor(sqrt(N))):
term:=1.:
pi:=29.00:
for n from 1 to N do
	term:=term*(-1)*mul(i,i=6*(n-1)+1..6*n)/((2880)^3*n^6):
	pi:=pi+term*(5418*n^2+693*n+29):
od:
evalf(sqrt(128*sqrt(5)/pi)):
end:


###################  2   ###################

RK4:=proc(f,x,y,x0,y0,x1,h) local N,k1,k2,k3,k4,n,L:
N:=floor((x1-x0)/h):
L:=[y0]:
for n from 1 to N do
	k1:=subs({x=x0+(n-1)*h,y=L[n]}, f):
	k2:=subs({x=x0+(n-1+1/2)*h,y=L[n]+1/2*h*k1},f):
	k3:=subs({x=x0+(n-1+1/2)*h,y=L[n]+1/2*h*k2},f):
	k4:=subs({x=x0+n*h,y=L[n]+h*k3},f):
	
	L:=[op(L),L[n]+h/6*(k1+2*k2+2*k3+k4)]:
od:
L:
end:



##########################  3  #######################
EulerH:=proc(f,x,y,x0,y0,x1,Di) local n,h,sol:
Digits:=Di+10:

sol:=dsolve({diff(y(x),x)=subs(y=y(x),f),y(0)=y0}):
if sol=null then
	RETURN(FAIL):
fi:
sol:=rhs(sol):
h:=x1-x0:
n:=2:
while abs(evalf(Euler(f,x,y,x0,y0,x1,h)[n]-subs(x=x1,sol)))>10^(-Di) do
	h:=h/2:
	n:=2*n-1:
od:
h:
end:

RK4H:=proc(f,x,y,x0,y0,x1,Di) local n,h,sol:
Digits:=Di+10:

sol:=dsolve({diff(y(x),x)=subs(y=y(x),f),y(0)=y0}):
if sol=null then
	RETURN(FAIL):
fi:
sol:=rhs(sol):
h:=x1-x0:
n:=2:
while abs(evalf(RK4(f,x,y,x0,y0,x1,h)[n]-subs(x=x1,sol)))>10^(-Di) do
	h:=h/2:
	n:=2*n-1:
od:
h:
end:

impEulerH:=proc(f,x,y,x0,y0,x1,Di) local n,h,sol:
Digits:=Di+10:

sol:=dsolve({diff(y(x),x)=subs(y=y(x),f),y(0)=y0}):
if sol=null then
	RETURN(FAIL):
fi:
sol:=rhs(sol):
h:=x1-x0:
n:=2:
while abs(evalf(impEuler(f,x,y,x0,y0,x1,h)[n]-subs(x=x1,sol)))>10^(-Di) do
	h:=h/2:
	n:=2*n-1:
od:
h:
end:


##########   CLASS STUFFF  ###############

#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: