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

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

Help:=proc():
    print(` Heq3(U,x,h,k,x1,t1) , Compare(U,t,x,h,k,x1,t1) , Dirichlet(f,x,y,h,k,x1,y1) `):
end:

##PROBLEM 1##

#Heq3(U,x,h,k,x1,t1): Crank-Nicolson
#inputs an expression U of the variable x
#two small numbers h and k (the mesh sizes in x and t resp.)
#and a final time t1
#outputs a list of lists such that
#L[j][n] is an APPROXIMATION to the value of u(x1, t1)
#of the temperature at x=j*h and t=n*k
#u(x,t) is the solution of the 1D-Heat equation
#u_{xx}=u_t, u(x,0)=U(x), u(0,t)=0, u(1,t)=0
Heq3:=proc(U,x,h,k,x1,t1) local T,j,n,N,J,r,eq,var:
    N:=t1/k:
    J:=1/h:
    r:=k/h^2:
    
    for j from 0 to J do
        T[j,0]:=subs(x=j*h,U):
    od:

    for n from 0 to N-1 do
        var:={seq(T[j,n+1], j=0..J)}:
        eq:={T[0,n+1]=0, T[J,n+1]=0}:
        eq:= eq union {seq((2+2*r)*T[j,n+1]-r*T[j-1,n+1]-r*T[j+1,n+1]=(2-2*r)*T[j,n]+r*T[j-1,n]+r*T[j+1,n],j=1..J-1)}:
        var:=solve(eq,var):
        for j from 0 to J do
            T[j,n+1]:=subs(var,T[j,n+1]):
        od:
    od:
    #T[i,j] is an approximation to u(i*h,k*j)
    return T[x1/h,t1/k]:
end:

##PROBLEM 2##

#Compare(U,t,x,h,k,x1,t1): inputs an expression U in the variables
#t and x, that you know beforehand satisfies the Heat Equation and
#the above boundary and initial conditions (it should return FAIL
#if it does not!), and ouputs, a list with FOUR values
#[The Exact Value, OutputOfHeEq1, OutputOfHeEq2,OutputOfHeEq3] 
#and an integer in the set {1,2,3} indicating who got closest
#to the exact value.
Compare:=proc(U,t,x,h,k,x1,t1) local L, aarghs, closest, dist, i, d:
    L:=[0,0,0,0]:
    L[1]:=evalf(subs({x=x1, t=t1}, U)):
    aarghs:=subs(t=0, U),x,h,k,x1,t1:
    L[2]:=evalf(Heq1(aarghs)):
    L[3]:=evalf(Heq2(aarghs)):
    L[4]:=evalf(Heq3(aarghs)):
    closest:=1:
    dist:=abs(L[2]-L[1]):
    for i from 3 to 4 do
        d:=abs(L[i]-L[1]):
        if d < dist then
            dist:=d:
            closest:=i-1:
        fi:
    od:
    return L, closest:
end:

##PROBLEM 3##

#Dirichlet(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
Dirichlet:=proc(f,x,y,h,k,x1,y1) local T, i, j, var, eq:
    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
    eq:=eq union {seq(seq(T[i][j]=(k^2*(T[i+1][j]+T[i-1][j])+h^2*(T[i][j+1]+T[i][j-1])-h^2*k^2*subs({x=i*h,y=j*k}, f))/(2*h^2+2*k^2), j=1..1/k-1), i=1..1/h-1)}:
    var:=solve(eq, var):
    T[x1/h][y1/k]:=subs(var, T[x1/h][y1/k]):
    return T[x1/h][y1/k]:
end:

#evalf(Dirichlet(x^2+y^2,x,y,1/10,1/10,1/2,1/2)); returned -4.184843553, which I know is wrong
#The other two took too long