#hw17.txt; Due March 30, 2014; Frank Wagner
Help:=proc(): print(` SIrec(L,Ini,f,m,n)  `): print(` dIBVP(L,Ini,Fini,f,n,N) `):
print(` GR(p,N) `): print(` LGR(p,N) `): end:


###################
#####Problem 1#####
###################
#SIrec(L,Ini,f,m,n) inputs lists L and Ini, a function
#f in m, and a number n and outputs the n-th term of the
#solution of the linear recurrence equation with
#constant coefficients
#a(n)=L[1]*a(n-1)+L[2]*a(n-2)+ ... + L[-1]*a(n-nops(L))+f(n)
#and initial conditions
#a(0)=Ini[1], ..., a(nops(Ini)-1)=Ini[-1]

SIrec:=proc(L,Ini,f,m,n) local i:

if nops(L)<>nops(Ini) then
print(Ini, L, `should be lists of the SAME length `):
 RETURN(FAIL):
fi:

if n<0 then
 0:
elif n<nops(Ini) then
 Ini[n+1]:
else
 add(L[i]*SIrec(L,Ini,f,m,n-i),i=1..nops(L))+subs(m=n,f):
fi:


end:


###################
#####Problem 2#####
###################
#dIBVP(L,Ini,Fini,f,n,N) inputs lists L, Ini, and Fini,
#a function f in n, and a number N and outputs the N-th
#term of the solution of the linear recurrence equation
#with constant coefficients
#a(n)=L[1]*a(n-1)+L[2]*a(n-2)+ ... + L[-1]*a(n-nops(L))+f(n)
#and initial conditions
#a(0)=Ini[1], ..., a(nops(Ini)-1)=Ini[-1]
#a(N)=Fini[1], a(N-1)=Fini[2], ...

dIBVP:=proc(L,Ini,Fini,f,n,N) local eq,var,a,i,j:

if nops(L)<>nops(Ini)+nops(Fini) then
 RETURN(FAIL):
fi:

var:={seq(a[i],i=0..N)}:

eq:={seq(a[i]=Ini[i+1],i=0..nops(Ini)-1),
     seq(a[N-i]=Fini[i+1],i=0..nops(Fini)-1)}:

eq:=
eq union
{seq(a[i]-add(L[j]*a[i-j],j=1..nops(L))-subs(n=i,f),
i=nops(L)..N)}:

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
else
subs(var,[seq(a[i],i=0..N)]):
fi:

end:



###################
#####Problem 3#####
###################
#dBVP process needed for this part#

#dBVP(L,Ini,Fini,N): dBVP: Discrete boundary value problem
#The list of length N+1
#giving the values of a(n) for n=0 to n=N
#to the solution of the linear recurrence equation
#with constant coefficients
#a(n)=L[1]*a(n-1)+L[2]*a(n-2)+ ... + L[-1]*a(n-nops(L))
#a(0)=Ini[1], ..., a(nops(Ini)-1)=Ini[-1].
#a(N)=Fini[1], a(N-1)=Fini[2], ..
dBVP:=proc(L,Ini,Fini,N) local eq,var,a,i:

if nops(L)<>nops(Ini)+nops(Fini) then
 RETURN(FAIL):
fi:

var:={seq(a[i],i=0..N)}:

eq:={seq(a[i]=Ini[i+1],i=0..nops(Ini)-1),
     seq(a[N-i]=Fini[i+1],i=0..nops(Fini)-1)}:

eq:=
eq union
{seq(a[i]-add(L[j]*a[i-j],j=1..nops(L)), 
i=nops(L)..N)}:


var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
else
subs(var,[seq(a[i],i=0..N)]):
fi:

end:

#end dBVP process#



#GR(p,N) inputs a number p between 0 and 1 and a
#positive integer N and outputs a list of length
#N+1 so that G[i+1] gives you the probability that
#you exit the casino a winner if you enter with
#i dollars, repeatedly play a game where the
#probability of winning 1 dollar is p and the probability
#of losing 1 dollar is 1-p, and leave when you either have
#N dollars or nothing.

GR:=proc(p,N) local L,Ini,Fini:

#Let G be the list of winning probabilities, with
#G[i+1] being the probability of leaving with N dollars
#given you enter with i dollars.
#Clearly, G[1]=0 and G[N+1]=1.
#If we start with i dollars with 0<i<N, the probability 
#we win the first game is p, and so the probability
#we win the first hand and leave a winner is p*G[i+2].
#Meanwhile, the probability we lose the first hand is 1-p,
#and so the probability we lose the first hand and leave
#a winner is (1-p)*G[i]. So,
#G[i+1]=p*G[i+2]+(1-p)*G[i]
#(1/p)*G[i+1]=G[i+2]+((1-p)/p)*G[i]
#G[i+2]=(1/p)*G[i+1]+((p-1)/p)*G[i]

L:=[1/p,(p-1)/p]:
Ini:=[0]:
Fini:=[1]:

dBVP(L,Ini,Fini,N):

end:




###################
#####Problem 4#####
###################
#LGR(p,N) inputs a number p between 0 and 1 and a
#positive integer N and outputs a list of length
#N+1 so that G[i+1] gives you the expected duration
#of the game if you enter a casino with i dollars
#and repeatedly play a game where the probability
#of winning 1 dollar is p and the probability of
#losing 1 dollar is 1-p, and leave when you either have
#N dollars or nothing.

LGR:=proc(p,N) local Ini,Fini,f,L:

#Let G be the list of expected durations, with
#G[i+1] being the expected duration of a game
#when you enter with i dollars.
#Clearly, G[1]=0 and G[N+1]=0 since you are 
#already at an endpoint.
#If we start with i dollars with 0<i<N, the 
#probability we win the first game is p, and
#the expected duration if we win the first hand 
#is one more than the expected duration of 
#entering the game with i+1 dollars (since we 
#played a hand already).
#Meanwhile, the probability of losing the first 
#hand is 1-p, and the expected duration if we lose
#the first hand is one more than the expected 
#duration of entering the game with i-1 dollars.
#So,
#G[i+1]=p*(G[i+2]+1)+(1-p)*(G[i]+1)
#G[i+1]=p*G[i+2]+p+(1-p)*G[i]+1-p
#G[i+1]=p*G[i+2]+(1-p)*G[i]+1
#(1/p)*G[i+1]=G[i+2]+((1-p)/p)*G[i]+(1/p)
#G[i+2]=(1/p)*G[i+1]-((p-1)/p)*G[i]-(1/p)

Ini:=[0]:
Fini:=[0]:
f:=-1/p:
L:=[1/p,(p-1)/p]:

dIBVP(L,Ini,Fini,f,n,N):

end:



###################
#####Problem 5#####
###################
#A few iterations of the results of GR(1/2,N) and LGR(1/2,N)
#show a clear pattern that we can use to conjecture a closed form
#for the result.
#
#GR(1/2,N) looks to be the sequence:
#[seq(i/N,i=0..N)], or p=i/N when entering with i dollars
#
#
#LGR(1/2,N) appears to be the sequence:
#
#with(combinat)
#seq((N-1)*i - 2*nops(choose(i,2)),i=0..N)
#
#where nops(choose(i,2)) is the iC2 or (i)
#                                      (2)
#i.e. i!/(2!*(i-2)!)