#hw13.txt; March 6, 2014; Frank Wagner
Help:=proc(): print(` OtherJesus(n)  `): print(` RK41(f,x,y,x0,y0,h) `): print(` RK4(f,x,y,x0,y0,x1,h) `):
print(` EulerH(f,x,y,x0,y0,x1,Di) `): print(` impEuler(f,x,y,x0,y0,x1,Di) `): 
print(`RK4H(f,x,y,x0,y0,x1,Di) `):  end:


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

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

c:=1:
s:=29:

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

evalf(sqrt(128*sqrt(5)/s)):

end:



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

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/6*(k1+2*k2+2*k3+k4):

end:


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:




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

#old stuff needed for this problem#

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


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:=proc(f,x,y,x0,y0,h) local ystar:

ystar:=Euler1(f,x,y,x0,y0,h):

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

end:


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



EulerH:=proc(f,x,y,x0,y0,x1,Di) local F,A,d,h,sol:

F:=dsolve({diff(y(x),x)=subs(y=y(x),f),y(x0)=y0}):
sol:=eval(y(x),F):
A:=evalf(evalf(subs(x=x1,sol),Di),Di*10):

h:=x1-x0:
d:=evalf(abs(A-Euler(f,x,y,x0,y0,x1,h)[nops(Euler(f,x,y,x0,y0,x1,h))]),Di*10):

 while d>=10^(-Di) do
  h:=h/10:
  d:=evalf(abs(A-Euler(f,x,y,x0,y0,x1,h)[nops(Euler(f,x,y,x0,y0,x1,h))]),Di*10):
 od:


 

h:
end:


impEulerH:=proc(f,x,y,x0,y0,x1,Di) local F,A,d,h,sol:

F:=dsolve({diff(y(x),x)=subs(y=y(x),f),y(x0)=y0}):
sol:=eval(y(x),F):
A:=evalf(subs(x=x1,sol),Di*10):

h:=x1-x0:
d:=evalf(abs(A-impEuler(f,x,y,x0,y0,x1,h)[nops(impEuler(f,x,y,x0,y0,x1,h))]),Di*10):

 while d>=10^(-Di) do
  h:=h/10:
  d:=evalf(abs(A-impEuler(f,x,y,x0,y0,x1,h)[nops(impEuler(f,x,y,x0,y0,x1,h))]),Di*10):
 od:
 

h:
end:


RK4H:=proc(f,x,y,x0,y0,x1,Di) local F,A,d,h,sol:

F:=dsolve({diff(y(x),x)=subs(y=y(x),f),y(x0)=y0}):
sol:=eval(y(x),F):
A:=evalf(subs(x=x1,sol),Di*10):

h:=evalf(x1-x0):
d:=evalf(abs(A-RK4(f,x,y,x0,y0,x1,h)[nops(RK4(f,x,y,x0,y0,x1,h))]),Di*10):

 while d>=10^(-Di) do
  h:=h/10:
  d:=evalf(abs(A-RK4(f,x,y,x0,y0,x1,h)[nops(RK4(f,x,y,x0,y0,x1,h))]),Di*10):
 od:
 

h:
end:


#My original attempt at these programs included something analogous to the bisection method in order
#to find the correct value of h. However, that proved to be too time-consuming.
#Even with this simplistic method which only returns the greatest value of h that is a power of 10,
#EulerH and impEulerH take far too long, while RK4H is quick enough that I suppose it could make some
#improvements.
#In fact, all EulerH's and impEulerH's took longer than I was willing to wait.
#What could we change to make this quicker?
#
#
#RK4H(y,x,y,0,1,x1,10) returned .001
#
#RK4H(1+x+y,x,y,0,1,x1,10) returned .001
#
#RK4H(1+x*y,x,y,0,1,x1,10) returned .01