#April 17, 2014, continuing Rayleigh-Ritz
#C23.txt
Help:=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) , BSpline(d,Lx,i)`):
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:

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

end: