#Nathan Fox #Homework 13 #I give permission for this work to be posted online #Read procedures from class read(`C12.txt`): read(`C13.txt`): read(`C13a.txt`): Help:=proc(): print(` g(n) , OtherJesusPi(N) , RK41(f, x, y, x0, y0, h)`): print(` RK4(f, x, y, x0, y0, x1, h) `): print(` MethodH(Method,f,x,y,x0,y0,x1,Di) `): print(` EulerH(f,x,y,x0,y0,x1,Di) `): print(` impEulerH(f,x,y,x0,y0,x1,Di) `): print(` RK4H(f,x,y,x0,y0,x1,Di) `): print(` Prob3(f,x,y,x0,y0,x1,Di) `): end: ##PROBLEM 1## #coefficient ratio in JG formula g:=proc(n): -(6*n-1)*(6*n-2)*(6*n-3)*(6*n-4)*(6*n-5)/(3981312000*n^5): end: #Optimized JG formula OtherJesusPi:=proc(N) local ret, tot, n: ret:=29: tot:=1: for n from 1 to N do tot:=tot * g(n): ret:=ret + tot * (5418*n^2+693*n+29): od: sqrt(128*sqrt(5)/ret): end: ##PROBLEM 2## #RK41(f, x, y, x0, y0, 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*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): return y0+(h/6)*(k1+2*k2+2*k3+k4): end: #RK4(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:=RK41(f, x, y, xc, yc, h): xc:=xc+h: L:=[op(L), yc]: od: return L: end: ##PROBLEM 3## #MethodH(Method,f,x,y,x0,y0,x1,Di): #inputs an expression f in x, y such #that Maple's dsolve knows to find the exact solution to the #initial value problem #y'(x)=f(x,y(x)) , y(x0)=y0 , #(if it can't the procedure returns FAIL), and then finds out the #largest mesh size h for which Method gives an approximation #correct up to Di digits. #The h that is returned is always a dyadic fraction of abs(x1-x0) #Method can be a list of methods, in which case it returns a list #of methods MethodH:=proc(Method,f,x,y,x0,y0,x1,Di) local soln, H, M, Meth: soln:=DS(f, x, y, x0, y0): if soln = FAIL then print(`Maple isn't THAT smart!`): return FAIL: fi: soln:=evalf(subs(x=x1, soln), Di+1): if type(Method, list) then M:=Method: else M:=[Method]: fi: H:=[]: for Meth in Method do H:=[op(H), abs(x1-x0)]: while abs(Meth(f, x, y, x0, y0, x1, H[-1])[-1] - soln) > 10^(-Di) do H[-1]:=H[-1]/2: od: od: if type(Method, list) then return H: else return H[1]: fi: end: #EulerH(f,x,y,x0,y0,x1,Di): MethodH with Euler EulerH:=proc(f,x,y,x0,y0,x1,Di) return MethodH(Euler, f, x, y, x0, y0, x1, Di): end: #impEulerH(f,x,y,x0,y0,x1,Di): MethodH with impEuler impEulerH:=proc(f,x,y,x0,y0,x1,Di) return MethodH(impEuler, f, x, y, x0, y0, x1, Di): end: #RK4H(f,x,y,x0,y0,x1,Di): MethodH with RK4 RK4H:=proc(f,x,y,x0,y0,x1,Di) return MethodH(RK4, f, x, y, x0, y0, x1, Di): end: #Prob3(f,x,y,x0,y0,x1,Di): MethodH with [Euler, impEuler, RK4] Prob3:=proc(f,x,y,x0,y0,x1,Di) return MethodH([Euler, impEuler, RK4], f, x, y, x0, y0, x1, Di): end: #RK4H(y,x,y,0,1,x1,10) returned 1/128 #The other two took a really long time #impEulerH(y, x, y, 0, 1, 1, 5) returned 1/256 #EulerH(y, x, y, 0, 1, 1, 2) returned 1/256 #RK4H(1+x+y,x,y,0,1,1,10) took a really long time #RK4H(1+x+y,x,y,0,1,1,9) returned 1/128 #impEulerH(1+x+y,x,y,0,1,1,5) returned 1/512 #EulerH(1+x+y,x,y,0,1,1,3) returned 1/4096 #RK4H(1+x*y,x,y,0,1,1,10) took a really long time #RK4H(1+x*y,x,y,0,1,1,9) returned 1/64 #impEulerH(1+x*y,x,y,0,1,1,5) returned 1/256 #EulerH(1+x*y,x,y,0,1,1,2) returned 1/128