#Kafung Mok
#HW 17
#March 30, 2014


############################################## Class Notes March 27 ##################################################
#C17.txt: March 27, 2014
#Discrete Calculus
Help:=proc(): print(` GPi(N) , Dsin(x,h) , Srec(L,Ini,n)`): 
print(` dBVP(L,Ini,Fini,N) `):
end:

#GPi(N): the geometrical discrete Pi
#the number of lattice points (x,y) such
#that x^2+y^2<N divided by N^2
GPi:=proc(N) local co,x,y:

co:=0:

for x from 0 to N do
 for y from 1 while x^2+y^2<N^2 do 
  co:=co+1:
 od:
od:

(4*co+1)/N^2:

end:

#Pi is the smallest positive "real" number such that
#sin x=0,  y''(x)+y(x)=0, y(0)=0, y'(0)=1
#y(0)=0, (y(h))/h=1
#Dsin(x,h): the discrete analog of the sine function
#with mesh-size h
#Dsin(x,h): (y(x+h)-2*y(x)+y(x-h))/h^2 +y(x)=0
#(y(x+2*h)-2*y(x+h)+y(x))/h^2 +y(x+h)=0
##y(x+2*h)=(2-h^2)*y(x+h)-y(x)
###y(x)=(2-h^2)*y(x-h)-y(x-2*h))
Dsin:=proc(x,h) option remember:
if x=0 then
 0:
elif x=h then 
h:
else
 expand((2-h^2)*Dsin(x-h,h)-Dsin(x-2*h,h)):
fi:

 
end:

#Srec(L,Ini,n): 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))
#a(0)=Ini[1], ..., a(nops(Ini)-1)=Ini[-1]
Srec:=proc(L,Ini,n) local i:
option remember:

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]*Srec(L,Ini,n-i),i=1..nops(L)):
fi:


 
end:

#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 of Class Notes March 27 ###############################################


#############
# Problem 1 #
#############


#SIrec(L,Ini,f,n,N): has the same input as Srec but in addition 
#an expression, f, in the discerte variable n, that solves 
#Inhomogeneous linear recurrences 
#a(n)=L[1]*a(n-1)+L[2]*a(n-2)+ ... + L[-1]*a(n-nops(L))+f(n) 

SIrec:=proc(L,Ini,f,n,N) local i:
option remember:

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,n,N-i),i=1..nops(L))+subs(n=N,f):
fi:

 
end:




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


#dBVP(L,Ini,Fini,f,n,N): has the same input as dBVP(L,Ini,Fini,N)
#but in addition an expression, f, in discrete variable n.
dBVP:=proc(L,Ini,Fini,f,n,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))-subs(n=N,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 #
#############


#GR(p,N):inputs a real number p, between 0 and 1 (or symbolic), and a positive integer N, and outputting 
#a list of length N+1, let's call it G, such G[i+1] is your chance of exiting a winner in a casino with 
#the following rule: if you enter a casino with i dollars where at each step you win a dollar with 
#probability p and lose a dollar with probability 1-p, and you stop playing when you are either 
#broke (0 dollars), or got N dollars.
#Note:G[1]=0 and G[N+1]=1

GR:=proc(p,N) 
G:=dBVP([1/p,-(1-p)/p],[0],[1],N) :
end: