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