#hw23.txt; Due April 20, 2014; Frank Wagner
Help:=proc(): print(` Pw(Lf,x,Lx), Compare(p,q,v,x,x1) `):
print(` BSpline(d,Lx,i), RR3(p,q,f,x,Lx) `): end:


###################
#####Problem 1#####
###################
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(Lx)-1),
          x>=Lx[nops(Lx)],0):

end:




###################
#####Problem 2#####
###################

##old things from c23.txt##

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

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


##end old stuff from c23.txt##



Compare:=proc(p,q,v,x,x1) local u,f,y,i,N,A,B,C,D,ans:

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

A:=Eval1([seq(i/10,i=0..10)],RR(p,q,f,x,[seq(i/10,i=0..10)]),x,x1):
B:=Eval1([seq(i/20,i=0..20)],RR(p,q,f,x,[seq(i/20,i=0..20)]),x,x1):
C:=Eval1([seq(i/30,i=0..30)],RR(p,q,f,x,[seq(i/30,i=0..30)]),x,x1):
D:=Eval1([seq(i/40,i=0..40)],RR(p,q,f,x,[seq(i/40,i=0..40)]),x,x1):
ans:=evalf(subs(x=x1,u)):

[A,B,C,D,ans]:

end:


#Compare(x^2,1+x^3,sin(x),x,.513); returns
##[.1236866956, .1227627904, .1226467519, .1227246810, .1226153759]
#
#Compare(x^2,1+x^3,x^7+1,x,.513); returns
##[.2540307694, .2524677030, .2522576259, .2522113970, .2521669733]
#
#Compare(x^2+3,sin(x),exp(x),x,.542); returns
##[.4233556174, .4264487959, .4265383022, .4266448089, .4268274813]



##################
####Problem 3#####
##################

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


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

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

var:=solve(eq,var):

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

end:


#BSpline(3,[-2,-1,0,1,2],2); returns
#the same piecewise function as on page
#page 596 of the packet.



###################
#####Problem 4#####
###################

RR3:=proc(p,q,f,x,Lx) local S,Sj,n,h,i,phi,Lu,Lu1,var,F,eq:

S:=BSpline(3,[-2,-1,0,1,2],2):
n:=nops(Lx)-2:
h:=1/(n+1):
for i from 0 to n+1 do
 Sj[i]:=subs(x=(x-Lx[i+1])/h,S):
od:

phi[0]:=Sj[0]-4*subs(x=(x+h)/h,S):
phi[1]:=Sj[1]-subs(x=(x+h)/h,S):
for i from 2 to n-1 do
 phi[i]:=Sj[i]:
od:
phi[n]:=Sj[n]-subs(x=(x-(n+2)*h)/h,S):
phi[n+1]:=Sj[n+1]-4*subs(x=(x-(n+2)*h)/h,S):

Lu:=add(c[i]*phi[i],i=0..n+1):

var:={ seq( c[i], i=0..n+1)}:

Lu1:=diff(Lu,x):

F:=p*Lu1^2+q*Lu^2-2*f*Lu:

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

var:=solve(eq,var):

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

end:



#I'm not sure why this does not work,
#but it returns "division by zero"
#with everything I tried