#C21.txt, April 10, 2014 Help:=proc() :print(`DisOp1(Oper,D1,Sten,Sh,h,ord)`): print(`DisOp(Oper,D1,Sten,Sh,h)`): print(`BVP3(Oper,D1,f,x,h1,Ini,Fini,x1)`): print(`BVP5(Oper,D1,f,x,h1,Ini,Fini,x1)`): print(`PDisOp1(Oper,D1,D2,Sten,Shx,Shy,h,ord)`): print(`PDisOp(Oper,D1,D2,Sten,Shx,Shy,h, MaxOrd)`): end: #old stuff #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: #end old stuff #[0,h1,2*h1,3*h1, i*h1=x1,,,, N*h1] #BVP3(Oper,D1,f,x,h1,x1): inputs #(i) Oper: a second-order const. coeff. diff. oper #phrased in terms of the diff. oper. D1 #(ii) an expression in the cont. variable x #(iii) Sh, a symbol for the shift operator u(x)->u(x+h) #(h is a symbol) #(iv): a small NUMERIC mesh-size h1, #(v): x1 a number between 0 and 1 ( a multiple of h1) #OUTPUTS #an appx. to u(x1) where u(x) is a solution BVP #Oper(D1)u(x)=f(x), f(0)=Ini, f(1)=Fini #using the Central 3-point discretation of Oper #with mesh size h1 BVP3:=proc(Oper,D1,f,x,h1,Ini,Fini,x1) local T,i,eq,var,N,Oper1,i1,h,Sh: N:=1/h1: if degree(Oper,D1)<>2 then RETURN(FAIL): fi: var:={seq(T[i],i=0..N)}: eq:={T[0]=Ini,T[N]=Fini}: Oper1:=DisOp(Oper,D1,{-1,0,1},Sh,h)[1]: for i from 1 to N-1 do eq:=eq union {subs(h=h1,add(coeff(Oper1,Sh,i1)*T[i+i1],i1=-1..1))= subs(x=i*h1,f)}: od: var:=solve(eq,var): if var=NULL then RETURN(FAIL): fi: subs(var,T[x1/h1]): end: #BVP5(Oper,D1,f,x,h1,x1): inputs #(i) Oper: a second-order const. coeff. diff. oper #phrased in terms of the diff. oper. D1 #(ii) an expression in the cont. variable x #(iii) Sh, a symbol for the shift operator u(x)->u(x+h) #(h is a symbol) #(iv): a small NUMERIC mesh-size h1, #(v): x1 a number between 0 and 1 ( a multiple of h1) #OUTPUTS #an appx. to u(x1) where u(x) is a solution BVP #Oper(D1)u(x)=f(x), f(0)=Ini, f(1)=Fini #using the Central 5-point discretation of Oper #with mesh size h1 BVP5:=proc(Oper,D1,f,x,h1,Ini,Fini,x1) local T,i,eq,var,N,Oper1,i1,h,Sh: N:=1/h1: if degree(Oper,D1)<>2 then RETURN(FAIL): fi: var:={seq(T[i],i=0..N)}: eq:={T[0]=Ini,T[N]=Fini}: Oper1:=DisOp(Oper,D1,{-2,-1,0,1,2},Sh,h)[1]: for i from 2 to N-2 do eq:=eq union {subs(h=h1,add(coeff(Oper1,Sh,i1)*T[i+i1],i1=-2..2))= subs(x=i*h1,f)}: od: for i from 1 to 1 do Oper1:=DisOp(Oper,D1,{-1,0,1,2,3},Sh,h)[1]: eq:=eq union {subs(h=h1,add(coeff(Oper1,Sh,i1)*T[i+i1],i1=-1..3))= subs(x=i*h1,f)}: od: for i from N-1 to N-1 do Oper1:=DisOp(Oper,D1,{-3,-2,-1,0,1},Sh,h)[1]: eq:=eq union {subs(h=h1,add(coeff(Oper1,Sh,i1)*T[i+i1],i1=-3..1))= subs(x=i*h1,f)}: od: var:=solve(eq,var): if var=NULL then RETURN(FAIL): fi: subs(var,T[x1/h1]): end: ###PDO #D1(f)=diff(f,x), Shx(f): f(x,y)->f(x+h,y) #D2(f)=diff(f,y), Shy(f): f(x,y)->f(x,y+h) #PDisOp1(Oper,D1,D2,Sten,Shx,Shy,h,ord): #inputs a partial differential operator with constant coeff. #a set of 2D lattice points, Sten, symbols for the shift operator #Shx u(x,y)->u(x+h,y), and Shy #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 #PDisOp1(D1^2+D2^2,D1,D2,{[0,0],[-1,0],[1,0],[0,-1],[0,1]}, # Shx,Shy, h,3); should give #(1/Shx+1/Shy+Shx+Shy-4)/h^2 PDisOp1:=proc(Oper,D1,D2,Sten,Shx,Shy,h,ord) local eq,var,a,i, j,Oper1, Oper1a,var1,i1: option remember: Oper1:=add(a[i]*Shx^i[1]*Shy^i[2], i in Sten): var:={seq(a[i], i in Sten)}: Oper1a:=subs({Shx=exp(h*D1), Shy=exp(h*D2)}, Oper1-Oper): Oper1a:=taylor(Oper1a,D1=0,ord+1): Oper1a:=add(coeff(Oper1a,D1,i1)*D1^i1,i1=0..ord): Oper1a:=taylor(Oper1a,D2=0,ord+1): Oper1a:=add(coeff(Oper1a,D2,i1)*D2^i1,i1=0..ord): eq:={seq(seq(coeff(coeff(Oper1a,D1,i),D2,j),j=0..ord-i), i=0..ord)}: var:=solve(eq,var): if var=NULL then RETURN(FAIL): fi: var1:={}: for i from 1 to nops(var) do if op(1,op(i,var))=op(2,op(i,var)) then var1:=var1 union {op(2,op(i,var))}: fi: od: Oper1:=subs(var,Oper1): Oper1:=subs({seq(var1[i]=0,i=1..nops(var1))},Oper1): if Oper1=0 then RETURN(FAIL): else RETURN(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 PDisOp:=proc(Oper,D1,D2,Sten,Shx,Shy,h,MaxOrd) local ord,try1,ord1: for ord from 1 to MaxOrd do try1:= PDisOp1(Oper,D1,D2,Sten,Shx,Shy,h,ord): if try1<>FAIL then for ord1 from ord to MaxOrd while PDisOp1(Oper,D1,D2,Sten,Shx,Shy,h,ord1)<>FAIL do od: RETURN( PDisOp1(Oper,D1,D2,Sten,Shx,Shy,h,ord1-1)): fi: od: FAIL: end: