#Nathan Fox #Homework 23 #I give permission for this work to be posted online #Read procedures from class read(`C23.txt`): Help:=proc(): print(` Pw(Lf,x,Lx) , Compare(p,q,v,x,x1,L:=[10,20,30,40]) , BSpline(d,Lx,x,i) , TestFS(c,Lx,x) , RR3(p,q,f,x,Lx) `): end: ##PROBLEM 1## #Pw(Lf,x,Lx): inputs a piecwise function in our format, where Lf #is a list of expressions in the variable x, and Lx is a list of #increasing numbers, and outputs it in Maple's notation. #Anything out of bounds should be 0 #For example #Pw([x,x^2],x,[0,1,2]) #should return #piecewise(x<0, 0, x>=0 and x<1, x, 1<=x and x<2, x^2, x>=2, 0) Pw:=proc(Lf, x, Lx) local i: if nops(Lf) + 1 <> nops(Lx) then return FAIL: fi: return piecewise(x=Lx[i] and x=Lx[nops(Lx)], 0): end: ##PROBLEM 2## #Compare(p,q,v,x,x1,L): 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 ranging over L, and evluates them at x1. The output is a #list of length nops(L)+1 with the last entry being the #exact value. #NOTE: L defaults to [10,20,30,40] Compare:=proc(p, q, v, x, x1, L:=[10,20,30,40]) local u, f, N, i: u:=v*x*(1-x): f:=-diff(p*diff(u, x), x) + q*u: return [seq(subs(x=x1, RR(p, q, f, x, [seq(i/N, i=0..N)])[floor(x1*N)+1]), N in L), subs(x=x1, u)]: end: #evalf(Compare(x^2,1+x^3,sin(x),x,.513)); returned [0.123687, 0.122767, 0.122639, 0.12265444, 0.1226153759] #evalf(Compare(x^2,1+x^3,x^7+1,x,.513)); returned [0.2540307694, 0.2524677026, 0.2522576258, 0.2522113977, 0.2521669733] #evalf(Compare(x^2+3,sin(x),exp(x),x,.542)); took too long ##PROBLEM 3## #BSpline(d,Lx,i): inputs a positive 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 xxd, and continuous to order d-1 #(d+1)*(d+1) unknown numbers, #(d+2)*d equations #and f(Lx[i+1])=1 #NOTE: x needed to be an argument, so I made it one BSpline:=proc(d, Lx, x, i) local a, S, i1, j1, var, eq, diff2: diff2:=proc(u, x, n) if n > 0 then return diff(u, x\$n): else return u: fi: end: if nops(Lx) <> d+2 then return FAIL: fi: #Degree d polynomials S:=[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)}: #Value conditions eq:={subs(x=Lx[1], S[1]) = 0, subs(x=Lx[d+2], S[d+1]) = 0, subs(x=Lx[i+1], S[i]) = 1}: #Match conditions eq:=eq union {seq(seq(subs(x=Lx[i1], diff2(S[i1], x, j1)) = subs(x=Lx[i1], diff2(S[i1-1], x, j1)), i1=2..d+1), j1=0..d-1)}: #End conditions eq:=eq union {seq(subs(x=Lx[1], diff2(S[1], x, j1)) = 0, j1=0..d-1)}: eq:=eq union {seq(subs(x=Lx[d+2], diff2(S[d+1], x, j1)) = 0, j1=0..d-1)}: var:=solve(eq, var): return subs(var, S): end: ##PROBLEM 4## #TestFS(c,Lx,x): a typlical linear combination of the Si TestFS:=proc(c, Lx, x) local i, S, Si: S:=BSpline(3, [-2, -1, 0, 1, 2], x, 2); Si:=proc(i) local j: return [seq(subs(x=((x-Lx[i])/(Lx[i+1]-Lx[i])), S[j]), j=1..nops(S))]: end: return expand(add(c[i] * Si(Lx, x, i), i=1..nops(Lx)-2)), {seq(c[i], i=1..nops(Lx)-2)}: end: #RR3(p,q,f,x,Lx): analog of RR(p,q,f,x,Lx), where the phii(x) #are replaced by those defined on p. 597 of BF, using the cubic #B-spline, S(x) RR3:=proc(p,q,f,x,Lx) local c, Lu, F, eq, var: Lu:=TestFS(c,Lx,x): var:=Lu[2]: Lu:=Lu[1]: F:=EvalFunct(p,q,f,x,Lx,Lu): eq:={seq(diff(F, var[i]), i=1..nops(var))}: var:=solve(eq, var): if var = NULL then return FAIL: fi: return subs(var, Lu): end: