#!/usr/local/bin/maple
# -*- maplev -*-

# Nathaniel Shar
# HW 11 - 2 Mar 2014
# Experimental Mathematics

# From c11.txt:

#C11.txt: Feb. 27, 2014
Help:=proc(): print(` Orb(g,x0,K1,K2,d1)`):
print(` OrbF(f,x,g,x0,K1,K2,d1) `):
end:

#Orb(g,x0,K1,K2,d1): running the logistic map
#x->g*x*(1-x) K1+K2 times, starting at x0,  and displaying
#the last K2 terms
#For example:
#Orb(2.5,.5,10000,1000);
Orb:=proc(g,x0,K1,K2,d1) local L,i,x1:
    x1:=x0:
    L:=[]:

    for i from 1 to K1+K2+1 do
        L:=[op(L),x1]:
        x1:=g*x1*(1-x1):
    od:

    evalf(L[K1+1..K1+K2],d1):

end:


#OrbF(f,x,g,x0,K1,K2,d1): running the general logistic map
#x->g*f K1+K2 times and displaying
#the last K2 terms
#For example:
#OrbF(2.5,.5,10000,1000);
OrbF:=proc(f,g,x0,K1,K2,d1) local L,i,x1:
    x1:=x0:
    L:=[]:

    for i from 1 to K1+K2+1 do
        L:=[op(L),x1]:
        x1:=g*f(x1):
    od:

    map(x->truncate(x, d1), L[K1+1..K1+K2]):

end:

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

# So, this problem is hard. This isn't a perfect solution -- it can
# really screw up and give nonsensical results if Digits isn't set
# high enough, and it's possible that the intervals will be wider than
# epsilon if some strange transitional period occurs. And it isn't
# very smart about recognizing the pattern of periods, so it can get
# confused sometimes by sporadic low-period trajectories.

# The good news is it does seem to work most of the time anyway,
# especially if you set maxmul to a low value. That makes it take
# longer, so you have to be patient.

truncate := proc(val, digits):
    return floor(val*10^digits)/10^digits:
end:

OrbitSize := proc(f, g, x0, K1, K2, d1) option remember:
    return nops(convert(OrbF(f, g, x0, K1, K2, d1), set)):
end:

# Note: maxmul should be a power of 2 if you want the intervals to
# have lengths that are integer multiples of epsilon.
NumerBifOrbF := proc(f, x0, startup, K, precision, n, epsilon, amax := 10, maxmul := 64) local i, a, Stack, T, lasta, pa, lower, upper, rv:
    T := table();
    a := 0:
    Stack := stack[new](maxmul * epsilon);
    while not stack[empty](Stack) do:
        lasta := a:
        a := stack[pop](Stack):
        pa := OrbitSize(f, a, x0, startup, K, precision):
        print(a, pa, lasta, T[lasta]):
        T[a] := pa:
        # If we got a change, go back and search at a finer mesh size.
        if T[a] <> T[lasta] and a-lasta > epsilon then:
            stack[push](a, Stack):
            stack[push](a-(a-lasta)/2, Stack):
            stack[push](lasta, Stack):
            a := lasta:
        else:
            if stack[empty](Stack) and a <= amax then:
                if T[a] = K then:
                    printf("Warning: Stopping at %f due to apparent chaos\n", a):
                    # Don't go on -- further data is probably
                    # useless.
                else:

                    stack[push](a+maxmul*epsilon, Stack):
                fi:
            fi:
        fi:
    od:
    rv := []:
    for i from 0 to n-1 do:
        lower := GetMax(select(x->evalb(T[x]=2^i), [indices(T, nolist)]), x->x)[1]:
        upper := GetMin(select(x->evalb(T[x]=2^(i+1)), [indices(T, nolist)]), x->x)[1]:
        rv := [op(rv), [lower, upper]]:
    od:
end:

GetMax := proc(S, f) local champion, record, new, newval, i:
    if S = {} then:
        return FAIL:
    fi:
    champion := S[1]:
    record := f(champion):
    for i from 2 to nops(S) do:
        new := S[i]:
        newval := f(new):
        if newval > record then:
            champion := new:
            record := newval:
        fi:
    od:
    return [champion, record]:
end:

GetMin := proc(S, f) local champion, record, new, newval, i:
    if S = {} then:
        return FAIL:
    fi:
    champion := S[1]:
    record := f(champion):
    for i from 2 to nops(S) do:
        new := S[i]:
        newval := f(new):
        if newval < record then:
            champion := new:
            record := newval:
        fi:
    od:
    return [champion, record]:
end:

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

# These values are sensitive to Digits. I did them with Digits = 20
# but I am not sure if this was enough.

# NumerBifOrbF(x->x^2*(1-x),.5,10000,1000,10,4,.005, 7, 256)
# = [[5.33, 5.335], [5.76, 5.765], [5.86, 5.865], [5.88, 5.885]]

# NumerBifOrbF(x->x*(1-x)^2,.5,10000,1000,10,4,.005, 7, 256)
# = [[4.00, 4.005], [4.995, 5.005], [5.23, 5.24], [5.285, 5.290]]

# NumerBifOrbF(x->x^2*(1-x)^2,.5,10000,1000,10,4,.02, 14)
# = [[10.40, 10.42], [11.80, 11.84], [12.14, 12.16], [12.22, 12.24]]

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

# I don't know how to do this problem, but AllCV at least gives some
# candidates (16/3 for the first one and 4 for the second one) that
# seem to work. It blows up horribly for the 3rd one and doesn't give
# explicit results at all. Even fsolve didn't work, but looking at a
# plot suggests that the first bifurcation should be between 11.6 and
# 11.7. Obviously that doesn't match with the data. So I don't know if
# this is the right approach at all.

FixedPoints := proc(f) local x:
    return [solve(f(x)=x, x)]:
end:

CriticalValues := proc(f, pt, param) local x:
    return [solve(abs(subs({x=pt}, diff(f(x), x))) = 1, param)]:
end:

AllCV := proc(f, param):
    return map(x->CriticalValues(f, x, param), FixedPoints(f)):
end:

#############
# Problem 4 #
#############

InMandel := proc(c, K) local t, i:
    t := c:
    for i from 1 to K do:
        t := t^2 + c:
        # These "10"s could be set a lot lower. I think 2 would
        # work. I'm using a higher value for safety.
        if Re(t) > 10 or Im(t) > 10 then:
            return false:
        fi:
    od:
    return true:
end:

Mandel := proc(res, size, K) local i, j, L:
    L := []:
    for i from -size to size by res do:
        for j from -size to size by res do:
            print(i, j):
            if InMandel(i + I*j, K) then:
                L := [op(L), [i, j]]:
            fi:
        od:
    od:
    return L:
end:

# MS.jpg was produced with res = 0.01, K = 30.