#hw23.txt.  Katie McKeon.  April 18, 2014.

read(`C22.txt`):
read(`C23.txt`):
read(`hw22.txt`):


Help23:=proc()
print(`Pw(Lf,x,Lx)  Compare(p,q,v,x,x1) BSpline(d,Lx,i) RR3(p,q,f,x,Lx)`): 
end:


#Pw(Lf,x,Lx) inputs a list of functions Lf in variable x
#and a partition Lx and converts to the piecewise function
#in maple
Pw:=proc(Lf,x,Lx) local i,L:

L:=[x<Lx[1],0]:
for i from 1 to nops(Lf) do
	L:=[op(L),x>=Lx[i] and x<Lx[i+1], Lf[i]]:
od:
L:=[op(L),x>Lx[-1],0]:

piecewise(op(L)):

end:


#Compare(p,q,v,x,x1) inputs expressions p,q, v in x, and a number x1 
between #0 and 1 and outputs the exact solution to the boundary value problem 
#-(p(x)*y'(x))'+ q(x)*y(x)=f(x) , 0<=x<=1, y(0)=0, y(1)=0
#followed by the four approximations given by RR with Lx=[seq(i/N,i=0..N)]
#for N=10,20,30,40
Compare:=proc(p,q,v,x,x1) local u, f, i,N,Lx,Lf,L:

u:=v*x*(1-x):
f:=-diff(p*diff(u,x),x)+q*u:

L:=[evalf(subs(x=x1, u))]:

for N from 10 while N<=40 do
	Lx:=[seq(i/N,i=0..N)]:
	Lf:=RR(p,q,f,x,Lx):
	L:=[op(L),Eval1(Lx,Lf,x,x1)]:	

	N:=N+10:
od:

L:

end:



#Data
#Compare(x^2,1+x^3,sin(x),x,.513);
#	[.1226153759, .1236858405, .1227883722, .1226664385]
#Compare(x^2,1+x^3,x^7+1,x,.513);
#	[.2521669733, .2540307694, .2524619719, .2522423878]
#Compare(x^2+3,sin(x),exp(x),x,.542); 
#	[.4268274813, .4233556228, .4260753215, .4265211453]





#BSpline(d,Lx,i), inputs a pos. 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 xxd, and continuous to order d-1 (d+1)*(d+1) unknown
#numbers, (d+2)*d equations and f(Lx[i+1])=1 
BSpline:=proc(d,Lx,i) local a,Lu,i1,j1,var,eq:

if nops(Lx)<>d+2 then
  RETURN(FAIL):
fi:

Lu:=[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)}:

eq:={subs(x=Lx[1],Lu[1]),
seq(subs(x=Lx[1],diff(Lu[1],x$i1)),i1=1..d-1)}:

eq:= eq union {subs(x=Lx[d+2],Lu[d+1]),
seq(subs(x=Lx[d+2],diff(Lu[d+1],x$i1)),i1=1..d-1)}:

eq:= eq union {seq(subs(x=Lx[j1+1],Lu[j1+1])=subs(x=Lx[j1+1],Lu[j1]),j1=1..d)}:

eq:= eq union {seq(seq(subs(x=Lx[j1+1],diff(Lu[j1+1],x$i1))=
		subs(x=Lx[j1+1],diff(Lu[j1],x$i1)),j1=1..d),
		i1=1..d-1)}:

eq:= eq union {subs(x=Lx[i+1],Lu[i+1])=1}:

var:=solve(eq,var):
if var = NULL then
	RETURN(FAIL):
fi:

subs(var, Lu):

end:










###########In Progress###############


#RR3(p,q,f,x,Lx): implements the Rayleigh-Ritz method
#using cubic splines rather than piecewise functions
#to appx. the solution of the BVP ODE
#-(p(x)*y'(x))' + q(x)*y(x)=f(x), 0<=x<=1, u(0)=0, u(1)=0
RR3:=proc(p,q,f,x,Lx) local c, Lu, S,h,n, F,Lu1, eq,var:
      
S:=BSpline(3,[-2,-1,0,1,2],2);
S:=Pw(S,x,[-2,-1,0,2,2]):   
n:=nops(Lx)-2:
h:=1/(n+1):

Lu:=TestF2(S,n,h,c,Lx,x):

    var:=Lu[2]:
    Lu:=Lu[1]:
    
   F:=EvalFunct(p,q,f,x,Lx,Lu):

    eq:=evalf({seq(diff(F,var[i]),i=1..nops(var))}):
    
    var:=solve(eq,var):

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



end:

#TestF2(S,n,h,c,Lx,x) returns a linear combination of the cubic spline basis
TestF:=proc(c,Lx,x) local i:
   expand(add(c[i]*phii2(S,n,h,Lx,x,i),i=1..nops(Lx)-1)), 
      { seq( c[i], i=1..nops(Lx)-1)}: 
   end:

   phii2:=proc(S,n,h,Lx,x,i) local j :
   if i>=2 and i<=n-1 then
	[0$(i-2),seq(subs(x=x/h,S[j]),j=1..4),0$(nops(Lx)-i-3)]:
   elif i=0 then
	[S[3]-4*subs(x=x/h, 0$nops(Lx)-3]:
   elif i=1 then

   elif i=n then

   elif i=n+1 then
	[0$nops(Lx)-3,]:
   fi:



end: