Mathematica Commands for Multiplying, Conjugating and Finding Resultants of Partial Differential Operators

Alex Kasman

December 2009 (modified Oct 2012)

DEFINITIONS: Partial Differential Operator: Multiplication and Conjugation

maxdimen=10;(*maximum number of variables*)DXp[vect_]:=Module[{i},Product[DX[i]^(vect[[i]]),{i,1,Length[vect]}]]

DXv[alpha_]:=Module[{i},Table[DX[i],{i,1,Length[alpha]}].alpha]

pdocoef[L_,vect_]:=Module[{xx},If[Union[vect]=={0},L/.DX[i_]->0,Coefficient[Collect[((Expand[L]/.DXp[vect]->xx)/.DX[i_]->0),xx],xx]]]

ptm[coef_,vect_,M_]:=Module[{i,sf},For[i=1;sf=M,i<Length[vect]+1,i=i+1,sf=ppm[i,vect[[i]],sf]];
coef sf]

ptm[0,vect_,M_]:=0

ppm2[i_,1,M_]:=D[M,X[i]]+M DX[i]
ppm[i_,0,M_]:=M

ppm[i_,n_,M_]:=(ppm2[i,n,M])
ppm2[i_,n_,M_]:=ppm2[i,n-1,ppm2[i,1,M]]

mtable[M_,dimen_]:=Module[{max,i,j,k,sf},For[i=1,i≤dimen+1,i=i+1,max[i]=Exponent[Collect[M,DX[i]],DX[i]]];
For[i=1;sf=Table[j[k],{k,1,dimen}],i<dimen+1,i=i+1,sf=Table[sf,{j[i],0,max[i]}]];
Flatten[sf,dimen-1]]

pdoexp[ll_,n_]:=Module[{outL,j},(For[j=0;outL=1,j<n,j=j+1,outL=pdomult[ll,outL]];outL)]


pdomult[L_,M_]:=Module[{i,j,LM,lst,dimen},dimen=findn[{L,M}];
lst=mtable[L,dimen];
For[i=1;LM=0,i≤Length[lst],i=i+1,LM=LM+ptm[pdocoef[L,lst[[i]]],lst[[i]],M]];
pdosimp[LM]]

pdotermapply[coef_,vect_,M_]:=Module[{i,sf},For[i=1;sf=M,i<Length[vect]+1,i=i+1,sf=pdopowapply[i,vect[[i]],sf]];
coef sf]

pdotermapply[0,vect_,M_]:=0

pdopowapply[i_,n_,M_]:=(D[M,{X[i],n}])

pdoapply[L_,M_]:=Module[{i,j,LM,lst,dimen},dimen=findn[L];lst=mtable[L,dimen];
For[i=1;LM=0,i≤Length[lst],i=i+1,LM=LM+pdotermapply[pdocoef[L,lst[[i]]],lst[[i]],M]];
Simplify[LM]]

pdosimp[0]:=0

pdosimp[L_]:=Module[{lst,i,sf,dimen},dimen=findn[L];Print["Simplifying..."];
lst=mtable[L,dimen];
For[i=1;sf=0,i<Length[lst]+1,i=i+1,If[MatchQ[ToString[pdocoef[L,lst[[i]]]],"0"],Print["."],Print["Simplifying ",lst[[i]]," term."];
sf=sf+Simplify[pdocoef[L,lst[[i]]]] DXp[lst[[i]]]]];
sf]

pdoconj[P_,K_]:=Module[{L,lco,LK,KP,i,j,lst,lst2,deg,jj,slv,hh,dimen},dimen=findn[{P,K}];
lst=mtable[P,dimen];
L=Sum[lco[lst[[i]]] DXp[lst[[i]]],{i,1,Length[lst]}];
deg=pdodeg[P];
For[i=1,i<Length[lst]+1,i=i+1,If[Sum[lst[[i]][[j]],{j,1,dimen}]>deg-1,jj[lst[[i]]]=0;
L=L/.lco[lst[[i]]]->pdocoef[P,lst[[i]]],jj[lst[[i]]]=100]];
KP=pdomult[K,P];
LK=pdomult[L,K];
lst2=mtable[LK,dimen];
For[i=Length[lst],i>0,i=i-1,If[jj[lst[[i]]]≠0,For[j=Length[lst2],j>0,j=j-1,hh=pdocoef[LK,lst2[[j]]];
If[StringMatchQ[ToString[D[hh,lco[lst[[i]]]]],"0"],Print["."],slv=Solve[pdocoef[LK,lst2[[j]]]==pdocoef[KP,lst2[[j]]],lco[lst[[i]]]][[1]];
j=-5;
L=L/.slv;
LK=LK/.slv]]]];
L=pdosimp[L]]



pdodeg[L_]:=Module[{transL,i,j,deg,const,xx},transL=L/.DX[i_]->const[i] xx;
deg=Exponent[Collect[transL,xx],xx];
deg]

findn[L_]:=Max[Union[{1},Table[Which[Length[Position[L,DX[i]]]>0,i,True,0],{i,1,maxdimen}]]]

(* added Oct 2012 for Robert Kragler *)

mtable[M_,K_,dimen_]:=Module[{max,i,j,k,sf},For[i=1,i≤dimen+1,i=i+1,max[i]=Exponent[Collect[M,DX[i]],DX[i]]-Exponent[Collect[K,DX[i]],DX[i]]];
For[i=1;sf=Table[j[k],{k,1,dimen}],i<dimen+1,i=i+1,sf=Table[sf,{j[i],0,max[i]}]];
Flatten[sf,dimen-1]]

pdofactor[P_,K_]:=Module[{L,lco,LK,KP,i,j,lst,lst2,deg,jj,slv,hh,dimen},dimen=findn[{P,K}];
lst=mtable[P,K,dimen];
L=Sum[lco[lst[[i]]] DXp[lst[[i]]],{i,1,Length[lst]}];
deg=pdodeg[P];
For[i=1,i<Length[lst]+1,i=i+1,jj[lst[[i]]]=100];
LK=pdomult[L,K];
lst2=mtable[LK,dimen];
For[i=Length[lst],i>0,i=i-1,If[jj[lst[[i]]]≠0,For[j=Length[lst2],j>0,j=j-1,hh=pdocoef[LK,lst2[[j]]];
If[StringMatchQ[ToString[D[hh,lco[lst[[i]]]]],"0"],Print["."],slv=Solve[pdocoef[LK,lst2[[j]]]==pdocoef[P,lst2[[j]]],lco[lst[[i]]]][[1]];
j=-5;
L=L/.slv;
LK=LK/.slv]]]];
L=pdosimp[L]]

DEFINITIONS: Resultants

Omega[d_]:=Module[{i,out,currdeg,currvar,wlist},wlist=Tuples[Table[i,{i,0,d}],numvars];
i=1;
out={};
While[i≤Length[wlist],If[Total[wlist[[i]]]≤d,out=Append[out,wlist[[i]]]];
i=i+1];
out]

mmono[omega_]:=Module[{i},Product[DX[i]^omega[[i]],{i,1,Length[omega]}]]

getdeg[operator_]:=Module[{temp,coef},Exponent[Collect[operator/.DX[i_]->coef[i] temp,temp],temp]]

Cf[operator_,term_]:=Coefficient[operator,term]/.DX[i_]->0
Cf[operator_,1]:=operator/.DX[i_]->0

vectorize[operator_,deg_]:=Module[{i,j,omega},omega=Omega[deg];
Table[Cf[operator,mmono[omega[[i]]]],{i,1,Length[omega]}]]

unvectorize[vec_,deg_]:=Module[{omega,i},omega=Omega[deg];
vec.Table[mmono[omega[[i]]],{i,1,Length[omega]}]]

predecessor[omegalist_,i_]:=Module[{j,k,elt,search,found,test},elt=omegalist[[i]];
k=0;
found=-1;
While[found<0,k=k+1;
search=elt;
search[[k]]=Max[search[[k]]-1,0];
test=Position[omegalist,search][[1,1]];
If[test<i,found=test]];
{found,k}]

compose[operator_,i_]:=Expand[Simplify[D[operator,X[i]]+operator DX[i]]]

mrows[bigN_,operator_]:=Module[{out,vec,omegalist,i,l,pred,rem},l=getdeg[operator];
Print["making rows for ",operator];
out={vectorize[operator,bigN]};
i=2;
omegalist=Omega[bigN-l];
While[i≤Length[omegalist],pred=predecessor[omegalist,i];
rem=vectorize[compose[unvectorize[out[[pred[[1]]]],bigN],pred[[2]]],bigN];
out=Append[out,rem];
i=i+1];
out]

mrows2[bigN_,operator_,ii_]:=Module[{out,tempout,mcol,vec,omegalist,i,l,pred,rem},l=getdeg[operator];
Print["making rows for ",operator];
mcol={};
out={vectorize[operator,bigN]};
mcol=Append[mcol,OP[ii]];
i=2;
omegalist=Omega[bigN-l];
While[i≤Length[omegalist],mcol=Append[mcol,OP[ii] mmono[omegalist[[i]]]];
pred=predecessor[omegalist,i];
rem=vectorize[compose[unvectorize[out[[pred[[1]]]],bigN],pred[[2]]],bigN];
out=Append[out,rem];
i=i+1];
tempout=Transpose[out];
tempout[[1]]=mcol;
Transpose[tempout]]

Oresultantmatrix[operators_,rows_]:=Module[{i,j,kbigN,rowlist,rem,M,out},numvars=Length[operators]-1;
out={};
bigN=Sum[getdeg[operators[[i]]],{i,1,numvars+1}]-numvars;
Print["bigN=",bigN];
M=Flatten[Table[mrows2[bigN,operators[[i]]-mu[i],i],{i,1,numvars+1}],1];
Print[MatrixForm[M]];
rem=Simplify[Det[Table[M[[rows[[i]]]],{i,1,Length[rows]}]]];
Table[pdosimp[rem/.{OP[i]->1,OP[m_]->0}],{i,1,numvars+1}]]

pdosimp[L_,n_]:=Module[{liszt,i,sofar},Print["Simplifying..."];
liszt=maketablefor[L,n];
For[i=1;sofar=0,i<Length[liszt]+1,i=i+1,If[MatchQ[ToString[pdocoef[L,liszt[[i]]]],"0"],Print["."],verbose["Simplifying ",liszt[[i]]," term."];
sofar=sofar+Simplify[pdocoef[L,liszt[[i]]]] DXpower[liszt[[i]]]]];
sofar]

maketablefor[M_,dimen_]:=Module[{max,i,j,k,sofar},For[i=1,i≤dimen+1,i=i+1,max[i]=Exponent[Collect[M,DX[i]],DX[i]]];
For[i=1;sofar=Table[j[k],{k,1,dimen}],i<dimen+1,i=i+1,sofar=Table[sofar,{j[i],0,max[i]}]];
Flatten[sofar,dimen-1]]

pdocoef[L_,zerovect_]:=L/.DX[i_]->0/;Union[zerovect]=={0}

pdocoef[L_,vect_]:=Module[{rm},Coefficient[Collect[((Expand[L]/.DXpower[vect]->rm)/.DX[i_]->0),rm],rm]]

DXpower[vect_]:=Module[{i},Product[DX[i]^(vect[[i]]),{i,1,Length[vect]}]]

resultantmatrix[operators_]:=Module[{i,j,bigN,rowlist,rem,M,out},numvars=Length[operators]-1;
out={};
bigN=Sum[getdeg[operators[[i]]],{i,1,numvars+1}]-numvars;
M=Flatten[Table[mrows[bigN,operators[[i]]-mu[i]],{i,1,numvars+1}],1]]

resultantmatrix[operators_,rows_]:=Module[{i,j,kbigN,rowlist,rem,M,out},numvars=Length[operators]-1;
out={};
bigN=Sum[getdeg[operators[[i]]],{i,1,numvars+1}]-numvars;
M=Flatten[Table[mrows[bigN,operators[[i]]-mu[i]],{i,1,numvars+1}],1];
Simplify[Det[Table[M[[rows[[i]]]],{i,1,Length[M[[1]]]}]]]]

Examples

Be sure the execute the section of definitions above before starting.

The notation for differential operators is that X[i] is a variable and DX[i] is the corresponding derivative.

L=DX[1]+X[2] DX[2]
Q=DX[2]-X[1]+X[3] DX[3]

the pdomult[] command multiplies the operators as "composition":

LQ=pdomult[L,Q]

QL=pdomult[Q,L]

The pdosimp[] command simplifies an operator and in particular reveals that LQ minus QL is not zero (composition is not commutative):

pdosimp[LQ-QL]

pdoapply[] applies an operator to a function

pdoapply[L,Exp[X[1] K[1]+X[2]K[2]]]

If KoL factors as KoL=QoK then pdoconj[L,K] returns Q = K o L o K^-1

pdoconj[DX[4]^2+X[2]DX[2],DX[4]-1/X[4]]

If L factors as L=QoK then pdofactor[L,K] returns Q :

prod=pdomult[Cos[X[1]]DX[1]+DX[2],Sin[X[2]]DX[1]+X[1]^2DX[2]];

pdofactor[prod,Sin[X[2]]DX[1]+X[1]^2DX[2]]

The same notation of operators is used when computing resultants:

L1=DX[1]^2-2/X[1]^2+DX[2]^2-2/X[2]^2
L2=DX[1]^3-3/X[1]^2 DX[1]+3/X[1]^3
L3=DX[2]^3-3/X[2]^2 DX[2]+3/X[2]^3

resultantmatrix[{L1,L2,L3}]

resultantmatrix[{L1,L2,L3},{1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,26,27,28,29,31,32}]