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

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

Help:=proc():
    print(` Pw(Lf,x,Lx) , Compare(p,q,v,x,x1,L:=[10,20,30,40]) , BSpline(d,Lx,x,i) , TestFS(c,Lx,x) , RR3(p,q,f,x,Lx) `):
end:

##PROBLEM 1##

#Pw(Lf,x,Lx): inputs a piecwise function in our format, where Lf
#is a list of expressions in the variable x, and Lx is a list of
#increasing numbers, and outputs it in Maple's notation.
#Anything out of bounds should be 0
#For example
#Pw([x,x^2],x,[0,1,2])
#should return
#piecewise(x<0, 0, x>=0 and x<1, x, 1<=x and x<2, x^2, x>=2, 0)
Pw:=proc(Lf, x, Lx) local i:
    if nops(Lf) + 1 <> nops(Lx) then
        return FAIL:
    fi:
    return piecewise(x<Lx[1], 0, seq(op([x>=Lx[i] and x<Lx[i+1],
        Lf[i]]), i=1..nops(Lf)), x>=Lx[nops(Lx)], 0):
end:

##PROBLEM 2##

#Compare(p,q,v,x,x1,L): inputs expressions p, q, v in x, and a
#number x1 between 0 and 1, then constructs u(x):=v(x)*x*(1-x) and
#f(x):= -(p(x)*u'(x))'+ q(x)*u(x)
#(so we know a priori the exact solution of the BVP
#-(p(x)*y'(x))'+ q(x)*y(x)=f(x) , 0<=x<=1, y(0)=0, y(1)=0
#namely u(x))
#and then uses
#RR(p,q,f,x,[seq(i/N,i=0..N)]):
#with N ranging over L, and evluates them at x1. The output is a
#list of length nops(L)+1 with the last entry being the
#exact value.
#NOTE: L defaults to [10,20,30,40]
Compare:=proc(p, q, v, x, x1, L:=[10,20,30,40]) local u, f, N, i:
    u:=v*x*(1-x):
    f:=-diff(p*diff(u, x), x) + q*u:
    return [seq(subs(x=x1, RR(p, q, f, x,
        [seq(i/N, i=0..N)])[floor(x1*N)+1]), N in L),
        subs(x=x1, u)]:
end:

#evalf(Compare(x^2,1+x^3,sin(x),x,.513)); returned [0.123687, 0.122767, 0.122639, 0.12265444, 0.1226153759]
#evalf(Compare(x^2,1+x^3,x^7+1,x,.513)); returned [0.2540307694, 0.2524677026, 0.2522576258, 0.2522113977, 0.2521669733]
#evalf(Compare(x^2+3,sin(x),exp(x),x,.542)); took too long

##PROBLEM 3##

#BSpline(d,Lx,i): inputs a positive integer d, and a list
#Lx[x0,...,x(d+1)]
#of numbers (increasing), outputs the B-spline of degree d
#on these points
#i.e. it is 0 for x<x0 and x>xd, and continuous to order d-1
#(d+1)*(d+1) unknown numbers,
#(d+2)*d equations
#and f(Lx[i+1])=1
#NOTE: x needed to be an argument, so I made it one
BSpline:=proc(d, Lx, x, i) local a, S, i1, j1, var, eq, diff2:
    diff2:=proc(u, x, n)
        if n > 0 then
            return diff(u, x$n):
        else
            return u:
        fi:
    end:
    if nops(Lx) <> d+2 then
        return FAIL:
    fi:
    #Degree d polynomials
    S:=[seq(add(a[i1,j1]*x^j1, j1=0..d), i1=0..d)]:
    var:={seq(seq(a[i1,j1], j1=0..d), i1=0..d)}:
    #Value conditions
    eq:={subs(x=Lx[1], S[1]) = 0, subs(x=Lx[d+2], S[d+1]) = 0,
        subs(x=Lx[i+1], S[i]) = 1}:
    #Match conditions
    eq:=eq union {seq(seq(subs(x=Lx[i1], diff2(S[i1], x, j1)) =
        subs(x=Lx[i1], diff2(S[i1-1], x, j1)), i1=2..d+1),
        j1=0..d-1)}:
    #End conditions
    eq:=eq union {seq(subs(x=Lx[1], diff2(S[1], x, j1)) = 0,
        j1=0..d-1)}:
    eq:=eq union {seq(subs(x=Lx[d+2], diff2(S[d+1], x, j1)) = 0,
        j1=0..d-1)}:
    
    var:=solve(eq, var):
    return subs(var, S):
end:

##PROBLEM 4##

#TestFS(c,Lx,x): a typlical linear combination of the Si
TestFS:=proc(c, Lx, x) local i, S, Si:
    S:=BSpline(3, [-2, -1, 0, 1, 2], x, 2);
    Si:=proc(i) local j:
        return [seq(subs(x=((x-Lx[i])/(Lx[i+1]-Lx[i])), S[j]),
            j=1..nops(S))]:
    end:
    return expand(add(c[i] * Si(Lx, x, i), i=1..nops(Lx)-2)), {seq(c[i], i=1..nops(Lx)-2)}:
end:

#RR3(p,q,f,x,Lx): analog of RR(p,q,f,x,Lx), where the phii(x)
#are replaced by those defined on p. 597 of BF, using the cubic
#B-spline, S(x)
RR3:=proc(p,q,f,x,Lx) local c, Lu, F, eq, var:
    Lu:=TestFS(c,Lx,x):
    var:=Lu[2]:
    Lu:=Lu[1]:
    F:=EvalFunct(p,q,f,x,Lx,Lu):
    eq:={seq(diff(F, var[i]), i=1..nops(var))}:
    var:=solve(eq, var):
    if var = NULL then
        return FAIL:
    fi:
    return subs(var, Lu):
end: