#Ross Berkowitz  hw26.txt
#Posting Permission Granted

Help:=proc():
print("ExpD2(M,N)"):
end:

ExpD2:=proc(N,M) local a,eq,var,i,j:

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

eq:={seq(a[i][0],i=0..N),seq(a[0][j],j=0..M),
	 seq(a[i][M],i=0..N),seq(a[N][j],j=0..M),
		seq(seq(a[i][j]=1+1/4*(	a[i-1][j]+
								a[i][j-1]+
								a[i+1][j]+
								a[i][j+1]
							),j=1..M-1)
						,i=1..N-1)
	}:
	
var:=solve(eq,var):
if var=NULL then
	RETURN(FAIL):
fi:
subs(var,[seq( [seq(a[i][j],j=0..M)], i=0..N)]):
end:
















###################   OLD STUFF  ###################

#C26.txt, Gambler's ruin, April 28, 2014 
Help26:=proc(): print(`Die1(p), GR(i,N,p), ExpD(N,p)`): 
print(`ExpDf(N,p) , PGF(N,p,t), PGFf(N,p,t), GuessP(L,x) `):
end:

#Die1: inputs p rational between 0 and 1
#and outputs 1 with prob. p and 0 with prob. 1-p
Die1:=proc(p) local a,b,ra:
a:=numer(p):
b:=denom(p):
ra:=rand(1..b)():
if ra<=a then
 1:
else
 -1:
fi:

end:


#GR(i,N,p): simulates a gambler's ruin session
#staring with i dollars, ending with 0 or N dollars
#and with prob. of winning a dollar =p
GR:=proc(i,N,p) local x,j:

x:=i:

for j from 0 while (x>0 and x<N) do
 x:=x+Die1(p):
od:

x,j:

end:

#ExpD(N,p): Let f(i) be the expected duration of the
#game starting at i, outputs the vector
#[f(1),f(2),..., f(N-1)] of expected durations with
#prob. p of winning one dollar
#f(i)=p*f(i+1)+(1-p)*f(i-1)+1, i=1..N-1, f(0)=0, f(N)=0
ExpD:=proc(N,p) local eq,var,f,i,var1:
var:=[seq(f[i],i=1..N-1)]:
eq:={seq(f[i]=p*f[i+1]+(1-p)*f[i-1]+1,i=1..N-1)}:
eq:=subs({f[0]=0,f[N]=0},eq):
var1:=solve(eq,convert(var,set)):
subs(var1,var):
end:


#ExpD(N,p): fast version;
#Let f(i) be the expected duration of the
#game starting at i, outputs the vector
#[f(1),f(2),..., f(N-1)] of expected durations with
#prob. p of winning one dollar
#f(i)=(f(i-1)-(1-p)*f(i-2)-1)/p, i=1..N-1, f(0)=0, f(N)=0
ExpDf:=proc(N,p) local f,i,a,a1:
f[1]:=a:
f[0]:=0:
for i from 2 to N do
f[i]:=normal((f[i-1]-(1-p)*f[i-2]-1)/p):
od:

a1:=solve(f[N]=0,a):

[seq(subs(a=a1,f[i]),i=1..N-1)]:

end:


#PGF(N,p,t): Let f(i,t) be the prob. generating function of the
#a rational function in t whose Maclaurin expansion's
#coeff. of t^j
#gives the exact prob. of finishing in j steps
#game starting at i, outputs the vector
#[f(1),f(2),..., f(N-1)] of expected durations with
#prob. p of winning one dollar
#f(i,t)=t*(p*f(i+1,t)+(1-p)*f(i-1,t)), i=1..N-1, f(0)=0, f(N)=0
PGF:=proc(N,p,t) local eq,var,f,i,var1:
var:=[seq(f[i],i=1..N-1)]:
eq:={seq(f[i]=t*(p*f[i+1]+(1-p)*f[i-1]),i=1..N-1)}:
eq:=subs({f[0]=1,f[N]=1},eq):
var1:=solve(eq,convert(var,set)):
subs(var1,var):
end:


#PGFf(N,p,t): fast version;
#f(i,t)=t*(p*f(i+1,t)+(1-p)*f(i-1,t))
#f(i,t)= (f(i-1,t)/t-(1-p)*f(i-2,t))/p
PGFf:=proc(N,p,t) local f,i,a,a1:
f[1]:=a:
f[0]:=1:
for i from 2 to N do
f[i]:= normal((f[i-1]/t-(1-p)*f[i-2])/p):
od:

a1:=solve(f[N]=1,a):

[seq(subs(a=a1,f[i]),i=1..N-1)]:

end:

# begin Added after class
GuessP1:=proc(L,x,d) local i,P,a,var,eq:


if nops(L)-d<2 then
RETURN(FAIL):
fi:

P:=add(a[i]*x^i,i=0..d):

var:={seq(a[i],i=0..d)}:
eq:={seq(subs(x=i,P)=L[i],i=1..nops(L))}:

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
fi:

subs(var,P):

end:

#GuessP(L,x): guesses a polynomial expression in x for the list L  such that L[i]=P(i)
GuessP:=proc(L,x)  local d,gu:

for d from 0 to nops(L)-3 do
gu:=GuessP1(L,x,d):
if gu<>FAIL then
  RETURN(gu):
fi:
od:

FAIL:

end:
#end Added after class