#hw18.txt; Due April 6, 2014; Frank Wagner
Help:=proc(): print(` SIrecS(L,Ini,f,m,N), dBVPIc(L,Ini,Fini,f,m,N) `):
print(` P1dClever(f,x,h,mu0,mu1), GlEclever(f,x,h,mu0,mu1) `):
print(` EmpiricalStability(f,x,mu0,mu1,ListH) `): end:


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

##SIrec program from MY hw17.txt##

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


##end MY SIrec program##


#SIrecS(L,Ini,f,m,N) inputs a linear recurrence equation with
#constant coefficents, coded by L, a list of initial values Ini,
#a function f in m, and outputs the LIST of length N+1 such that
#L[i+1]=SIrec(L,Ini,f,m,i)
SIrecS:=proc(L,Ini,f,m,N) local n:

[ seq(SIrec(L,Ini,f,m,n),n=0..N) ]:
end:




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

#dBVPIc(L,Ini,Fini,f,m,N): a (hopefully) clever way to do
#dBVPI(...)
#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))+f(n)
#a(0)=Ini[1], ..., a(nops(Ini)-1)=Ini[-1]
#a(N)=Fini[1], a(N-1)=Fini[2], ...
dBVPIc:=proc(L,Ini,Fini,f,m,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:=SIrecS(L,Ini1,f,m,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:



###################
#####Problem 3#####
###################

#P1dClever(f,x,h,mu0,mu1) approximates the 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
#and h the mesh size.
#h^2*f(x)~u(x+h)-2*u(x)+u(x-h) -> u(x+h)=2*u(x)-u(x-2*h)+h^2*f(x)
#a(n)=2*a(n-1)-a(n-2)+h^2*f(n-1)
P1dClever:=proc(f,x,h,mu0,mu1) local Ini,Fini,Li,fy:

Li:=[2,-1]:
Ini:=[mu0]:
Fini:=[mu1]:

fy:=h^2*f:
fy:=subs(x=(x-1)*h,fy):

dBVPIc(Li,Ini,Fini,fy,x,1/h):

end:



###################
#####Problem 4#####
###################

##P1ch program from c18.txt##

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

##end P1ch program##


#GlEclever(f,x,h,mu0,mu1) inputs the same things as P1ch and P1dClever
#but gives the global error of the 1d Poisson eq.
GlEclever:=proc(f,x,h,mu0,mu1) local R1,A1:

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

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

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


end:



###################
#####Problem 5#####
###################

#EmpiricalStability(f,x,mu0,mu1,ListH) finds the sequence of ratios of
#G1EClever(f,x,h,mu0,mu1)/h^2
#for all h in ListH
EmpiricalStability:=proc(f,x,mu0,mu1,ListH):

[seq(GlEclever(f,x,ListH[i],mu0,mu1)/ListH[i]^2,i=1..nops(ListH))]:

end:


#EmpiricalStability(x^4,x,0,1,[1/10,1/20,1/30,1/40]); returns
##[977/25000, 188461/4800000, 29887/759375, 8061/204800] or
##[0.03908, 0.0392627083, 0.03935736626, 0.03936035156]
#EmpiricalStability(exp(x),x,0,1,[1/10,1/20,1/30,1/40]); returns
##[0.017526126, 0.017647883, 0.017651, 0.017648592]
#EmpiricalStability(1/(3-x),x,0,1,[1/10,1/20,1/30,1/40]);
##made Maple crash.
##Even GlEClever(1/(3-x),x,1/20,0,1) took far too long to execute.