#hw23.txt. Katie McKeon. April 18, 2014. read(`C22.txt`): read(`C23.txt`): read(`hw22.txt`): Help23:=proc() print(`Pw(Lf,x,Lx) Compare(p,q,v,x,x1) BSpline(d,Lx,i) RR3(p,q,f,x,Lx)`): end: #Pw(Lf,x,Lx) inputs a list of functions Lf in variable x #and a partition Lx and converts to the piecewise function #in maple Pw:=proc(Lf,x,Lx) local i,L: L:=[x=Lx[i] and xLx[-1],0]: piecewise(op(L)): end: #Compare(p,q,v,x,x1) inputs expressions p,q, v in x, and a number x1 between #0 and 1 and outputs the exact solution to the boundary value problem #-(p(x)*y'(x))'+ q(x)*y(x)=f(x) , 0<=x<=1, y(0)=0, y(1)=0 #followed by the four approximations given by RR with Lx=[seq(i/N,i=0..N)] #for N=10,20,30,40 Compare:=proc(p,q,v,x,x1) local u, f, i,N,Lx,Lf,L: u:=v*x*(1-x): f:=-diff(p*diff(u,x),x)+q*u: L:=[evalf(subs(x=x1, u))]: for N from 10 while N<=40 do Lx:=[seq(i/N,i=0..N)]: Lf:=RR(p,q,f,x,Lx): L:=[op(L),Eval1(Lx,Lf,x,x1)]: N:=N+10: od: L: end: #Data #Compare(x^2,1+x^3,sin(x),x,.513); # [.1226153759, .1236858405, .1227883722, .1226664385] #Compare(x^2,1+x^3,x^7+1,x,.513); # [.2521669733, .2540307694, .2524619719, .2522423878] #Compare(x^2+3,sin(x),exp(x),x,.542); # [.4268274813, .4233556228, .4260753215, .4265211453] #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 xxd, 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[j1+1],Lu[j1+1])=subs(x=Lx[j1+1],Lu[j1]),j1=1..d)}: eq:= eq union {seq(seq(subs(x=Lx[j1+1],diff(Lu[j1+1],x\$i1))= subs(x=Lx[j1+1],diff(Lu[j1],x\$i1)),j1=1..d), i1=1..d-1)}: eq:= eq union {subs(x=Lx[i+1],Lu[i+1])=1}: var:=solve(eq,var): if var = NULL then RETURN(FAIL): fi: subs(var, Lu): end: ###########In Progress############### #RR3(p,q,f,x,Lx): implements the Rayleigh-Ritz method #using cubic splines rather than piecewise functions #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 RR3:=proc(p,q,f,x,Lx) local c, Lu, S,h,n, F,Lu1, eq,var: S:=BSpline(3,[-2,-1,0,1,2],2); S:=Pw(S,x,[-2,-1,0,2,2]): n:=nops(Lx)-2: h:=1/(n+1): Lu:=TestF2(S,n,h,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: #TestF2(S,n,h,c,Lx,x) returns a linear combination of the cubic spline basis TestF:=proc(c,Lx,x) local i: expand(add(c[i]*phii2(S,n,h,Lx,x,i),i=1..nops(Lx)-1)), { seq( c[i], i=1..nops(Lx)-1)}: end: phii2:=proc(S,n,h,Lx,x,i) local j : if i>=2 and i<=n-1 then [0\$(i-2),seq(subs(x=x/h,S[j]),j=1..4),0\$(nops(Lx)-i-3)]: elif i=0 then [S[3]-4*subs(x=x/h, 0\$nops(Lx)-3]: elif i=1 then elif i=n then elif i=n+1 then [0\$nops(Lx)-3,]: fi: end: