#!/usr/local/bin/maple # -*- maplev -*- # Nathaniel Shar # HW 13 - 9 Mar 2014 # Experimental Mathematics ############# # Problem 1 # ############# f1 := proc(n) option remember: if n = 0 then: return 1: else: return -pochhammer(6*n-5, 6)/(2880^3*n^6) * f1(n-1): end: end: OtherJesusPi := proc(n) option remember: return sqrt(128*sqrt(5)/add(f1(i)*(5418*i^2 + 693*i + 29), i=0..n)): end: # It actually seems closer to 6 digits per term. ############# # Problem 2 # ############# Euler := proc(f, x0, y0, mesh, x1) local y1, i: y1 := y0: for i from x0 to x1 by mesh while i < x1 do: y1 := y1 + mesh*f(i, y1): od: return y1: end: ImprovedEuler := proc(f, x0, y0, mesh, x1) local y1, i, ystar: option remember: y1 := y0: for i from x0 by mesh while i < x1 do: ystar := y1 + mesh*f(i, y1): y1 := y1 + mesh*(f(i, y1) + f(i+mesh, ystar))/2: od: return y1: end: RK4 := proc(f, x0, y0, mesh, x1) local y1, i, k: option remember: y1 := y0: for i from x0 by mesh while i < x1 do: k[1] := f(i, y1): k[2] := f(i + mesh/2, y1 + mesh*k[1]/2): k[3] := f(i + mesh/2, y1 + mesh*k[2]/2): k[4] := f(i + mesh, y1 + mesh*k[3]): y1 := y1 + mesh/6*(k[1] + 2*k[2] + 2*k[3] + k[4]): od: return y1: end: ############# # Problem 3 # ############# FindMeshSize := proc(ApproxMethod, f, x0, y0, x1, digits) local x, y, soln, val, MAXSTEP: MAXSTEP := 18: Digits := max(10, 2*digits): soln := dsolve({diff(y(x), x) = f(x, y(x)), y(x0) = y0}): if soln = NULL then: return FAIL: fi: y := t->subs(x=t, op(2, soln)): val := y(x1): return 2^(-bsearchForMax([seq(MAXSTEP-i, i=0..MAXSTEP)], proc(t): return evalf(log[10](abs(ApproxMethod(f, x0, y0, evalf((x1-x0)/2^t), x1) - val))) < -digits: end proc )): end: # Assumes list is sorted and predicate is nondecreasing. bsearchForMax := proc(L, predicate) local n, s: n := nops(L): if predicate(L[n]) then: return L[n]: elif n <= 1 then: return L[1]: else: s := floor((n-1)/2)+1: if predicate(L[s]) then: return bsearchForMax(L[s+1..n], predicate): else: return bsearchForMax(L[1..s], predicate): fi: fi: end: # Results: #1: # Euler: 10 digits? funny joke # ImprovedEuler: 2^-16 # Runge-Kutta: 2^-6 #2: # Euler: 10 digits? A monkey at a typewriter will take less time to # type the correct 10 digits. # ImprovedEuler: 2^-16 # Runge-Kutta: 2^-7 #3: # Euler: 10 digits? And Congress did something good for the country # too, right? # ImprovedEuler: 2^-14 # RungeKutta: 2^-6