#C21.txt, April 10, 2014



Help:=proc() :print(`DisOp1(Oper,D1,Sten,Sh,h,ord)`): 
print(`DisOp(Oper,D1,Sten,Sh,h)`):
print(`BVP3(Oper,D1,f,x,h1,Ini,Fini,x1)`):
print(`BVP5(Oper,D1,f,x,h1,Ini,Fini,x1)`):
print(`PDisOp1(Oper,D1,D2,Sten,Shx,Shy,h,ord)`):
print(`PDisOp(Oper,D1,D2,Sten,Shx,Shy,h, MaxOrd)`):
end:

#old stuff

#D1(f)=diff(f,x), Sh(f): f(x)->f(x+h)
#DisOp1(Oper,D1,Sten,Sh,h,ord):
#inputs a differential operator with constant coeff.
#a set of integers, Sten, a symbol for the shift operator
#Sh u(x)->u(x+h), a symbol for h, the mesh size
#and a pos. intetger ord
#Outputs: a recurrence operator (a Laurent poly in Sh)
#whose powers are in the Sten
#For example
#DisOp1(D1,D1,{0,1},Sh,h,1); should give
#(Sh-1)/h, 
#DisOp1(D1^2,D1,{-1,0,1},Sh,h,2); should give
#(Sh-2+1/Sh)/h^2
DisOp1:=proc(Oper,D1,Sten,Sh,h,ord) local eq,var,a,i, Oper1, Oper1a:
option remember:
Oper1:=add(a[i]*Sh^i, i in Sten):
var:={seq(a[i], i in Sten)}:
Oper1a:=subs(Sh=exp(h*D1),Oper1)-Oper:

Oper1a:=taylor(Oper1a,D1=0,ord+1):

eq:={seq(coeff(Oper1a,D1,i),i=0..ord)}:

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
else
 RETURN(subs(var,Oper1)):
fi:


end:


#inputs a differential operator with constant coeff.
#a set of integers, Sten, a symbol for the shift operator
#Sh u(x)->u(x+h), a symbol for h, the mesh size

#Outputs: a recurrence operator (a Laurent poly in Sh)
#whose powers are in the Sten, with the highest possible
#order, along with that order
#For example
#DisOp(D1,D1,{0,1},Sh,h); should give
#(Sh-1)/h, 
#DisOp1(D1^2,D1,{-1,0,1},Sh,h,2); should give
#(Sh-2+1/Sh)/h^2
DisOp:=proc(Oper,D1,Sten,Sh,h) local ord,Oper1:

for ord from 0 while DisOp1(Oper,D1,Sten,Sh,h,ord)<>FAIL do od:
DisOp1(Oper,D1,Sten,Sh,h,ord-1),ord-1:
end:

#end old stuff
#[0,h1,2*h1,3*h1, i*h1=x1,,,, N*h1]
#BVP3(Oper,D1,f,x,h1,x1): inputs
#(i) Oper: a second-order const. coeff. diff. oper
#phrased in terms of the diff. oper. D1
#(ii) an expression in the cont. variable x
#(iii) Sh, a symbol for the shift operator u(x)->u(x+h)
#(h is a symbol)
#(iv): a small NUMERIC mesh-size h1,
#(v): 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 3-point discretation of Oper
#with mesh size h1
BVP3:=proc(Oper,D1,f,x,h1,Ini,Fini,x1)
local T,i,eq,var,N,Oper1,i1,h,Sh:

N:=1/h1:

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,{-1,0,1},Sh,h)[1]:

for i from 1 to N-1 do

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

od:

var:=solve(eq,var):

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

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



end:


#BVP5(Oper,D1,f,x,h1,x1): inputs
#(i) Oper: a second-order const. coeff. diff. oper
#phrased in terms of the diff. oper. D1
#(ii) an expression in the cont. variable x
#(iii) Sh, a symbol for the shift operator u(x)->u(x+h)
#(h is a symbol)
#(iv): a small NUMERIC mesh-size h1,
#(v): 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 5-point discretation of Oper
#with mesh size h1
BVP5:=proc(Oper,D1,f,x,h1,Ini,Fini,x1)
local T,i,eq,var,N,Oper1,i1,h,Sh:

N:=1/h1:

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,{-2,-1,0,1,2},Sh,h)[1]:

for i from 2 to N-2 do

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

od:

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

od:

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

od:




var:=solve(eq,var):

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

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



end:


###PDO



#D1(f)=diff(f,x), Shx(f): f(x,y)->f(x+h,y)
#D2(f)=diff(f,y), Shy(f): f(x,y)->f(x,y+h)
#PDisOp1(Oper,D1,D2,Sten,Shx,Shy,h,ord):
#inputs a partial differential operator with constant coeff.
#a set of 2D lattice points, Sten, symbols for the shift operator
#Shx u(x,y)->u(x+h,y), and Shy 
#a symbol for h, the mesh size
#and a pos. intetger ord
#Outputs: a recurrence operator (a Laurent poly in Sh)
#whose powers are in the Sten
#For example
#PDisOp1(D1^2+D2^2,D1,D2,{[0,0],[-1,0],[1,0],[0,-1],[0,1]},
# Shx,Shy, h,3); should give
#(1/Shx+1/Shy+Shx+Shy-4)/h^2

PDisOp1:=proc(Oper,D1,D2,Sten,Shx,Shy,h,ord) 
local eq,var,a,i, j,Oper1, Oper1a,var1,i1:
option remember:

Oper1:=add(a[i]*Shx^i[1]*Shy^i[2], i in Sten):

var:={seq(a[i], i in Sten)}:

Oper1a:=subs({Shx=exp(h*D1), Shy=exp(h*D2)},  Oper1-Oper):

Oper1a:=taylor(Oper1a,D1=0,ord+1):

Oper1a:=add(coeff(Oper1a,D1,i1)*D1^i1,i1=0..ord):

Oper1a:=taylor(Oper1a,D2=0,ord+1):

Oper1a:=add(coeff(Oper1a,D2,i1)*D2^i1,i1=0..ord):



eq:={seq(seq(coeff(coeff(Oper1a,D1,i),D2,j),j=0..ord-i),
i=0..ord)}:

var:=solve(eq,var):


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

var1:={}:

for i from 1 to nops(var) do

if op(1,op(i,var))=op(2,op(i,var)) then
 var1:=var1 union {op(2,op(i,var))}:
fi:

od:


Oper1:=subs(var,Oper1):

Oper1:=subs({seq(var1[i]=0,i=1..nops(var1))},Oper1):

if Oper1=0 then
 RETURN(FAIL):
else
RETURN(Oper1):
fi:
end:


#inputs a differential operator with constant coeff.
#a set of integers, Sten, a symbol for the shift operator
#Sh u(x)->u(x+h), a symbol for h, the mesh size

#Outputs: a recurrence operator (a Laurent poly in Sh)
#whose powers are in the Sten, with the highest possible
#order, along with that order
#For example
#DisOp(D1,D1,{0,1},Sh,h); should give
#(Sh-1)/h, 
#DisOp1(D1^2,D1,{-1,0,1},Sh,h,2); should give
#(Sh-2+1/Sh)/h^2
PDisOp:=proc(Oper,D1,D2,Sten,Shx,Shy,h,MaxOrd) local ord,try1,ord1:

for ord from 1 to MaxOrd do
try1:= PDisOp1(Oper,D1,D2,Sten,Shx,Shy,h,ord):

if try1<>FAIL then

for ord1 from ord to MaxOrd while PDisOp1(Oper,D1,D2,Sten,Shx,Shy,h,ord1)<>FAIL do od:

RETURN( PDisOp1(Oper,D1,D2,Sten,Shx,Shy,h,ord1-1)):

fi:
od:

FAIL:

end: