#hw21.txt.  April 12, 2014.  Katie McKeon

read(`hw20.txt`):
read(`C21.txt`):

Help:=proc()
print(`BVPk(Oper,k,D1,f,x,h1,Ini,Fini,x1)`):
print(`Compare(Oper,D1,f,x,h1,Ini,Fini,x1)`):
print(`SuperDirichlet(f,x,y,h,k,x1,y1) `):
end:



#BVPk(Oper,k,D1,f,x,h1,Ini,Fini,x1): inputs
#(i) Oper: a second-order const. coeff. diff. oper
#phrased in terms of the diff. oper. D1
#(ii) an integer k representing the stencil {-k..0..k}
#(iii) an expression f in the cont. variable x
#(iv) Sh, a symbol for the shift operator u(x)->u(x+h)
#(h is a symbol)
#(v): a small NUMERIC mesh-size h1,
#(vi): x1 a number between 0 and 1 ( a multiple of h1)
#OUTPUTS
#an appx. to u(x1) where u(x) is a solution BVP
#Oper(D1)u(x)=f(x), f(0)=Ini, f(1)=Fini
#using the Central 2k+1-point discretation of Oper
#with mesh size h1
BVPk:=proc(Oper,k,D1,f,x,h1,Ini,Fini,x1)
local T,i,eq,var,N,Oper1,i1,h,Sh,Sten:

N:=1/h1:
Sten:={seq(i,i=-k..k)}:

if degree(Oper,D1)<>2 then
 RETURN(FAIL):
fi:

var:={seq(T[i],i=0..N)}:

eq:={T[0]=Ini,T[N]=Fini}:

Oper1:=DisOp(Oper,D1,Sten,Sh,h)[1]:

for i from k to N-k do

eq:=eq union 
{subs(h=h1,add(coeff(Oper1,Sh,i1)*T[i+i1],i1 in Sten))=
               subs(x=i*h1,f)}:

od:

for i from 1 to k-1 do
Sten:={seq(i1,i1=-i..2*k-i)}:
Oper1:=DisOp(Oper,D1,Sten,Sh,h)[1]:
eq:=eq union 
{subs(h=h1,add(coeff(Oper1,Sh,i1)*T[i+i1],i1 in Sten))=
               subs(x=i*h1,f)}:

od:

for i from N-k+1 to N-1 do
Sten:={seq(i1,i1=-2*k+N-i..N-i)}:
Oper1:=DisOp(Oper,D1,Sten,Sh,h)[1]:
eq:=eq union 
{subs(h=h1,add(coeff(Oper1,Sh,i1)*T[i+i1],i1 in Sten))=
               subs(x=i*h1,f)}:


od:

var:=solve(eq,var):

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

subs(var,T[x1/h1]):

end:




#Compare(Oper,D1,f,x,h1,Ini,Fini,x1) outputs the exact value u(x1)
#where u(x) is the solution to the BVP Oper u(x)=f(x) with initial
#values at 0 Ini and end values at 1 Fini along with the quadruple
#of estimations for u(x1) using BVPk with k=1,2,3,4
Compare:=proc(Oper,D1,f,x,h1,Ini,Fini,x1) local k:

[BVPexact(Oper,D1,f,x,0,1,[Ini],[Fini],x1), seq(BVPk(Oper,k,D1,f,x,h1,Ini,Fini,x1),k=1..4)]

end:


#DATA
#evalf(Compare(D1^2+4*D1+11,D1,x^5,x,1/20,1,2,1/2))
#[11.74972904, 11.94283371, 11.75012878, 11.74971554, 11.74972885]

#evalf(Compare(D1^2-6*D1+7,D1,exp(x),x,1/20,2,3,1/2)) 
#[3.635999802, 3.642310609, 3.635884948, 3.636000947, 3.635999797]





#PDisOp(D1^2+D2^2,D1,D2,{seq(seq([i1, j1],i1=-2..2),j1=-2..2)},Shx,Shy,h,50) returns 
Stencil:=1/(2*h^2*Shx^2)+1/(2*h^2*Shy^2)+1/(6*h^2*Shx^2*Shy^2)-
1/(2*h^2*Shx^2*Shy)-1/(2*h^2*Shx*Shy^2)+(2/3)*Shx/(h^2*Shy)+
(2/3)*Shy/(h^2*Shx)+(1/12)*Shx^2/(h^2*Shy^2)+(1/12)*Shy^2/(h^2*Shx^2)-
(1/6)*Shy^2/(h^2*Shx)-(1/6)*Shx^2/(h^2*Shy)-(1/3)*Shx/(h^2*Shy^2)-
(1/3)*Shy/(h^2*Shx^2)+4/(3*h^2*Shx*Shy)-4/h^2+Shy/h^2+Shx/h^2:
#We factor out 1/(12*h^2) and rewrite this stencil in SuperDirichlet


#SuperDirichlet(f,x,y,h,k,x1,y1) inputs a function f(x,y) and approximates
#u(x1,y1) (where u is the solution to the Dirichlet problem)
#using the larger stencil {seq(seq([i1, j1],i1=-2..2),j1=-2..2)} for
#interior points
SuperDirichlet:=proc(f,x,y,h,k,x1,y1)local T,j,n,N,J,r,eq,var:

J:=1/h:
N:=1/h:

for j from 0 to J do
T[j,0]:=0: T[j,N]:=0:
od:

for n from 0 to N do
T[0,n]:=0: T[J,n]:=0:
od:

var:={seq(seq(T[j,n], j = 1..J-1), n = 1..N-1)}:

eq:={seq(seq(subs({x = h*j, y = h*n},f)*12*h^2 =

6*T[j-2, n]+6*T[j, n-2]+2*T[j-2, n-2]-6*T[j-2, n-1]-6*T[j-1, n-2]+8*T[j+1, n-1]+8*T[j-1, n+1]+T[j+2, n-2]+T[j-2, n+2]-2*T[j-1, n+2]-2*T[j+2, n-1]-4*T[j+1, n-2]-4*T[j-2, n+1]+16*T[j-1, n-1]-48*T[j, n]+12*T[j, n+1]+12*T[j+1, n],

j = 2..J-2), n = 2..N-2)}:



#Using the five point stencil for the boundary points
j:=1:
eq:=eq union {seq(subs({x = h*j, y = h*n},f) =

(T[j+1,n]-2*T[j,n]+T[j-1,n])/h^2 + (T[j,n+1]-2*T[j,n]+T[j,n-1])/h^2,

n = 1..N-1)}:


j:=J-1:
eq:=eq union {seq(subs({x = h*j, y = h*n},f) =

(T[j+1,n]-2*T[j,n]+T[j-1,n])/h^2 + (T[j,n+1]-2*T[j,n]+T[j,n-1])/h^2,

n = 1..N-1)}:

n:=1:
eq:=eq union {seq(subs({x = h*j, y = h*n},f) =

(T[j+1,n]-2*T[j,n]+T[j-1,n])/h^2 + (T[j,n+1]-2*T[j,n]+T[j,n-1])/h^2,

j = 2..N-2)}:

n:=N-1:
eq:=eq union {seq(subs({x = h*j, y = h*n},f) =

(T[j+1,n]-2*T[j,n]+T[j-1,n])/h^2 + (T[j,n+1]-2*T[j,n]+T[j,n-1])/h^2,

j = 2..N-2)}:


var:=solve(eq, var):

subs(var, T[x1/h, y1/k]):

end:

#DATA
#evalf(subs({x = 1/2, y = 1/2}, sin(Pi*x)*sin(Pi*y)));
#	1
#evalf(SuperDirichlet(-2*sin(Pi*x)*Pi^2*sin(Pi*y), x, y, 1/10, 1/10, 1/2, 1/2));
#	0.9998528807
#evalf(Dirichlet(-2*sin(Pi*x)*Pi^2*sin(Pi*y), x, y, 1/10, 1/10, 1/2, 1/2));
#	1.008265418

#f := sin(Pi*x)*sin(Pi*y)*exp(x*y):
#evalf(subs({x = 1/2, y = 1/2}, f));
#	1.284025417
#evalf(SuperDirichlet(diff(f, x, x)+diff(f, y, y), x, y, 1/10, 1/10, 1/2, 1/2));
#	1.283893417
#evalf(Dirichlet(diff(f, x, x)+diff(f, y, y), x, y, 1/10, 1/10, 1/2, 1/2));
#	1.293297089