#!/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