#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[i] and 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]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