#hw23.txt, 4/19/14, Anthony Zaleski

Help:=proc():
HelpC():
print(`Pw(Lf,x,Lx),Compare(p,q,v,x,x1), BSpline(d,Lx,i),RR3(p,q,f,x,Lx),phii3:=proc(Lx,x,i),TestF3(c,Lx,x),PlotPiecewise(Lx,S,x,Ly:=0)`):
end:

##########Problem 1: Converting piecewise formats##########

#Pw(Lf,x,Lx): Inputs a piecewise function in our
#notation, outputs it using Maple's notation.

Pw:=proc(Lf,x,Lx) local i:
piecewise(x<Lx[1],0,seq(op([x>=Lx[i] and x<Lx[i+1],Lf[i]]),i=1..nops(Lf)),x>=Lx[-1],0):
end:

##########Problem 2: Testing Rayleigh-Ritz##########

#Compare(p,q,v,x,x1): 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=10,20,30,40, and evlautes them at x1. 
#The output is a list of length 5 with the last entry being the exact value.

Compare:=proc(p,q,v,x,x1) local u,f,L,N:
u:=v*x*(1-x):
f:= -diff((p*diff(u,x)),x)+ q*u: #print(f);
[seq(Eval1([seq(i/10/N,i=0..10*N)],RR(p,q,f,x,[seq(i/10/N,i=0..10*N)]),x,x1),N=1..4), evalf(subs(x=x1,u))]:
end:

###Output:
#Compare(x^2,1+x^3,sin(x),x,.513);
#[.1236849191, .1227644544, .1226438106, .1226878035, .1226153759]
#Compare(x^2,1+x^3,x^7+1,x,.513);
#[.2540307694, .2524677030, .2522576259, .2522113970, .2521669733]
#Compare(x^2+3,sin(x),exp(x),x,.542);
#[.4233556116, .4264488324, .4265384158, .4266447313, .4268274813]


##########Problem 3: BSpline##########

#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 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
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)}:

#Continuity at the interior points:
eq:=eq union {seq(subs(x=Lx[i1],Lu[i1-1])=subs(x=Lx[i1],Lu[i1]),i1=2..d+1)}:

#Continuity of derivatives at interior points:
for j1 from 1 to d-1 do
	eq:=eq union {seq(subs(x=Lx[i1],diff(Lu[i1-1],x$j1))=subs(x=Lx[i1],diff(Lu[i1],x$j1)),i1=2..d+1)}:
od:

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

var:=solve(eq,var):

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

subs(var,Lu):

end:

#Test:
#Lu:=BSpline(3,Lx,2): 
#Lu2:=expand([(2-x)^3-4*(1-x)^3-6*x^3+4*(1+x)^3,(2-x)^3-4*(1-x)^3-6*x^3,(2-x)^3-4*(1-x)^3,(2-x)^3]/4):
#expand(Lu-Lu2);
#[0, 0, 0, 0]


##########Problem 4: Rayleigh-Ritz with phi's composed of splines##########


#First, we change phii and TestF:
phii3:=proc(Lx,x,i) local L0,Lx0,S0,h,Lx1,j,S1,a,b, Lxi,Si:
if i=0 then
	L0:=Lx[2..3]-[Lx[1]$2]:
	Lx0:=[Lx[1]$5]+[-L0[2],-L0[1],0,L0[1],L0[2]]:
	S0:=BSpline(3,Lx0,2)[3..4]:

	h:=L0[1]: 
	Lx1:=[Lx[1]$5]+[seq(-h+j*h,j=-2..2)]:
	S1:=BSpline(3,Lx1,2)[4]:
	a:=subs(x=Lx[1],S0[1]): b:=subs(x=Lx[1],S1):
	RETURN([S0[1]-a/b*S1,S0[2],0$nops(Lx)-3]):
fi:
if i=1 then
	L0:=Lx[2..4]-[Lx[1]$3]:
	Lx0:=[Lx[1]$5]+[-L0[1],0,L0[1],L0[2],L0[3]]:
	S0:=BSpline(3,Lx0,2)[2..4]:

	h:=L0[1]: 
	Lx1:=[Lx[1]$5]+[seq(-h+j*h,j=-2..2)]:
	S1:=BSpline(3,Lx1,2)[4]:
	a:=subs(x=Lx[1],S0[1]): b:=subs(x=Lx[1],S1):
	RETURN([S0[1]-a/b*S1,S0[2],S0[3],0$nops(Lx)-4]):
fi:
if i=nops(Lx)-1 then
	L0:=Lx[-3..-2]-[Lx[-1]$2]:
	Lx0:=[Lx[-1]$5]+[L0[1],L0[2],0,-L0[1],-L0[2]]:
	S0:=BSpline(3,Lx0,2)[1..2]:

	h:=L0[2]: 
	Lx1:=[Lx[-1]$5]+[seq(-h-j*h,j=-2..2)]:
	S1:=BSpline(3,Lx1,2)[1]:
	a:=subs(x=Lx[-1],S0[2]): b:=subs(x=Lx[-1],S1):
	RETURN([0$nops(Lx)-3,S0[1],S0[2]-a/b*S1]):

fi:
if i=nops(Lx)-2 then
	L0:=Lx[-4..-2]-[Lx[-1]$3]:
	Lx0:=[Lx[-1]$5]+[L0[1],L0[2],L0[3],0,-L0[3]]:
	S0:=BSpline(3,Lx0,2)[1..3]:

	h:=L0[3]: 
	Lx1:=[Lx[-1]$5]+[seq(-h-j*h,j=-2..2)]:
	S1:=BSpline(3,Lx1,2)[1]:
	a:=subs(x=Lx[-1],S0[3]): b:=subs(x=Lx[-1],S1):
	RETURN([0$nops(Lx)-4,S0[1],S0[2],S0[3]-a/b*S1]):

fi:

Lxi:=Lx[i-1..i+3]:
Si:=BSpline(3,Lxi,2):

[0$i-2,op(Si),0$nops(Lx)-i-3]:

end:


#TestF3(c,Lx,x): a typical inear combination of the phi_i
TestF3:=proc(c,Lx,x) local i:
expand(add(c[i]*phii3(Lx,x,i),i=1..nops(Lx)-2)), 
     { seq( c[i], i=1..nops(Lx)-2)}:
end:

 
#################

   #RR3(p,q,f,x,Lx): implements the Rayleigh-Ritz method, using B-splines
   #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, F, eq,var:
   
    Lu:=TestF3(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:






#############################
##########OLD STUFF##########

#April 17, 2014, continuing Rayleigh-Ritz
#C23.txt
HelpC:=proc():print(` EvalFunct(p,q,f,x,Lx,Lu):`): 
print(` EvalFunct(1,0,2*Pi^2*sin(Pi*x),x,[0,1/2,1],[u1(x),u2(x)]); `):
print(` RR(p,q,f,x,Lx)`):
end:

#old stuff
 Help22:=proc(): print(` EL(L,Y,Y1,y,t), Eval1(Lx,Lf,x,x1), Mul1a(p,Lf) `): 
  print(`Mul1b(Lf,Lg) , phii(Lx,x,i), TestF(c,Lx,x)`):
  end:

  #Eval1(Lx,Lf,x,x1): Given a partition of [0,1]=[x0,x1,x2,...,xn,x_{n+1}]
  #and corresponding list of expressions Lf of n+1 expressions
  #where the function between xi and x_{i+1} is Lf[i+1]
  #EVALUATE it at x1
  Eval1:=proc(Lx,Lf,x,x1) local i:

    if x1<0 or x1>1 then
       RETURN(0):
    fi:

   for i from 1 while Lx[i]<x1 do od:
      i:=i-1:
    
   subs(x=x1,Lf[i]): 

  end:

  #Mul1a(p,Lf): multiplies the global function p(x) by the piecewise function Lf
   Mul1a:=proc(p,Lf)  local i:
    [seq(p*Lf[i],i=1..nops(Lf))]:
    end:
    
  #Mul1b(Lf,Lg)
  Mul1b:=proc(Lf,Lg) local i:
   [seq(Lf[i]*Lg[i],i=1..nops(Lf))]:   
  end:
  
   phii:=proc(Lx,x,i) :
   [0$(i-1),(x-Lx[i])/(Lx[i+1]-Lx[i]),(Lx[i+2]-x)/(Lx[i+2]-Lx[i+1]), 0$(nops(Lx)-i-2)]:
   end:


#end old stuff

   #TestF(c,Lx,x): a typical inear combination of the phi_i
    TestF:=proc(c,Lx,x) local i:
     expand(add(c[i]*phii(Lx,x,i),i=1..nops(Lx)-2)), 
      { seq( c[i], i=1..nops(Lx)-2)}:
     end:

   #EvalFunct(p,q,f,x,Lx,Lu): inputs expressions p,q,f
   #in the variable x, a partition of [0,1], Lx, and
   #a piece-wise function Lu, in x
   # 
   #outputs the value of (11.37)
   #INT( p(x)*u'(x)^2+q(x)*u(x)^2-2*f(x)*u(x),x=0..1); 
   
   EvalFunct:=proc(p,q,f,x,Lx,Lu) local F,Lu1,i:
    
   Lu1:=diff(Lu,x):

   F:=expand(p*Mul1b(Lu1,Lu1)+ q*Mul1b(Lu,Lu)-2*f*Lu):
   
   add(int(F[i],x=Lx[i]..Lx[i+1]), i=1..nops(Lx)-1):
  

   end:

   #RR(p,q,f,x,Lx): implements the Rayleigh-Ritz method
   #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
   RR:=proc(p,q,f,x,Lx) local c, Lu, F, eq,var:
   
    Lu:=TestF(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:

#PlotPiecewise:=proc(Lx,S,x,Ly:=0): Inputs array Lx and array S (eg a spline)
#of expressions in x.  Plots the piecewise-defined S.  If Ly is inputted, also 
#plots the discrete [Lx,Ly] points (useful for checking an interpolation)
PlotPiecewise:=proc(Lx,S,x,Ly:=0) local i:
if Ly=0 then
	plots[multiple](plot,seq([S[i],x=Lx[i]..Lx[i+1]],i=1..nops(S))):
else 
	plots[multiple](plot,seq([S[i],x=Lx[i]..Lx[i+1]],i=1..nops(S)),[Lx,Ly,style=point,symbolsize=20]):
fi:
end:

##########END OLD##########