#March 31, 2014, C18.txt difference methods
Help:=proc(): print(` Srec(L, Ini, n), dBVP(L,Ini,Fini,N) `): 
print(`SrecS(L,Ini,N) , dBVPc(L,Ini,Fini,N), P1ch(f,x,h,mu0,mu1)`):
print(`P1d(f,x,h,mu0,mu1), GlE(f,x,h,mu0,mu1) `):
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
 expand(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 old stuff


#SrecS(L,Ini,N): inputs a linear recurrence equation
#with constant coefficients, coded by L
#and a list of initial values, Ini
#outputs the LIST of length N+1 such that
#L[i+1]=Srec(L,Ini,n)
SrecS:=proc(L,Ini,N) local n:
[ seq(Srec(L,Ini,n),n=0..N)]:
end:

#dBVPc(L,Ini,Fini,N): a (hopefully) clever way
#to do dBVP(...)
#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], ..
dBVPc:=proc(L,Ini,Fini,N) local eq,var,a,i,Ini1, Li:

Ini1:=[op(Ini),seq(a[i],i=1..nops(Fini))]:

var:={seq(a[i],i=1..nops(Fini))}:

Li:=SrecS(L,Ini1,N):

#Li[N+1]=Fini[1], Li[N]=Fini[2], Li[N-1]=Fini[3]:
eq:={seq(Li[N+2-i]=Fini[i], i=1..nops(Fini))}:

var:=solve(eq,var):

subs(var,Li):


end:






#P1ch(f,x,h,mu0,mu1): the exact solution of the 1-D Poisson eq.
#u''(x)=f(x), u(0)=mu0, u(1)=mu1  
#where f is an expression in the variable x
P1ch:=proc(f,x,h,mu0,mu1) local u,c1,c2,eq,var:
u:=int(int(f,x)+c1,x)+c2:
var:={c1,c2}:
eq:={subs(x=0,u)=mu0, subs(x=1,u)=mu1}:
var:=solve(eq,var):

u:=subs(var,u):


[seq(subs(x=i*h,u),i=0..1/h)]:

end:

#P1d(f,x,h,mu0,mu1): approximate solution of the 1-D Poisson eq.
#u''(x)=f(x), u(0)=mu0, u(1)=mu1  
#where f is an expression in the variable x
#h the mesh size, and
#u''(x) ~ u(x+h)-2*u(x)+u(x-h)
P1d:=proc(f,x,h,mu0,mu1) local u,Li,var,eq:

Li:=[seq(u[i],i=0..1/h)]:

var:={seq(u[i],i=0..1/h)}:

eq:={ seq(u[i+1]-2*u[i]+u[i-1]=h^2*subs(x=i*h,f), i=1..1/h-1)}:

eq:=eq union {u[0]=mu0, u[1/h]=mu1}:

var:=solve(eq,var):

subs(var,Li):

end:



#Global error in using  a mesh size h
#to solve the 1D Poisson eq.
#u''(x)=f(x), u(0)=mu0, u(1)=mu1
GlE:=proc(f,x,h,mu0,mu1) local R1, A1:

R1:=P1ch(f,x,h,mu0,mu1):

A1:=P1d(f,x,h,mu0,mu1):

max(seq(abs(R1[i]-A1[i]),i=1..nops(R1))):


end: