#C20a.txt: getting better and better discrete versions #of differential operators, April 7, 2014 Help:=proc() :print(`DisOp1(Oper,D1,Sten,Sh,h,ord)`): print(`DisOp(Oper,D1,Sten,Sh,h)`):end: #D1(f)=diff(f,x), Sh(f): f(x)->f(x+h) #DisOp1(Oper,D1,Sten,Sh,h,ord): #inputs a differential operator with constant coeff. #a set of integers, Sten, a symbol for the shift operator #Sh u(x)->u(x+h), a symbol for h, the mesh size #and a pos. intetger ord #Outputs: a recurrence operator (a Laurent poly in Sh) #whose powers are in the Sten #For example #DisOp1(D1,D1,{0,1},Sh,h,1); should give #(Sh-1)/h, #DisOp1(D1^2,D1,{-1,0,1},Sh,h,2); should give #(Sh-2+1/Sh)/h^2 DisOp1:=proc(Oper,D1,Sten,Sh,h,ord) local eq,var,a,i, Oper1, Oper1a: option remember: Oper1:=add(a[i]*Sh^i, i in Sten): var:={seq(a[i], i in Sten)}: Oper1a:=subs(Sh=exp(h*D1),Oper1)-Oper: Oper1a:=taylor(Oper1a,D1=0,ord+1): eq:={seq(coeff(Oper1a,D1,i),i=0..ord)}: var:=solve(eq,var): if var=NULL then RETURN(FAIL): else RETURN(subs(var,Oper1)): fi: end: #inputs a differential operator with constant coeff. #a set of integers, Sten, a symbol for the shift operator #Sh u(x)->u(x+h), a symbol for h, the mesh size #Outputs: a recurrence operator (a Laurent poly in Sh) #whose powers are in the Sten, with the highest possible #order, along with that order #For example #DisOp(D1,D1,{0,1},Sh,h); should give #(Sh-1)/h, #DisOp1(D1^2,D1,{-1,0,1},Sh,h,2); should give #(Sh-2+1/Sh)/h^2 DisOp:=proc(Oper,D1,Sten,Sh,h) local ord,Oper1: for ord from 0 while DisOp1(Oper,D1,Sten,Sh,h,ord)<>FAIL do od: DisOp1(Oper,D1,Sten,Sh,h,ord-1),ord-1: end: #BVP(Oper,D1,f,x, a,b, Ini,Fini,h1,winS,x1): #uses discrete versions to order ord of #the diff. oper. with constant coeff. Oper #phrased in terms of D1 (diff. oper) #and a function f of x, and a number h1 (mesh size) #considers the BVP #Oper(D1) u(x)=f(x), f(a)=Ini[1], f'(a)=Ini[2], ... #f(b)=Fini[1], f'(b)=Fini[2], ... #output an apprx. to the number u(x1) BVP:=proc(Oper,D1,f,x,a,b,Ini,Fini,h1,winS,x1) local Oper1,Sten: if degree(Oper,D1)<>nops(Ini)+nops(Fini) then RETURN(FAIL): fi: end: