#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: