#hw23.txt Cole Franks 20/04/2014

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

##### Problem 1

Pw:=proc(Lf,x,Lx) local f, i:

f:=piecewise(Lx[1] <x and x <= Lx[2], Lf[1]):

for i from 3 to nops(Lx) do 

f:=piecewise(Lx[1] < x and x <= Lx[i-1], f, Lx[i-1] < x and x <=Lx[i], Lf[i-1]):



od:

convert(f, piecewise,x):


end:


####### Problem 2


Compare:=proc(p,q,v,x,x1) local f, i, u, N:

u:=v*(1-x)*x:

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

#print(f):

N:=10:

#RR(p,q,f,x,[seq(i/N,i=0..N)]):

[
	seq(
		eval(Pw(RR(p,q,f,x,[seq(i/N,i=0..N)]),x, [seq(i/N,i=0..N)]), x=x1), 
	N in {seq(10*j, j=1..4)}), 
	evalf(subs(x = x1, u))
]:









end:

#Compare(x^2,1+x^3,sin(x),x,.513);
#Compare(x^2,1+x^3,x^7+1,x,.513);
#Compare(x^2+3,sin(x),exp(x),x,.542);

 #  [0.1236847910, 0.1227645773, 0.1226418161, 0.1226996267, 0.1226153759]
#   [0.2540307694, 0.2524677030, 0.2522576259, 0.2522113970, 0.2521669733]
#   [0.4233556169, 0.4264488311, 0.4265383890, 0.4266447243, 0.4268274813]



####### problem 3

BSpline:=proc(d,Lx, i) local a, Lu, i1, var, eq, j:

if d <> nops(Lx) - 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:={seq(subs(x = Lx[1], diff(Lu[1], [x$i1])), i1=0..d-1)}:

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

#loop for the spots

for j from 2 to d+1 do 

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

od:

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

var:=solve(eq, var):

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

fi:

end:

##### Problem 4


splii:=proc(Lx,x,i) local f, c, Lshift:

	if 2 < i and i <= nops(Lx) - 2 then 
		
		Pw(BSpline(3, [seq(Lx[j], j =i-2..i+2)], 2), x, [seq(Lx[j], j =i-2..i	+2)]):

	elif i = 1 then 

		#Lshift is the new list of points for this spline

		Lshift:=[-Lx[3], -Lx[2], Lx[1], Lx[2], Lx[3]]:
		
		f:=Pw(BSpline(3, Lshift, 2), x, Lshift):

		#then subtract off the right amount to kill it at 0

		c:=eval(f, x=0)/eval(f, x=Lx[2]):

		convert(f - c*subs(x = x+Lx[2], f), piecewise, x):

	elif i = 2 then 
		
		Lshift:=[-Lx[2], Lx[1], Lx[2], Lx[3], Lx[4]]:

		f:=Pw(BSpline(3, Lshift, 2), x, Lshift):

		c:=eval(f, x=0)/eval(f, x=Lx[3]):

		convert(f - c*subs(x = x+Lx[3], f), piecewise, x):

	elif i = nops(Lx) then 

		Lshift:=[Lx[-3], Lx[-2], Lx[-1], 2*Lx[-1] - Lx[-2], 2*Lx[-1] - Lx[-3]]:

		#print(Lshift):

		f:=Pw(BSpline(3, Lshift, 2), x, Lshift):

		c:=eval(f, x=Lx[-1])/eval(f, x=Lx[-2]):

		convert(f - c*subs(x = x-(Lx[-1]-Lx[-2]), f), piecewise, x):

	elif i = nops(Lx)-1 then 

		Lshift:=[Lx[-4], Lx[-3], Lx[-2], Lx[-1], 2*Lx[-1] - Lx[-2]]:

		

		f:=Pw(BSpline(3, Lshift, 2), x, Lshift):

		c:=eval(f, x=Lx[-1])/eval(f, x=Lx[-3]):

		convert(f - c*subs(x = x-(Lx[-1]-Lx[-3]), f), piecewise, x):

	fi:

end:

TestF3:=proc(c,Lx,x) local i, f:

	f:=convert(add(c[i]*splii(Lx,x,i),i=1..nops(Lx)), piecewise, x):

	
     convert(piecewise(x<0, 0, x>=0 and x < Lx[-1], f, x >= Lx[-1] ,0), piecewise, x), 
      { seq( c[i], i=1..nops(Lx))}:
     end:


EvalFunct3:=proc(p,q,f,x,Lu) local F,Lu1,i:

	#now Lu is a piecewise function
    
   Lu1:=diff(Lu,x):

   F:=expand(p*Lu1*Lu1+ q*Lu*Lu-2*f*Lu):
   
   int(F, x=0..1):
  

   end:

RR3:=proc(p,q,f,x,Lx) local c, Lu, F, eq,var, y:
   
    Lu:=TestF3(c,Lx,x):

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


    eq:=evalf({seq(diff(F,y),y in var)}):


    
    var:=solve(eq,var):

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