#Nathan Fox
#Homework 21
#I give permission for this work to be posted online

#Read procedures from class
read(`C21.txt`):
read(`hw20.txt`):
read(`C19.txt`):

Help:=proc():
    print(` BVPk(Oper,k,D1,f,x,h1,Ini,Fini,x1) , Compare(Oper,D1,f,x,h1,Ini,Fini,x1,K:=4) , PDisOp1(Oper,D1,D2,Sten,Shx,Shy,h,k,ord) , PDisOp(Oper,D1,D2,Sten,Shx,Shy,h,k:=false) , SuperDirichlet(f,x,y,h,k,x1,y1) `):
end:

##PROBLEM 1##

#BVPk(Oper,k,D1,f,x,h1,Ini,Fini,x1): generalizes
#BVP3(Oper,D1,f,x,h1,Ini,Fini,x1) and
#BVP5(Oper,D1,f,x,h1,Ini,Fini,x1)
#by using, in the middle the stencil
#{-k,-(k-1), ..., k-1, k}
#and adjusting near the boundary such that 
#BVPk(Oper,1,D1,f,x,h1,Ini,Fini,x1) and
#BVPk(Oper,2,D1,f,x,h1,Ini,Fini,x1) 
#output the same thing as BVP3(Oper,D1,f,x,h1,Ini,Fini,x1) and
#BVP5(Oper,D1,f,x,h1,Ini,Fini,x1) respectively.
BVPk:=proc(Oper,k,D1,f,x,h1,Ini,Fini,x1)
 local T, i, i1, eq, var, N, Oper1, Sh, h:
    if degree(Oper, D1) <> 2 then
        return FAIL:
    fi:
    
    N:=1/h1:
    var:={seq(T[i], i=0..N)}:
    eq:={T[0]=Ini, T[N]=Fini}:
    
    Oper1:=DisOp(Oper,D1,{seq(i, i=-k..k)},Sh,h)[1]:
    
    for i from k to N-k do
        eq:=eq union {add(subs(h=h1, coeff(Oper1, Sh, i1))*T[i+i1], i1=-k..k)=subs(x=i*h1, f)}:
    od:
    
    for i from 1 to k-1 do
        Oper1:=DisOp(Oper, D1, {seq(i1, i1=-i..-1+2*k)},Sh,h)[1]:
        eq:=eq union {add(subs(h=h1, coeff(Oper1, Sh, i1))*T[i+i1], i1=-i..-1+2*k)=subs(x=i*h1, f)}:
    od:
    
    for i from N-k+1 to N-1 do
        Oper1:=DisOp(Oper, D1, {seq(i1, i1=N-i-2*k..N-i)},Sh,h)[1]:
        eq:=eq union {add(subs(h=h1, coeff(Oper1, Sh, i1))*T[i+i1], i1=N-i-2*k..N-i)=subs(x=i*h1, f)}:
    od:
    
    var:=solve(eq, var):
    
    if var = NULL then
        return FAIL:
    fi:
    
    return subs(var, T[x1/h1]):
end:

##PROBLEM 2##

#Compare(Oper,D1,f,x,h1,Ini,Fini,x1,K:=4): outputs the (K+1)-tuple
#[Exact Answer, Using k=1, Using k=2, ..., Using k=K] 
#with BVPk(Oper,k,D1,f,x,h1,Ini,Fini,x1) for k=1..K
#Note: K defaults to 4
Compare:=proc(Oper,D1,f,x,h1,Ini,Fini,x1,K:=4) local k:
    return [BVPexact(Oper,D1,f,x,0,1,[Ini],[Fini],x1), seq(BVPk(Oper,k,D1,f,x,h1,Ini,Fini,x1), k=1..K)]:
end:

#evalf(Compare(D1^2+4*D1+11,D1,x^5,x,1/20,1,2,1/2)); returned [11.74972904, 11.94283371, 11.75012878, 11.74971582, 11.74972885]
#evalf(Compare(D1^2-6*D1+7,D1,exp(x),x,1/20,2,3,1/2)); returned [3.635999802, 3.642310610, 3.635884948, 3.636001020, 3.635999793]

##PROBLEM 3##

#SuperDirichlet(f,x,y,h,k,x1,y1): inputs an expression f in the
#variables x and y, small positive numbers h and k, and numbers
#x1, y1 between 0 and 1 (that must be multiples of h and k
#respectively), and outputs a finite-difference approximation to
#u(x1,y1), using mesh size h in the x direction, and mesh-size k
#in the y-direction, and u(x,y) is the solution of the
#Poisson equation
#diff(u,x,x) + diff(u,y,y) = f
#in the unit square [0,1]x[0,1], that vanishes on the four sides
#of the unit square.
#NOTE: I expect that this procedure will be painfully slow
#NOTE: Either this doesn't work at all, or it's extremely unstable
#(using test function x*(1-x)*y*(1-y), which has Laplacian
#2*x^2-2*x+2*y^2-2*y)
#NOTE: My version of PDisOp takes k as an argument; I have included it at the bottom of this file for reference
SuperDirichlet:=proc(f,x,y,h,k,x1,y1)
local T, i, j, var, eq, Oper, BOperx, BOperX, BOpery, BOperY, BOperxy, BOperxY, BOperXy, BOperXY, Shx, Shy, operget, oper, n, expr, dx, dy:
    
    Oper:=PDisOp(D1^2+D2^2,D1,D2,{seq(seq([i,j],i=-2..2),j=-2..2)},Shx,Shy,h,k):
    BOperx:=PDisOp(D1^2+D2^2,D1,D2,{seq(seq([i,j],i=-1..3),j=-2..2)},Shx,Shy,h,k):
    BOperX:=PDisOp(D1^2+D2^2,D1,D2,{seq(seq([i,j],i=-3..1),j=-2..2)},Shx,Shy,h,k):
    BOpery:=PDisOp(D1^2+D2^2,D1,D2,{seq(seq([i,j],i=-2..2),j=-1..3)},Shx,Shy,h,k):
    BOperY:=PDisOp(D1^2+D2^2,D1,D2,{seq(seq([i,j],i=-2..2),j=-3..1)},Shx,Shy,h,k):
    BOperxy:=PDisOp(D1^2+D2^2,D1,D2,{seq(seq([i,j],i=-1..3),j=-1..3)},Shx,Shy,h,k):
    BOperXy:=PDisOp(D1^2+D2^2,D1,D2,{seq(seq([i,j],i=-3..1),j=-1..3)},Shx,Shy,h,k):
    BOperxy:=PDisOp(D1^2+D2^2,D1,D2,{seq(seq([i,j],i=-1..3),j=-1..3)},Shx,Shy,h,k):
    BOperXY:=PDisOp(D1^2+D2^2,D1,D2,{seq(seq([i,j],i=-3..1),j=-3..1)},Shx,Shy,h,k):
    
    operget:=proc(i, j)
        if i = 1 and j = 1 then
            return BOperxy:
        elif i = 1 and j = 1/k - 1 then
            return BOperxY:
        elif i = 1/h - 1 and j = 1 then
            return BOperXy:
        elif i = 1/h - 1 and j = 1/k - 1 then
            return BOperXY:
        elif i = 1 then
            return BOperx:
        elif i = 1/h - 1 then
            return BOperX:
        elif j = 1 then
            return BOpery:
        elif j = 1/k - 1 then
            return BOperY:
        else
            return Oper:
        fi:
    end:
    
    var:={seq(seq(T[i][j], j=0..1/k), i=0..1/h)}:
    #Boundary conditions
    eq:={seq(op([T[i][0]=0, T[i][1/k]=0]), i=0..1/h)}:
    eq:=eq union {seq(op([T[0][j]=0, T[1/h][j]=0]), j=0..1/k)}:
    #Interior conditions
    for i from 1 to 1/h - 1 do
        for j from 1 to 1/k - 1 do
            oper:=[op(operget(i, j)[1])]:
            expr:=0:
            for n from 1 to nops(oper) do
                dx:=degree(oper[n], Shx):
                dy:=degree(oper[n], Shy):
                expr:=expr + coeff(oper[n], Shx^dx*Shy*dy)*T[i+dx][j+dy]:
            od:
            eq:=eq union {subs({x=i*h,y=j*k}, f)=expr}:
        od:
    od:
    var:=solve(eq, var):
    T[x1/h][y1/k]:=subs(var, T[x1/h][y1/k]):
    return T[x1/h][y1/k]:
end:

##NOTE: This is insane, and I couldn't get it to work.  I'm not sure which procedure doesn't work (SuperDirichlet, PDisOp, or PDisOp1), but something is going horribly wrong

##MY VERSION OF PDisOp##

#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+k)
#PDisOp1(Oper,D1,D2,Sten,Shx,Shy,h,k,ord):
#inputs a partial differential operator with constant coefficients
#a set of 2D lattice points, Sten, symbols for the shift operator
#Shx u(x,y)->u(x+h.u) and Shy u(x,y)->u(x,y+h), symbols for h and k, the mesh sizes,
#and a positive integer ord
#outputs: a recurrence operator (a Laurent polynomial in Shx and Shy)
#whose powers are in the set Sten
#For example
#PDisOp1(D1^2+D2^2, D1, D1, {[0,0], [-1,0], [1,0], [0,-1],[0,1]}, Shx, Shy, h, h, 3); should give
#(1/Shx+1/Shy+Shx+Shy-4)/h^2
PDisOp1:=proc(Oper, D1, D2, Sten, Shx, Shy, h, k, ord) local eq, var, a, i, j, Oper1, Oper1a:
    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(k*D2)}, Oper1)-Oper:
    print(Oper1a):
    
    Oper1a:=taylor(Oper1a, D1=0, ord+1):
    Oper1a:=convert(Oper1a, polynom):
    Oper1a:=taylor(Oper1a, D2=0, ord+1):
    Oper1a:=convert(Oper1a, polynom):
    
    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:
    
    return subs(var, Oper1):
end:

#PDisOp(Oper,D1,D2,Sten,Shx,Shy,h,k:=false):
#inputs a partial differential operator with constant coefficients
#a set of lattice points, Sten, symbols for the shift operators
#Shx u(x,y)->u(x+h,y) and Shy u(x,y)->u(x,y+k), and symbols for h and k, the mesh sizes
#outputs: a recurrence operator (a Laurent polynomial in Shx and Shy)
#with the highest possible order
#whose powers are in the set Sten
#also returns the order
#Note if k receives false, uses h for k
PDisOp:=proc(Oper, D1, D2, Sten, Shx, Shy, h, k:=false) local ord, Oper1, Oper2, k2:
    if k = false then
        k2:=h:
    else
        k2:=k:
    fi:
    for ord from 1 while true do
        Oper1:=PDisOp1(Oper, D1, D2, Sten, Shx, Shy, h, k2, ord):
        if Oper1 = FAIL then
            return Oper2, ord-1:
        fi:
        Oper2:=Oper1:
    od:
end: