##################################################################################################### # Algorithms for Linear Partial Differential Operators. Focus on Darboux and Invariant theory. # By Ekaterina Shemyakova. #################################################################################################################### # # Highlights of the package #################################################################################################################### # # ginv: ginv(L) returns a set of generating differential invariants for L with respect to the gauge transformations # ginv(L,M) returns a set of generating joint differential invariants for L and M with respect to the gauge transformations # # Laplace_solver # Solve by Laplace method by < or = n interations. # L is either operator in the form D_xy + a D_x + b D_y + c or # L is the corresponding PDE u_xy + a u_x + b u_y + c u = 0 # one may need to spesify the third argument func, the dependable variable with respect to which is our PDE ##################################################################################################### # 1. Basic operations: # # LPDO__set_vars( vars::list ) set a list of variables. # # LPDO__create( ) create an LPDO of the degree zero (that is all the coefficients are zeros). # LPDO__create0(a00) create an LPDO of degree 0 with the free coefficient equaled to a00. # LPDO__create1:= proc(a10,a01,a00) create LPDO a10*Dx + a01*Dy + a00. # LPDO__create2:= proc(a20,a11,a02,a10,a01,a00) create LPDO a20*Dxx + a11*Dxy + a02*Dyy+a10*Dx + a01*Dy + a00. # Create_LPDO_from_coeff:= proc(a::list) given list "a" of the coefficients, creates an LPDO. # # LPDO__degree(L) find degree of LPDO L. # LPDO__Symbol(L) computes the symbol of L in variables X and Y, bivariate case only. # LPDO__simplify(a) simplify LPDO coefficients and decrease degree. # LPDO__expand(a) expand LPDO coefficients and decrease degree # LPDO__factor(a) factor LPDO coefficients and decrease degree # LPDO__print(a) pretty print. # LPDO__add(L1::array,L2::array) add LPDOs L1 and L2. # LPDO__inc(L1::array,L2::array) add LPDOs L1 and L2 and save the result as L1 (L1 += L2). # LPDO__set_value(L::array,f,lst::list) set the coefficient at D_{lst} in LPDO L equaled to f. # LPDO__add_value(L::array,f,lst::list) add to the coefficient at D_{lst} in LPDO L value f. # LPDO__value(L::array,lst::list) returns the coefficient at D_lst in LPDO L. # LPDO__minus(L1,L2) subtract L2 from L1. # LPDO__mult1(val,L::array) multiply LPDO L by function val on the left: # LPDO__mult(L1,L2) returns the result of the composition of operators L1 and L2. # LPDO__conj(L::array, f) retuns (1/f)*L*f. # LPDO__apply(L::array, f) apply LPDO L to function f. # LPDO__transpose(L::array) returns the formal transpose of LPDO L (a*D1^k1*D2^k2*D3^k3 -> (-D1)^k1*(-D2)^k2*(-D3)^k3 *a). ####################################################################################################### # 2. Invariants of LPDOs. # # LPDO__Lapl_h(L) Laplace h-invariant for L=Dxy+... and Dxx+Dyy+..., i.e. h s.t. L=(Dx+b)(Dy+a)+h # LPDO__Lapl_k(L) Laplace k-invariant for L=Dxy+... and Dxx+Dyy+..., i.e. k s.t. L=(Dy+a)(Dx+b)+k # LPDO__Lapl_m(L) h-k # LPDO__Inv_strict_parabolic_ord2(n,L) Invariants for operators of the form L=Dxx + aDx + b Dy + c , b <>0 # LPDO__Lapl_gen_h(p1,q1,p2,q2,a,b,c) invariant h for hyperbolic LPDOs Sym(L)= (p1 X + q1 Y)(p2 X + q2 Y) # LPDO__Lapl_gen_k(p1,q1,p2,q2,a,b,c) invariant k for hyperbolic LPDOs Sym(L)= (p1 X + q1 Y)(p2 X + q2 Y) # LPDO__Inv(n,a) computes the invariant number "n" for the operator "a" with the symbol XY(pX+qY), # n=1,..,7, and pXXX n=1..5 (4)(3) and pXXY n=1..5 # LPDO__Inv_Transposed(n,L) computes the invariant # number "n" for the operator "L^t" # with the symbol XY(pX+qY), XXY, XXX # in terms of the invariants of L # For the XXX case a02=0 and a11<>0 is denoted by I5=0 ######################################################################################################### # 3. Factorization property in terms of invariants. # # LPDO__obstacle(r1,s1,r2,s2,r3,s3,a::array) computes common obstacle to factorization of the type (S1)(S2)(S3). # LPDO__factor111(r1,s1,r2,s2,r3,s3,a::array) Si= ri X + si Y for the operator "a" with the symbol S1*S2*S3. # # factorizations of HYPERBOLIC LPDOs, that is with the symbol XYS, where S=pX+qY # LPDO__Is_Factorable(p,q,i1,i2,i3,i4,i5) outputs a set of sets of conditions s.t. input operator is factorable if and only if one of the sets of conditions holds # LPDO__Is_Factorable_via_Coeff(a) # LPDO__Fact_Cond_S_XY:= proc(p,q,I1,I2,I3,I4,I5) # LPDO__Fact_Cond_S_X_Y:= proc(p,q,i1,i2,i3,i4,i5) # LPDO__Fact_Cond_S_Y_X:= proc(p,q,i1,i2,i3,i4,i5) # LPDO__Fact_Cond_Y_X_S:= proc(p,q,i1,i2,i3,i4,i5) # LPDO__Fact_Cond_X_Y_S:= proc(p,q,i1,i2,i3,i4,i5) # LPDO__Fact_Cond_X_S_Y:= proc(p,q,i1,i2,i3,i4,i5) # LPDO__Fact_Cond_Y_S_X:= proc(p,q,i1,i2,i3,i4,i5) # LPDO__Fact_Cond_X_SY:= proc(p,q,i1,i2,i3,i4,i5) # LPDO__Fact_Cond_Y_SX:= proc(p,q,i1,i2,i3,i4,i5) # LPDO__Fact_Cond_SY_X:= proc(p,q,i1,i2,i3,i4,i5) # LPDO__Fact_Cond_XY_S:= proc(p,q,i1,i2,i3,i4,i5) # LPDO__Fact_Cond_SX_Y:= proc(p,q,i1,i2,i3,i4,i5) # # factorizations of NON-HYPERBOLIC LPDOs with the symbol XXY. In the rest of the cases we dont know sufficient conditions # LPDO__Fact_Cond_X_XY:= proc(i1,i2,i3,i4,i5) # LPDO__Fact_Cond_XY_X:= proc(i1,i2,i3,i4,i5) case i2=i3 works provided r[x] = i4+(1/4)*a11^2+(1/2)*a11[x]-r*a11+r^2 can be solved for r # LPDO__Fact_Cond_Y_XX:= proc(i1,i2,i3,i4,i5) # LPDO__Fact_Cond_XX_Y:= proc(i1,i2,i3,i4,i5) iff conditions # LPDO__Fact_Cond_X_X_Y:= proc(i1,i2,i3,i4,i5) iff conditions provided r1[x] = -i4-(1/4)*a11^2+(1/2)*a11[x]+r1*a11-r1^2 has solution w.r.t. r1. # # LPDO__Fact_Cond_X_XX_a02_0:= proc(i1,i2,i3,i4) # LPDO__Fact_Cond_Var2_Var1Var1:= proc(I1,I2,I3,I4,I5) Necessary and sufficient invariant conditions of # LPDOs of the eq. class defined by the values I1, ... ,I5 of the # generating set of invaraints to be factorable with the factorization type (Var2)(Var1^2) ######################################################################### # 4. Laplace transformations. # # Factorization_solver:=proc(L) solves PDEs u_xy+au_x + bu_y +cu=0 corresponding to factorizable L. Guarantees general solution. # Laplace_solver(n, L) returns solution of L by means of Laplace transformations method, n - max number of interations considered. # Laplace_solver(n,L,dependent_variable) # Solve by Laplace method by < or = n interations. # L is either operator in the form D_xy + a D_x + b D_y + c or # L is the corresponding PDE u_xy + a u_x + b u_y + c u = 0 # one may need to spesify the third argument, the dependable variable with respect to which is our PDE # LPDO__Lapl_trans(i,L) returns Li for i=+/-1 L=Dxy+... or L=Dxx + Dyy. # LPDO__Lapl_chain(n,L) for L=Dxy+... or L=Dxx + Dyy outputs [[k,h],[k_1,h1],...,[k_n,hn]] # ######################################################################### # 5. Darboux transformations theory. # # LPDO__Lapl_trans(i,L::array) returns Li for i=+/-1 L=Dxy+... or L=Dxx + Dyy. # LPDO__Lapl_chain(n,L) for L=Dxy+... or L=Dxx + Dyy outputs [[k,h],[k_1,h1],...,[k_n,hn]] # # Darboux := proc(L,z1,z2) by determinant formulae returns [M,L1,M1] # DarbouxX := proc(L,z1) returns [M,L1,M1], uses SOLVE -> can be a problem # DarbouxY := proc(L,z1) returns [M,L1,M1], uses SOLVE -> can be a problem # # LPDO__pl := proc(L, M) --- projection from the left, returns M-A*L without DxDy # LPDO__pr := proc(L, M) --- projection from the right, returns M-L*A without DxDy # # LPDO__pairs_invariants(L,M) finds pairs invariants w.r.t. gauge transforamtiosn of the pairs: m,h,q, R and w.r.t. to gauged evolution: I1, I2, and I3. # LPDO__DT_conditions_in_gauged_evolution_invariants:=proc(q,I2,I3) returns conditions for the existence of a Darboux transformation in terms of given invariants I1,I2,I3 of pair M, L # w.r.t. gauged evolution. # LPDO__pairs_invariants_evolution_via_gauged:=proc(m,h,q,R) finds pairs invariants w.r.t. to gauged evolution: I1, I2, and I3 in terms of # invariants of pairs w.r.t. gauge transformations of the pairs: m,h,q, R and # LPDO__joint_gauge_invariants:=proc(L::array,M::array) returns joint gauge invariants for operators L=D_xy + aD_x +bD_y +c and M=m_d D_x^d + ... m_{-d} D_y^d # M must be without any mixed derivatives, returns [m,h,[R_1,R_2, ... R_d],[R_{-1}, .., R_{-d}]] # LPDO__pairs_invariants3(L::array,M::array) returns the gauged evolution invariants for pair (L,M), where L=D_x D_x D_y + .. and M=D_x+m or M=D_y+m # ######################################################################################## # 6. Invertible Darboux transformations # # LPDO__XXY_DT_with_Dy:=proc(b,I2,I3,I4,r) # LPDO__XXY_DT_with_Dx:=proc(i1,i2,i3,i4,i5,r) case i2<>i3 only # LPDO__XYS_DT_with_Dx:=proc(i1,i2,i3,i4,i5,op::integer:=0) # LPDO__XYS_DT_with_Dy:=proc(i1,i2,i3,i4,i5) # LPDO__XYS_DT_with_Dx_Dy:=proc(i1,i2,i3,i4,i5) # # LPDO__a_step_of_Eq_and_Shift:=proc(L,M,N,L1,Minv,Ninv,A,G,S,T) performs a step of Equivalence-Shift procedure for L of arbitrarily form. Returns [Ln,Mn,Nn,L1n,Minvn,Ninvn,An,Gn]. # LPDO__a_step_of_Eq_and_Shift_basic:=proc(L,M,N,L1,S,T) performs a step of Equivalence-Shift procedure M -> M+S*L, N -> N+L1*S for L of arbitrarily form, # does not compute the auxilary operators, just the change in L and M, L1, N. Returns [Ln,Mn,Nn,L1n] # LPDO__data_for_invertible_DT_of_Type_I:=proc(C,M,f) returns full data for invertible DT for L=CM+f with symbol XXY: [L,M,N,L1,Minv,Ninv,A,G] # LPDO__a_step_of_Eq_and_Shift_for_Type_I:=proc(C,M,f,S,T) performs a step of Equivalence-Shift procedure for L with symbol XXY for type I so for L=CM+f # LPDO__DT_of_type_I:=proc(C,M,f) performs a step of invertible Darboux transformation of type I. Returns [N,L1,Minv,Ninv,A,G]. # LPDO__XYS_chain_of_first_order_invertible_DTs:=proc(chain,i1,i2,i3,i4,i5) computes the invariants in the requested branch of the chain. The branch is input by argument "chain" # as an array of numbers 1,2, and 3. E.g. [1,2,1,3,3,1] means that first one should apply a DT with X, then with Y, then with X, then with X+Y symbols. # Returns the list of the invariants. ####Internal################ LPDO__XYS_Chart_Plot_step_of_chain_of_first_order_invertible_DTs:=proc(a,i1,i2,i3,i4,i5) ######################## LPDO__XYS_Chart_Plot_step_of_chain_of_first_order_invertible_DTs:=proc(a,i1,i2,i3,i4,i5) ######################## LPDO__Charts_next_level:=proc(res) ######################## LPDO__clean_chart:=proc(l) # LPDO__XYZ_Chart:=proc(b,list_of_invs) computes (including plotting) b levels of DTs chart. Does not trace cycles. Returns list of lists: [L1,L2,L3] # L1: command display(L1, scaling=constrained,axes=Normal); will plot the chart # L2: list of points (D_y up, D_x right, D_x+D_y on the diagonal) # L3: list of lits of invariants (in correspondence with L2) ######################################################################################## # 7. Internal functions # # LPDO__nvars() find the number of the declared variables. # LPDO__size(L) number of monoms of an operator L. # LPDO__LPDO(a) convert from old LPDO. # LPDO__nmonoms(deg) Find the number of monomials of orders <= deg # LPDO__deg_0(nvar, nmo) # LPDO__create_by_degree(deg::integer);create an operator of the degree deg (all the coeff.=0) # LPDO__set_degree(a::array,deg::integer) set degree(a) >= deg # LPDO__multDi(m::integer, i::integer, k::integer); Let D[[m]] - diff.monom of index m. add D[i]^k to list D[[m]]; # LPDO__multDiii(f, m::integer, a::array); f*D[[m]]*a # LPDO__transpose1(am, i); # LPDO__Inv_XXX(n,a) # compute the invariant number "n" for an operator with the symbol pXXX, # n=1,2,3,4,5 if a02 neq 0 # n=1,..,4, if a02=0 # n=1,..,3, if a02=a11=0 # LPDO__Inv_XYS := proc(n::posint,a::array) compute the invariant number "n" for an operator with the symbol XY(pX+qY), n=1,..,7 # DIFF(f,lis::list) DIFF(f,[i1, ..., i_n])= D_{x1}^{i1} ... D_{xn}^{in} (f) # Subs_list() returns substitutions needed for computation of generating sets of invariants my means of moving frames # LPDO__Lapl_chain_step1(k,h) (k,h) -> (k_1,h1) for L=Dxy+... # LPDO__Lapl_chain_classic(n, L) for L=Dxy+... outputs [[k,h],[k_1,h1],...,[k_n,hn]] # LPDO__Lapl_gen_form(n, L) for L=Dxx + Dyy... outputs [[k,h],[k_1,h1],...,[k_n,hn]] # LPDO__Laplace_trans_classic(i,L::array) returns L1, if i=1, and L_{-1} if i=-1. # LPDO__Laplace_trans_gen_form(i,L::array) returns Li for L=Dxx + Dyy + aDx + bDy + c ###################################################################################################### # THE CODE STARTS ###################################################################################################### with(combinat): LPDO__warnings:=false: LPDO__vars:=[x,y]: #------------------------------------------------------- # Set a list of variables LPDO__set_vars := proc( vars::list ) global LPDO__vars; LPDO__vars := vars; end: #------------------------------------------------------- # Find the number of variables LPDO__nvars := proc( ) global LPDO__vars; if type(LPDO__vars, 'list') and nops(LPDO__vars)>0 then RETURN(nops(LPDO__vars)); fi; ERROR(`illegal LPDO__vars`); end: #----------------------------------------------------------- # Find the number of monomials of orders <= deg LPDO__nmonoms:= proc(deg); vectoint([deg+1, 0$(LPDO__nvars()-1)]); end: #----------------------------------------------------------- # LPDO__deg_0:= proc(nvar, nmo) local deg, nmo1; option remember; nmo1:=-1; for deg from 0 while nmo1<=nmo do nmo1:= vectoint([deg+1, 0$(nvar-1)]); if nmo1>=nmo then RETURN(deg); fi: od; ERROR(`illegal LPDO__deg_0`); end: #----------------------------------------------------------- LPDO__create:= proc( ); array(0..0, [0]); end: #----------------------------------------------------------- LPDO__create_by_degree:= proc(deg::integer) local nmo; nmo:=LPDO__nmonoms(deg); array(0..nmo-1, [0$nmo]); end: #----------------------------------------------------------- LPDO__create0:= proc(a00) local res; res:=LPDO__create(); res[0]:=a00; eval(res); end: #----------------------------------------------------------- # number of monoms of operator: LPDO__size:= proc(a::array); ArrayTools[Size](eval(a))[2]; # op(op(2, eval(a))[1])[2]+1; # that is the top bound of first dimension of array end: #---------------------------------------------------------------------- # find degree of operator: LPDO__degree:= proc(a::array); LPDO__deg_0(LPDO__nvars(), LPDO__size(a) ); end: #---------------------------------------------------------------------- # set degree(a) >= deg # # if degree(a) >= deg, then return source a; # if degree(a) < deg, expand array and return # LPDO__set_degree:= proc(a::array, deg::integer) local nmo, sz, res, i; nmo:=LPDO__nmonoms(deg); sz :=LPDO__size(a); if nmo<=sz then RETURN(a); fi; res:=array(0..nmo-1, [0$nmo]); for i from 0 to sz-1 do res[i]:=a[i]; od; eval(res); end: #---------------------------------------------------------------------- # simplify LPDO coefficients LPDO__simplify := proc(a::array) local i, im, deg, res; im:=-1; for i from LPDO__size(a)-1 by -1 while i>=0 do a[i]:=simplify(eval(a[i])); # a[i]:=expand(simplify(eval(a[i]))); if (a[i]<>0) and (im=-1) then im:=i; fi; od: # im - index of last, <>0, value; for deg from 0 do if LPDO__nmonoms(deg)-1 >= im then break; fi; od: res := LPDO__create_by_degree(deg); for i from 0 to im do res[i]:=a[i]; od; eval(res); end: #---------------------------------------------------------------------- # expand(simplify) LPDO coefficients LPDO__expand := proc(a::array) local i, im, deg, res; im:=-1; for i from LPDO__size(a)-1 by -1 while i>=0 do # a[i]:=simplify(eval(a[i])); a[i]:=expand(simplify(eval(a[i]))); if (a[i]<>0) and (im=-1) then im:=i; fi; od: # im - index of last, <>0, value; for deg from 0 do if LPDO__nmonoms(deg)-1 >= im then break; fi; od: res := LPDO__create_by_degree(deg); for i from 0 to im do res[i]:=a[i]; od; eval(res); end: #---------------------------------------------------------------------- # factor LPDO coefficients LPDO__factor := proc(a::array) local i, im, deg, res; im:=-1; for i from LPDO__size(a)-1 by -1 while i>=0 do # a[i]:=simplify(eval(a[i])); a[i]:=factor(eval(a[i])); if (a[i]<>0) and (im=-1) then im:=i; fi; od: # im - index of last, <>0, value; for deg from 0 do if LPDO__nmonoms(deg)-1 >= im then break; fi; od: res := LPDO__create_by_degree(deg); for i from 0 to im do res[i]:=a[i]; od; eval(res); end: #---------------------------------------------------------------------- # print LPDO__print := proc(a::array) local i, nv; #print(`LPDO__print, size=`, LPDO__size(a)); nv := LPDO__nvars(); for i from 0 to LPDO__size(a)-1 do if a[i]<>0 then print(inttovec(i,nv), a[i] ); fi; od: end: #---------------------------------------------------------------------- # (a1+a2) LPDO__add := proc(a1::array,a2::array) local nmo1, nmo2, i, res; nmo1 := LPDO__size(a1); nmo2 := LPDO__size(a2); if nmo1>= nmo2 then res:=copy(a1); for i from 0 to nmo2-1 do res[i]:=res[i]+a2[i]; od; else # if nmo1< nmo2 then res:=copy(a2); for i from 0 to nmo1-1 do res[i]:=res[i]+a1[i]; od; fi; LPDO__simplify(res); end: #---------------------------------------------------------------------- # a1 += a2 LPDO__inc := proc(a1::array,a2::array) local nmo1, nmo2, i, res; nmo1 := LPDO__size(a1); nmo2 := LPDO__size(a2); if nmo1>= nmo2 then for i from 0 to nmo2-1 do a1[i]:=a1[i]+a2[i]; od; else # if nmo1< nmo2 then a1 := LPDO__add(a1,a2); fi; RETURN(a1); end: #---------------------------------------------------------------------- # a += f*D[list] LPDO__add_value := proc(a::array,f,lst::list) local i, deg_lst, deg_a; if LPDO__nvars() <> nops(lst) then ERROR(`in LPDO__set_value illegal number of vars`); fi; deg_a := LPDO__degree(a); deg_lst := add(i,i=lst); if deg_lst > deg_a then a:=LPDO__set_degree(a,deg_lst); fi; i:=vectoint(lst); a[i]:=a[i]+f; RETURN(a); end: #---------------------------------------------------------------------- # a := f*D[list] LPDO__set_value := proc(a::array,f,lst::list) local i, deg_lst, deg_a; if LPDO__nvars() <> nops(lst) then ERROR(`in LPDO__set_value illegal number of vars`); fi; deg_a := LPDO__degree(a); deg_lst := add(i,i=lst); if deg_lst > deg_a then a:=LPDO__set_degree(a,deg_lst); fi; i:=vectoint(lst); a[i]:=f; RETURN(a); end: #---------------------------------------------------------------------- # a[ D[list] ] LPDO__value := proc(a::array,lst::list) local i; if LPDO__nvars() <> nops(lst) then ERROR(`in LPDO__value illegal number of vars`); fi; i:=vectoint(lst); if i>=LPDO__size(a) then RETURN(0); fi; a[i]; end: #---------------------------------------------------------------------- # (a1-a2) LPDO__minus := proc(a1::array,a2::array) local nmo1, nmo2, i, res; nmo1 := LPDO__size(a1); nmo2 := LPDO__size(a2); if nmo1>= nmo2 then res:=copy(a1); for i from 0 to nmo2-1 do res[i]:=res[i]-a2[i]; od; else # if nmo1< nmo2 then res:=copy(a2); for i from 0 to nmo1-1 do res[i]:=res[i]-a1[i]; od; fi; LPDO__simplify(res); end: #---------------------------------------------------------------------- # multiplication a by function val (left): return(val*a); LPDO__mult1 := proc(val,a) local i, res; res:=copy(a); for i from 0 to LPDO__size(a)-1 do res[i]:=val*a[i]; od; LPDO__simplify(res); end: #---------------------------------------------------------------------- # Let D[[m]] - diff.monom of index m. # add D[i]^k to list D[[m]]; LPDO__multDi:= proc(m::integer, i::integer, k::integer) local lst; if k=0 then RETURN(m); fi; lst:= inttovec(m,LPDO__nvars()); # D1^lst[1] * D2^lst[2] *... lst[i] := lst[i]+k; vectoint(lst); end: #---------------------------------------------------------------------- # D[i]^k * a LPDO__multDii:= proc(i::integer, k::integer, a::array) local res, m, k1, ind2; res := LPDO__create_by_degree( LPDO__degree(a)+k); for m from 0 to LPDO__size(a)-1 do ind2 := LPDO__multDi( m, i, k ); # k1=0: res[ind2] := res[ind2] + a[m]; for k1 from 1 to k do # k1>0 ind2 := LPDO__multDi( m, i, k-k1 ); res[ind2] := res[ind2] + binomial(k,k1)*diff(a[m], LPDO__vars[i]$k1); od; od; LPDO__simplify(res); end: #---------------------------------------------------------------------- # f*D[[m]]*a LPDO__multDiii:= proc(f, m::integer, a::array) local nv, inv, lst, res; nv := LPDO__nvars(); # number of variables lst := inttovec(m,nv); # list of derivativies of D[[m]] res := a; for inv from nv by -1 while inv>0 do if lst[inv]>0 then res := LPDO__multDii( inv, lst[inv], res ); fi; od; LPDO__mult1(f,res); end: #---------------------------------------------------------------------- # D[i]^k * a LPDO__mult:= proc(a::array, b::array) local res, res1, m, k1, ind2; res := LPDO__create_by_degree(LPDO__degree(a)+LPDO__degree(b)); for m from 0 to LPDO__size(a)-1 do if a[m]<>0 then LPDO__inc( res,LPDO__multDiii(a[m], m, b) ); fi; od; LPDO__simplify(res); end: #---------------------------------------------------------------------- # LPDO_ <-- LPDO LPDO__LPDO:= proc(a::array) local deg, res, i1, i2; if LPDO__nvars() <> 2 then ERROR(`LPDO__LPDO: nvars <> 2`); fi; deg := LPDO_degree(a); res := LPDO__create_by_degree(deg); for i1 from 0 to deg do for i2 from 0 to deg-i1 do LPDO__set_value( res, a[i1,i2], [i1,i2] ); od: od: LPDO__simplify(res); end: #---------------------------------------------------------------------- # (1/f)*a*f LPDO__conj:= proc(a::array, f) local ff; ff := LPDO__create0(f); LPDO__mult1(1/f, LPDO__mult(a,ff)); end: #---------------------------------------------------------------------- # a(f) LPDO__apply:= proc(a::array, f) local nv, res,m, am, lst, tmp, inv; nv := LPDO__nvars(); # number of variables res := 0; for m from 0 to LPDO__size(a)-1 do am := simplify(a[m]); if am<>0 then lst := inttovec(m,nv); # list of derivativies of D[[m]], D1^lst[1] * D2^lst[2] *... tmp := f; for inv from nv by -1 while inv>0 do if lst[inv]>0 then tmp := diff( tmp, LPDO__vars[inv]$lst[inv] ); fi; od; res := res + am*tmp; fi; od: simplify(res); end: #----------------------------------------------------------- LPDO__create1:= proc(a10,a01,a00) local res; if LPDO__nvars() <> 2 then ERROR(`LPDO__create1: nvars <> 2`); fi; res:=LPDO__create0(a00); LPDO__add_value(res, a10 , [1,0] ): LPDO__add_value(res, a01 , [0,1] ): # LPDO__print(res): eval(res); end: #----------------------------------------------------------- LPDO__create2:= proc(a20,a11,a02,a10,a01,a00) local res; if LPDO__nvars() <> 2 then ERROR(`LPDO__create2: nvars <> 2`); fi; res:=LPDO__create0(a00); LPDO__add_value(res, a10 , [1,0] ): LPDO__add_value(res, a01 , [0,1] ): LPDO__add_value(res, a20 , [2,0] ): LPDO__add_value(res, a11 , [1,1] ): LPDO__add_value(res, a02 , [0,2] ): # LPDO__print(res): eval(res); end: #---------------------------------------------------------------------- # one member: a*D1^k1*D2^k2*D3^k3 -> (-D1)^k1*(-D2)^k2*(-D3)^k3 *a LPDO__transpose1:= proc(am, i) local lst, z, dg, res, resa; lst := inttovec(i,LPDO__nvars()); # list of derivativies of D[[m]], D1^lst[1] * D2^lst[2] *... res := array(0..i, [0$(i+1)]); # print(777,am,i, res); dg := add(z,z=op(lst)); # degree res[i]:=(-1)^dg; # LPDO__print(res); resa:=LPDO__create0(am); # print(778); LPDO__print(resa); res:=LPDO__mult(res, resa); # print(779); LPDO__print(res); RETURN(eval(res)); end: #---------------------------------------------------------------------- # a*D1^k1*D2^k2*D3^k3 -> (-D1)^k1*(-D2)^k2*(-D3)^k3 *a LPDO__transpose:= proc(a::array) local m, res, am; res := LPDO__create(); for m from 0 to LPDO__size(a)-1 do am := simplify(a[m]); # print(780, m, am); if am<>0 then res := LPDO__add( res, LPDO__transpose1(am,m)); fi; od: RETURN(eval(simplify(res))); end: #---------------------------------------------------------------- # # invariant h for hyperbolic LPDOs Sym(L)= (p1 X + q1 Y)(p2 X + q2 Y) # LPDO__Lapl_gen_h:= proc(p1,q1,p2,q2,a,b,c) local resu,x,y; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; x:=LPDO__vars[1]; y:=LPDO__vars[2]; if p1=0 and q1=0 then ERROR(`LPDO__Lapl_gen_h: given operator has order one`); fi; if p2=0 and q2=0 then ERROR(`LPDO__Lapl_gen_h: given operator has order one`); fi; if p1=0 and p2=0 then ERROR(`LPDO__Lapl_gen_h: given operator is non hyperbolic`); fi; if q1=0 and q2=0 then ERROR(`LPDO__Lapl_gen_h: given operator is non hyperbolic`); fi; if p2<>0 and q1=(q2*p1/p2) then ERROR(`LPDO__Lapl_gen_h: given operator is non hyperbolic`); fi; if p2<>0 then resu:=1/(p1*q2-p2*q1)^2*c*p2^2*q1^2+1/(p1*q2-p2*q1)^2*c*p1^2*q2^2+1/(p1*q2-p2*q1)^2*p1*b^2*p2+2/(p1*q2-p2*q1)^2*p1^3*diff(q2,x)^2*p2-1/(p1*q2-p2*q1)^2*p1^2*q2^2*diff(a,x)-1/(p1*q2-p2*q1)^2*p2^2*q1^2*diff(b,y)-3/(p1*q2-p2*q1)^2*a*q2*q1^2*diff(p2,y)+2/(p1*q2-p2*q1)^2*b*q2*p1^2*diff(p2,x)-2/(p1*q2-p2*q1)^2*q2*p1^3*diff(p2,x)*diff(q2,x)+2/(p1*q2-p2*q1)^2*p1^2*diff(p2,x)^2*q2*q1+1/(p1*q2-p2*q1)^2*p1^2*q2^2*diff(q1,x)*diff(p2,y)+1/(p1*q2-p2*q1)^2*p1*q2^2*a*diff(p1,x)-1/(p1*q2-p2*q1)^2*q1*q2^2*diff(a,y)*p1-1/(p1*q2-p2*q1)^2*q1^2*q2^2*diff(p2,y)*diff(p1,y)+1/(p1*q2-p2*q1)^2*q1*q2^2*a*diff(p1,y)-2/(p1*q2-p2*q1)^2*q1^3*diff(p2,y)*p2*diff(q2,y)-1/(p1*q2-p2*q1)^2*p1*p2^2*diff(b,x)*q1+1/(p1*q2-p2*q1)^2*p1*p2^2*b*diff(q1,x)+1/(p1*q2-p2*q1)^2*p1^2*p2*diff(b,x)*q2-1/(p1*q2-p2*q1)^2*p1^2*p2^2*diff(q2,x)*diff(q1,x)+1/(p1*q2-p2*q1)^2*p2^2*q1*b*diff(q1,y)+1/(p1*q2-p2*q1)^2*p2*q1^2*q2*diff(a,y)-1/(p1*q2-p2*q1)^2*a*b*p1*q2-1/(p1*q2-p2*q1)^2*p2*q1^3*q2*diff(p2,`$`(y,2))-2/(p1*q2-p2*q1)^2*p1*p2*q2*q1^2*diff(p2,x,y)-2/(p1*q2-p2* q1)^2*p1^2*p2*q1*diff(q2,x,y)*q2-1/(p1*q2-p2*q1)^2*p2*q1^2*diff(q2,`$`(y,2))*p1*q2+2/(p1*q2-p2*q1)^2*p1^2*q2^2*q1*diff(p2,x,y)+1/(p1*q2-p2*q1)^2*p2^2*q1^2*diff(p1,y)*diff(q2,x)-1/(p1*q2-p2*q1)^2*a*b*p2*q1+2/(p1*q2-p2*q1)^2*a*p2*q1^2*diff(q2,y)-3/(p1*q2-p2*q1)^2*b*p2*p1^2*diff(q2,x)+2/(p1*q2-p2*q1)^2*p1*q1^2*diff(q2,y)^2*p2+1/(p1*q2-p2*q1)^2*q1^2*diff(p2,y)*b*p2-1/(p1*q2-p2*q1)^2*p1^2*p2*q2*diff(p2,`$`(x,2))*q1+2/(p1*q2-p2*q1)^2*p1*p2^2*q1^2*diff(q2,x,y)-2/(p1*q2-p2*q1)^2*p1*diff(p2,x)*p2*q1^2*diff(q2,y)+1/(p1*q2-p2*q1)^2*p1*p2^2*diff(p1,x)*diff(q2,x)*q1-1/(p1*q2-p2*q1)^2*p1*p2*b*diff(p1,x)*q2-1/(p1*q2-p2*q1)^2*p1*p2*q2*a*diff(q1,x)-1/(p1*q2-p2*q1)^2*p1*q2^2*q1*diff(p2,y)*diff(p1,x)-1/(p1*q2-p2*q1)^2*p1^2*p2*diff(q1,x)*diff(q2,y)*q2+1/(p1*q2-p2*q1)^2*p1^2*p2*q2*diff(p2,x)*diff(q1,x)+1/(p1*q2-p2*q1)^2*p1*p2*q1*diff(q2,y)*diff(p1,x)*q2+1/(p1*q2-p2*q1)^2*p2*q1*diff(b,y)*p1*q2+1/(p1*q2-p2*q1)^2*q1*q2^2*diff(q1,y)*diff(p2,y)*p1-1/(p1*q2-p2*q1)^2*p2*q1*b*diff(p1,y)*q2+2/(p1*q2-p2*q1)^2*a*p1*diff(q2,x)*p2*q1+1/(p1* q2-p2*q1)^2*a*p1*q2*q1*diff(q2,y)-2/(p1*q2-p2*q1)^2*q2*p1^2*diff(p2,x)*q1*diff(q2,y)+4/(p1*q2-p2*q1)^2*p1^2*diff(q2,x)*p2*q1*diff(q2,y)-3/(p1*q2-p2*q1)^2*p1*b*p2*q1*diff(q2,y)+2/(p1*q2-p2*q1)^2*p1*b*q2*q1*diff(p2,y)-2/(p1*q2-p2*q1)^2*p1*q2*q1^2*diff(p2,y)*diff(q2,y)+1/(p1*q2-p2*q1)^2*p1*diff(p2,x)*b*p2*q1+4/(p1*q2-p2*q1)^2*p1*diff(p2,x)*q2*q1^2*diff(p2,y)-2/(p1*q2-p2*q1)^2*p1^2*diff(p2,x)*diff(q2,x)*p2*q1+1/(p1*q2-p2*q1)^2*a*p1^2*q2*diff(q2,x)-1/(p1*q2-p2*q1)^2*p1*p2*q2*diff(p1,x)*diff(p2,x)*q1-2/(p1*q2-p2*q1)^2*c*p2*p1*q2*q1+2/(p1*q2-p2*q1)^2*q1^3*diff(p2,y)^2*q2+1/(p1*q2-p2*q1)^2*q2*a^2*q1+1/(p1*q2-p2*q1)^2*p2^2*q1^3*diff(q2,`$`(y,2))+1/(p1*q2-p2*q1)^2*p1^3*q2^2*diff(p2,`$`(x,2))-1/(p1*q2-p2*q1)^2*p2*q1*q2*a*diff(q1,y)-1/(p1*q2-p2*q1)^2*p2^2*q1*p1*diff(q2,x)*diff(q1,y)-1/(p1*q2-p2*q1)^2*p2*q1^2*q2*diff(p1,y)*diff(p2,x)-1/(p1*q2-p2*q1)^2*p2*q1*diff(q1,y)*diff(q2,y)*p1*q2+1/(p1*q2-p2*q1)^2*p2*q1*q2*p1*diff(p2,x)*diff(q1,y)+1/(p1*q2-p2*q1)^2*p2*q1^2*diff(q2,y)*diff(p1,y)*q2-3/(p1*q2-p2*q1)^2*a*q2*p1*diff(p2,x) *q1-2/(p1*q2-p2*q1)^2*q2*q1*diff(p2,y)*p1^2*diff(q2,x)+1/(p1*q2-p2*q1)^2*p1*p2*q2*diff(a,x)*q1-2/(p1*q2-p2*q1)^2*q1^2*diff(p2,y)*p1*diff(q2,x)*p2-1/(p1*q2-p2*q1)^2*p1^3*p2*diff(q2,`$`(x,2))*q2+1/(p1*q2-p2*q1)^2*p1^2*p2^2*diff(q2,`$`(x,2))*q1+1/(p1*q2-p2*q1)^2*q1^2*q2^2*diff(p2,`$`(y,2))*p1; else resu:=(c*q2*p1^2-a*b*p1+a*diff(q2,x)*p1^2+q1*a^2+a*q1*diff(q2,y)*p1-q2*p1^2*diff(a,x)+q2*p1*a*diff(p1,x)-q1*q2*diff(a,y)*p1+q1*q2*a*diff(p1,y))/q2/p1^2; fi; RETURN(simplify(eval(resu))); end: #---------------------------------------------------------------- # # invariants k for hyperbolic LPDOs Sym(L)= (p1 X + q1 Y)(p2 X + q2 Y) LPDO__Lapl_gen_k:= proc(p1,q1,p2,q2,a,b,c) local resu,x,y; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; x:=LPDO__vars[1]; y:=LPDO__vars[2]; if p1=0 and q1=0 then ERROR(`LPDO__Lapl_gen_h: given operator has order one`); fi; if p2=0 and q2=0 then ERROR(`LPDO__Lapl_gen_h: given operator has order one`); fi; if p1=0 and p2=0 then ERROR(`LPDO__Lapl_gen_h: given operator is non hyperbolic`); fi; if q1=0 and q2=0 then ERROR(`LPDO__Lapl_gen_h: given operator is non hyperbolic`); fi; if p2<>0 and q1=(q2*p1/p2) then ERROR(`LPDO__Lapl_gen_h: given operator is non hyperbolic`); fi; if p2<>0 then resu:= (c*p1^2*q2^2+p2^3*diff(p1,`$`(x,2))*q1^2-diff(a,x)*p2^2*q1^2+p1*p2*b^2-q2*diff(p2,y)*p1*b*q1+p1*b*q2^2*diff(p1,y)-2*p1*q2^3*diff(p1,y)*diff(q1,y)+a*p2^2*q1*diff(q1,x)-2*diff(p1,x)*p2^3*q1*diff(q1,x)+2*diff(p1,x)^2*p2^2*q1*q2+2*p1*p2*q2^2*diff(q1,y)^2+2*diff(p1,x)*p2^2*b*q1-p1^2*p2*diff(b,x)*q2-q2^2*p1^2*diff(b,y)+a*p2*q1*q2*diff(q1,y)+diff(a,x)*p2*q1*p1*q2+p1^2*p2*diff(q1,x)*diff(p2,x)*q2-p1*p2^2*diff(q2,x)*diff(q1,y)*q1+p1*p2^2*diff(q2,x)*diff(p1,x)*q1-p1*p2*diff(q2,x)*a*q1-diff(p2,x)*q2*diff(p1,y)*p2*q1^2-2*p1*p2^2*diff(q1,x)*q2*diff(p1,x)-2*p1*p2*diff(q1,x)*q2^2*diff(p1,y)+4*p1*p2^2*diff(q1,x)*q2*diff(q1,y)+p1*p2*b*q2*diff(p1,x)+q2*p1*p2*diff(q2,y)*diff(p1,x)*q1+4*diff(p1,x)*p2*q1*q2^2*diff(p1,y)-2*diff(p1,x)*p2^2*q1*q2*diff(q1,y)+2*q2*diff(p1,y)*p2*b*q1-2*c*p2*q1*p1*q2+q2*diff(p1,y)*diff(q2,y)*p2*q1^2-3*p1*p2^2*b*diff(q1,x)+q2^2*diff(a,y)*q1*p1+q2*p1^2*b*diff(q2,y)-a*b*p1*q2-3*p1*p2*b*q2*diff(q1,y)-q2^2*diff(p2,y)*p1*q1*diff(p1,x)+q2^2*diff(p2,y)*p1*q1*diff(q1,y)-q2*p1*p2*diff(q2,y)*diff(q1, y)*q1-q2*p1^2*p2*diff(q1,x)*diff(q2,y)-diff(p2,x)*p1*p2*b*q1-diff(p2,x)*p1*p2*q1*q2*diff(p1,x)+diff(p2,x)*p1*p2*q1*q2*diff(q1,y)-q2*diff(a,y)*p2*q1^2+c*p2^2*q1^2+2*p1*p2^3*diff(q1,x)^2-p1*p2^2*q2*diff(p1,`$`(x,2))*q1-2*p1*p2*q2^2*diff(p1,x,y)*q1-2*p1*p2^2*q2*diff(q1,x,y)*q1-q2^2*p1*p2*diff(q1,`$`(y,2))*q1+2*q2^3*diff(p1,y)^2*q1-a*p2*b*q1+q2^2*p1^2*diff(q1,x)*diff(p2,y)-2*p1*p2*q2^2*diff(p1,x)*diff(q1,y)-q2*p1*diff(q2,y)*a*q1+q2^2*diff(p1,`$`(y,2))*p2*q1^2+p1^2*p2^2*diff(q1,`$`(x,2))*q2+2*q2*diff(p1,x,y)*p2^2*q1^2-q2^3*p1*diff(p1,`$`(y,2))*q1-p1*p2^3*diff(q1,`$`(x,2))*q1+2*p1^2*p2*q2^2*diff(q1,x,y)+q2*p1*p2*diff(b,y)*q1-2*q2*diff(p1,y)*p2^2*q1*diff(q1,x)-2*q2^2*diff(p1,y)*p2*q1*diff(q1,y)+2*a*p2*diff(q1,x)*p1*q2-3*a*p2*q1*q2*diff(p1,x)-3*a*q1*q2^2*diff(p1,y)+2*a*p1*q2^2*diff(q1,y)-diff(p2,y)*q2^2*diff(p1,y)*q1^2+q2*diff(p2,y)*a*q1^2+diff(p2,x)*a*p2*q1^2+p1*p2^2*diff(b,x)*q1+p1^2*p2*b*diff(q2,x)-p1^2*p2^2*diff(q1,x)*diff(q2,x)+diff(q2,x)*diff(p1,y)*p2^2*q1^2+a^2*q1*q2+q2^3*p1^2*diff(q1,`$`(y,2)))/(-p2*q1+p1*q2) ^2; else resu:=-(-c*q2*p1^2+b*p1*a-b*p1*q2*diff(p1,y)-q1*a^2+3*q1*a*q2*diff(p1,y)-2*q1*q2^2*diff(p1,y)^2-2*diff(q1,y)*q2*p1*a+2*diff(q1,y)*q2^2*p1*diff(p1,y)+diff(b,y)*q2*p1^2-q1*q2*p1*diff(a,y)+q1*q2^2*p1*diff(p1,`$`(y,2))-q2^2*diff(q1,`$`(y,2))*p1^2-diff(q2,y)*b*p1^2+diff(q2,y)*q1*p1*a)/q2/p1^2; fi; RETURN(simplify(eval(resu))); end: #---------------------------------------------------------------- # # Laplace invariants for L=Dxy+... and S=Dxx+Dyy # checked # # LPDO__Lapl_h := proc(L::array) local a,b,c,x,y; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; if LPDO__degree(L) <> 2 then ERROR(`illegal degree of operator`); fi; x:=LPDO__vars[1]; y:=LPDO__vars[2]; a := LPDO__value(L, [1,0]); b := LPDO__value(L, [0,1]); c := LPDO__value(L, [0,0]); if LPDO__value(L, [2,0])=0 and LPDO__value(L, [1,1])=1 and LPDO__value(L, [0,2])=0 then RETURN(expand(simplify( a*b + diff(a,x)-c))); elif LPDO__value(L, [2,0])=1 and LPDO__value(L, [1,1])=0 and LPDO__value(L, [0,2])=1 then RETURN(simplify(eval( 1/4*b^2+1/4*a^2+(1/2*I)*(diff(b, x))+(1/2)*(diff(b,y))+1/2*diff(a,x)-1/2*I*diff(a,y)-c))); else ERROR(`operator should have the symbol of the form XY or X^2 + Y^2`); fi; print(11111); end: # LPDO__Lapl_k := proc(L::array) local a,b,c,x,y; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; if LPDO__degree(L) <> 2 then ERROR(`illegal degree of operator`); fi; x:=LPDO__vars[1]; y:=LPDO__vars[2]; a := LPDO__value(L, [1,0]); b := LPDO__value(L, [0,1]); c := LPDO__value(L, [0,0]); if LPDO__value(L, [2,0])=0 and LPDO__value(L, [1,1])=1 and LPDO__value(L, [0,2])=0 then RETURN(expand(simplify( -c+a*b + diff(b,y)))); elif LPDO__value(L, [2,0])=1 and LPDO__value(L, [1,1])=0 and LPDO__value(L, [0,2])=1 then RETURN(simplify(eval(-c+1/4*b^2+1/4*a^2+1/2*diff(a,x)+1/2*I*diff(a,y)-1/2*I*diff(b,x)+diff(b,y)/2))); else ERROR(`operator should have the symbol of the form XY or X^2 + Y^2`); fi; end: #-------------------------------------------------------------------------------------------------------- # h-k # checked # # LPDO__Lapl_m := proc(L::array) local a,b,c,x,y; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; if LPDO__degree(L) <> 2 then ERROR(`illegal degree of operator`); fi; x:=LPDO__vars[1]; y:=LPDO__vars[2]; a := LPDO__value(L, [1,0]); b := LPDO__value(L, [0,1]); c := LPDO__value(L, [0,0]); if LPDO__value(L, [2,0])=0 and LPDO__value(L, [1,1])=1 and LPDO__value(L, [0,2])=0 then RETURN(expand(simplify( diff(a,x) - diff(b,y) ))); elif LPDO__value(L, [2,0])=1 and LPDO__value(L, [1,1])=0 and LPDO__value(L, [0,2])=1 then RETURN(simplify(eval(-I*(diff(a, y)-(diff(b,x)))))); else ERROR(`operator should have the symbol of the form XY or X^2 + Y^2`); fi; end: #----------------------------------------------------------- # compute the invariant number "n" for an operator with the symbol # XY(pX+qY), n=1,..,7 LPDO__Inv_XYS := proc(n::posint,a::array) local p,q,a20,a11,a02,a10,a01,a00; if n>7 then ERROR(` LPDO__Inv: illegal index of invariant `); fi; if LPDO__nvars() <> 2 then ERROR(` LPDO__Inv: illegal number of vars`); fi; if LPDO__degree(a) <> 3 then ERROR(` LPDO__Inv: illegal degree of operator`); fi; if (LPDO__value(a,[3,0])<>0) or (LPDO__value(a,[0,3])<>0) then ERROR(` LPDO__Inv: the second argument should be an operator with the symbol of the form XY(pX+qY)`); fi; p :=LPDO__value(a,[2,1]); q :=LPDO__value(a,[1,2])/p; a20:=LPDO__value(a,[2,0])/p; a11:=LPDO__value(a,[1,1])/p; a02:=LPDO__value(a,[0,2])/p; a10:=LPDO__value(a,[1,0])/p; a01:=LPDO__value(a,[0,1])/p; a00:=LPDO__value(a,[0,0])/p; if (n=1) then RETURN(simplify(eval(2*p^3*(q^2*a20-1/2*q*a11+a02)))); elif (n=2) then RETURN(simplify(eval(-p^4*(q*diff(a02,LPDO__vars[2])-diff(q,LPDO__vars[2])*a02-diff(a20,LPDO__vars[1])*q^2)))); elif (n=3) then RETURN(simplify(eval(p^3*(-a20*a11+2*a20*diff(q,LPDO__vars[2])+a10+a20^2*q-diff(a11,LPDO__vars[2])+diff(a20,LPDO__vars[2])*q)))); elif (n=4) then RETURN(simplify(eval(p^3*q^2*(a01+1/q^2*a02^2-3/q^2*a02*diff(q,LPDO__vars[1]) -diff(a11,LPDO__vars[1])+1/q*diff(a02,LPDO__vars[1]) +1/q*diff(q,LPDO__vars[1])*a11-1/q*a02*a11)))); elif (n=5) then RETURN(simplify(eval(p^5*q*( a00+1/q*a20*a11*a02-1/q*a02*a10 -a01*a20-1/2*diff(a11,LPDO__vars[1],LPDO__vars[2])+diff(q,LPDO__vars[1])*diff(a20,LPDO__vars[2]) +2*diff(a20,LPDO__vars[1])/q*a02-diff(a20,LPDO__vars[1])*a11+diff(a20,LPDO__vars[1])*diff(q,LPDO__vars[2]) +2*diff(a20,LPDO__vars[1])*q*a20+diff(q,LPDO__vars[1],LPDO__vars[2])*a20)))); elif (n=6) then RETURN(eval(p)); elif (n=7) then RETURN(eval(LPDO__value(a,[1,2]))); else RETURN(`the first argument should be 1,..,7`); fi; end: #----------------------------------------------------------------- # # # computes invariants of L^t in terms of invariants of L, n=1,2,3,4,5 # for operators with the symbol XY(pX+qY) # LPDO__Inv_Transposed := proc(n,I1,I2,I3,I4,I5, a30, a21, a12, a03) local p,q; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; if a30=0 and a03=0 and a21<>0 and a12<>0 then p:=a21; q:=a12; if n=1 then RETURN(-2*q^2*diff(p,y)+2*diff(p,x)*q*p +I1+2*diff(q,y)*q*p-2*p^2*diff(q,x)); elif n=2 then RETURN(-diff(q,x,y)*p^2*q+diff(q,x)*diff(q,y) *p^2+diff(p,x,y)*p*q^2-diff(p,x) *diff(p,y)*q^2 -I2); elif n=3 then RETURN((2*q^3*diff(p,y)^2-I3*q^2+diff(q,`$`(y,2))*q^2*p^2 +q*p*diff(I1,y)-diff(q,y)*p*I1-q^3*p*diff(p,`$`(y,2)) -2*q*diff(p,y)*I1-2*diff(p,y)*diff(q,y)*q^2*p)/q^2 +2*p*I2/q^2); elif n=4 then RETURN((diff(p,`$`(x,2))*q^2*p^2-2*diff(q,x)*I1*p +diff(I1,x)*q*p-diff(p,x)*I1*q-q*p^3*diff(q,`$`(x,2)) -2*p^2*diff(q,x)*q*diff(p,x)-I4*p^2+2*p^3*diff(q,x)^2)/p^2 -2*q*I2/p^2); elif n=5 then RETURN((1/q*p^2*diff(p,y)+2/q^2*diff(q,y)*p^3)*I4 +(p*diff(q,x)+2*q*diff(p,x))*I3+(4/q*diff(q,x)*diff(p,y)*p +4/q^2*diff(q,x)*diff(q,y)*p^2-diff(p,x,y)*p +1/q*diff(p,x)*diff(q,y)*p-2/q*diff(q,x,y)*p^2 +3*diff(p,x)*diff(p,y))*I1+(-1/q*diff(q,y)*p^2 -2*diff(p,y)*p)*diff(I1,x)+(-diff(p,x)*p-2/q*diff(q,x)*p^2) *diff(I1,y)+p^3*diff(q,x)*diff(q,`$`(y,2))+I5 +2*q*diff(p,x)*diff(p,y)*diff(q,y)*p-2*p^2*diff(p,y) *diff(q,y)*diff(q,x)-q*p*diff(I3,x)-q*diff(q,x)*p^2 *diff(p,`$`(y,2))+p^2*diff(I1,x,y)-q*p^2*diff(q,`$`(y,2))*diff(p,x) +2*q*diff(q,x)*diff(p,y)^2*p+q^2*diff(p,x)*p *diff(p,`$`(y,2))-1/q*p^3*diff(I4,y)-2*q^2*diff(p,x)*diff(p,y)^2 +I2*(-2*p^3*q^2*diff(q,x)-p*q^2*I1+4*diff(p,y)*q^4*p -2*q^3*p^2*diff(p,x))/q^4/p +diff(I2,x)/q*p^2-p*diff(I2,y) ); fi; fi; # # operators of the symbol pXXY # if a30=0 and a12=0 and a03=0 and a21<>0 then p=a21; if n=1 then RETURN(-I1); elif n=2 then RETURN(-I2); elif n=3 then RETURN(I3-I2); elif n=4 then RETURN(I4-2*diff(I1,y)); elif n=5 then RETURN(diff(I3,x)+diff(I4,y)-diff(I1,`$`(y,2))-I5-diff(I2,x)/2); fi; fi; # # operators of the symbol pXXX # if a03=0 and a12=0 and a21=0 and a30<>0 then p=a30; # # consider the case a02=0 and a11<>0, which is denoted by introducing # artificial invariant i5=0 # if I5=0 then if n=1 then RETURN(-I1); elif n=2 then RETURN(I2-diff(I1,x)); elif n=3 then RETURN(I3-(diff(I1,x,y)*I1-diff(I1,x)*diff(I1,y))/I1); elif n=4 then RETURN(I3-I4+diff(I2,y)-diff(I1,x,y) +1/I1*diff(I1,x)*diff(I1,y) -1/I1*I2*diff(I1,y)); else ERROR(`LPDO__Inv_Transposed: there is only four invariants in this case, therefore, n=1,2,3,4 are the only allowed values`); fi; fi; fi; ERROR(`LPDO__Inv_Transposed: the procedures does not deal yet with such operators`) end: # #------------------------------------------------ # compute common obstacle to factorization of the type (S1)(S2)(S3), # Si= ri X + si Y for the operator a with the symbol S1*S2*S3 # outputs F1,F2,F3,obst, s.t. a=F1 F2 F3 + obst # # LPDO__obstacle:= proc(r1,s1,r2,s2,r3,s3,a::array,d) local F1,F2,F3,k1,k2,k3,Appr,obst,sol; if LPDO__nvars() <> 2 then ERROR(` LPDO__obstacle: illegal number of vars`); fi; if LPDO__degree(a) <> 3 then ERROR(` LPDO__obstacle: illegal degree of operator`); fi; F1:=LPDO__create1(r1,s1,k1): F2:=LPDO__create1(r2,s2,k2): F3:=LPDO__create1(r3,s3,k3): # LPDO__print(a); Appr:=LPDO__mult(F1,LPDO__mult(F2,F3)): # LPDO__print(Appr); obst:=LPDO__minus(a,Appr); if LPDO__degree(obst) > 2 then ERROR(` LPDO__obstacle: symbol operator <> S1*S2*S3`); fi; sol:=solve({obst[3],obst[4],obst[5]},{k1,k2,k3}); assign(sol); obst:=LPDO__simplify(obst): F1:=LPDO__simplify(F1): F2:=LPDO__simplify(F2): F3:=LPDO__simplify(F3): if d=1 then print(`The factors of the incomplete factorization in order. First factor:`); LPDO__print(F1); print(`Second factor:`); LPDO__print(F2); print(`Third factor:`); LPDO__print(F3); print(`The obstacle:`); LPDO__print(obst); fi; RETURN([eval(F1),eval(F2),eval(F3),eval(obst)]); end: #--------------------------------------------------- # One step of one of the two Laplace transformations # checked # returns L1, if i=1, and L_{-1} if i=-1. # LPDO__Laplace_trans_classic:= proc(i::integer,L::array) local a,b,c,h,k,L1,a1,b1,c1,x,y; if LPDO__nvars() <> 2 then ERROR(` LPDO__Laplace_trans_classic: illegal number of vars`); fi; if LPDO__degree(L) <> 2 then ERROR(` LPDO__Laplace_trans_classic: illegal degree of operator`); fi; if (i<>1 and i<>-1) then ERROR(`the first argument should be 1 or -1`); fi; h:=LPDO__Lapl_h(L); k:=LPDO__Lapl_k(L); a:=LPDO__value(L,[1,0]); b:=LPDO__value(L,[0,1]); c:=LPDO__value(L,[0,0]); x:=LPDO__vars[1]: y:=LPDO__vars[2]: if i=1 then if h=0 then ERROR(`This Laplace transformation cannot be performed as Laplace invariant h is 0`); else if LPDO__warnings=true then print(`This Laplace transformation can be performed iff Laplace invariant h is not zero. In your case it is`, simplify(eval(LPDO__Lapl_h(L)))): fi: a1:= a-diff(h,y)/h; b1:= b; c1:= a1*b1+diff(b,y) - h; L1:= LPDO__create2( 0, 1, 0, a1, b1, c1): RETURN(LPDO__simplify(eval(L1))); fi; fi; if i=-1 then if k=0 then ERROR(`This Laplace transformation cannot be performed as Laplace invariant k is 0`); else if LPDO__warnings=true then print(`This Laplace transformation can be performed iff Laplace invariant k is not zero. In your case it is`, simplify(eval(LPDO__Lapl_k(L)))): fi: a1:= a; b1:= b-diff(k,x)/k; c1:= a1*b1 + diff(a,x) - k: L1:= LPDO__create2( 0, 1, 0, a1, b1, c1): RETURN(LPDO__simplify(eval(L1))); fi; fi; end: #---------------------------------------------------------------------------------------- # # Laplace transformations for operators of the form S=Dxx + Dyy + aDx + bDy + c # i=1 stands for the first transforamtion # i=-1 stands for the second transforamtion # LPDO__Laplace_trans_gen_form:= proc(i::integer,S::array) local a,b,c,S1,a1,b1,c1; if LPDO__nvars() <> 2 then ERROR(` LPDO__Laplace_trans_gen_form: illegal number of vars`); fi; if LPDO__degree(S) <> 2 then ERROR(` LPDO__Laplace_trans_gen_form: illegal degree of operator`); fi; if (i<>1 and i<>-1) then ERROR(`the first argument should be 1 or -1`); fi; a:=LPDO__value(S,[1,0]); b:=LPDO__value(S,[0,1]); c:=LPDO__value(S,[0,0]); if i=1 then if LPDO__Lapl_h(S)=0 then ERROR(` LPDO__Laplace_trans(1): h=0, the transformation is not defined for this case`); fi; a1:= -(2*b*diff(b,LPDO__vars[1])+a*b^2-4*a*c-2*I*diff(b,`$`(LPDO__vars[1],2))+2*I*b*diff(b,LPDO__vars[2]) +a^3-4*I*diff(c,LPDO__vars[2])-4*diff(c,LPDO__vars[1])-2*diff(a,`$`(LPDO__vars[1],2)) -2*I*a*diff(b,LPDO__vars[1])-2*diff(a,`$`(LPDO__vars[2],2))-2*I*diff(b,`$`(LPDO__vars[2],2)) -2*a*diff(b,LPDO__vars[2])+4*I*a*diff(a,LPDO__vars[2]))/(-4*c+b^2+a^2-2*I*diff(b,LPDO__vars[1]) -2*diff(b,LPDO__vars[2])-2*diff(a,LPDO__vars[1])+2*I*diff(a,LPDO__vars[2])); b1:= -(b*a^2+b^3-4*I*b*diff(b,LPDO__vars[1])+4*I*diff(c,LPDO__vars[1])+2*I*diff(a,`$`(LPDO__vars[1],2)) +2*I*b*diff(a,LPDO__vars[2])-4*diff(c,LPDO__vars[2]) -2*diff(b,`$`(LPDO__vars[1],2))+2*I*diff(a,`$`(LPDO__vars[2],2))-4*b*c-2*I*a*diff(a,LPDO__vars[1]) -2*diff(b,`$`(LPDO__vars[2],2))+2*a*diff(a,LPDO__vars[2])-2*b*diff(a,LPDO__vars[1]))/(-4*c+b^2+a^2 -2*I*diff(b,LPDO__vars[1])-2*diff(b,LPDO__vars[2])-2*diff(a,LPDO__vars[1])+2*I*diff(a,LPDO__vars[2])); c1:= (-4*c^2+a*b*diff(b,LPDO__vars[2])*I+a*b*diff(b,LPDO__vars[1])+b*a*diff(a,LPDO__vars[2]) +a^2*diff(b,LPDO__vars[1])*I+b*diff(a,`$`(LPDO__vars[1],2))*I-I*a*diff(b,`$`(LPDO__vars[1],2)) -I*a*diff(b,`$`(LPDO__vars[2],2))-a*diff(a,`$`(LPDO__vars[1],2))+2*diff(b,LPDO__vars[1])^2 -b*diff(b,`$`(LPDO__vars[2],2))-b*diff(b,`$`(LPDO__vars[1],2))+b^2*diff(b,LPDO__vars[2]) -2*c*diff(b,LPDO__vars[2])+a^2*diff(a,LPDO__vars[1])+c*b^2+c*a^2-2*c*diff(a,LPDO__vars[1]) -2*a*diff(c,LPDO__vars[1])-2*b*diff(c,LPDO__vars[2])-a*diff(a,`$`(LPDO__vars[2],2)) -4*diff(b,LPDO__vars[1])*diff(a,LPDO__vars[2])-I*b*a*diff(a,LPDO__vars[1]) +2*diff(a,LPDO__vars[2])^2-2*I*a*diff(c,LPDO__vars[2])-I*b^2*diff(a,LPDO__vars[2]) -2*I*diff(b,LPDO__vars[1])*diff(a,LPDO__vars[1])+2*I*b*diff(c,LPDO__vars[1]) +2*I*diff(a,LPDO__vars[1])*diff(a,LPDO__vars[2])-2*I*diff(b,LPDO__vars[1])*diff(b,LPDO__vars[2]) +2*I*diff(b,LPDO__vars[2])*diff(a,LPDO__vars[2])-6*I*c*diff(b,LPDO__vars[1])+6*I*c*diff(a,LPDO__vars[2]) +b*diff(a,`$`(LPDO__vars[2],2))*I)/(-4*c+b^2+a^2-2*I*diff(b,LPDO__vars[1])-2*diff(b,LPDO__vars[2]) -2*diff(a,LPDO__vars[1])+2*I*diff(a,LPDO__vars[2])); S1:= LPDO__create2( 1, 0, 1, a1, b1, c1): RETURN(LPDO__simplify(S1)); fi; if i=-1 then if LPDO__Lapl_k(S)=0 then ERROR(` LPDO__Laplace_trans(-1): k=0, the transformation is not defined for this case`); fi; a1:= (-2*I*diff(b,`$`(LPDO__vars[1],2))+4*a*c+2*a*diff(b,LPDO__vars[2])+2*diff(a,`$`(LPDO__vars[2],2)) -2*I*diff(b,`$`(LPDO__vars[2],2))-a^3+4*diff(c,LPDO__vars[1])+2*diff(a,`$`(LPDO__vars[1],2)) -4*I*diff(c,LPDO__vars[2])-a*b^2-2*b*diff(b,LPDO__vars[1])-2*I*a*diff(b,LPDO__vars[1]) +4*I*a*diff(a,LPDO__vars[2])+2*I*b*diff(b,LPDO__vars[2]))/(-4*c+b^2+a^2 -2*diff(a,LPDO__vars[1])-2*I*diff(a,LPDO__vars[2])+2*I*diff(b,LPDO__vars[1])-2*diff(b,LPDO__vars[2])); b1:= -(-2*b*diff(a,LPDO__vars[1])-4*diff(c,LPDO__vars[2])+2*a*diff(a,LPDO__vars[2])+b^3-4*b*c-2*diff(b,`$`(LPDO__vars[1],2))-2*I*b*diff(a,LPDO__vars[2])-2*I*diff(a,`$`(LPDO__vars[2],2)) +2*I*a*diff(a,LPDO__vars[1])+4*I*b*diff(b,LPDO__vars[1])+b*a^2-4*I*diff(c,LPDO__vars[1]) -2*I*diff(a,`$`(LPDO__vars[1],2))-2*diff(b,`$`(LPDO__vars[2],2))) /(-4*c+b^2+a^2-2*diff(a,LPDO__vars[1])-2*I*diff(a,LPDO__vars[2])+2*I*diff(b,LPDO__vars[1]) -2*diff(b,LPDO__vars[2])); c1:= -(b*diff(a,`$`(LPDO__vars[1],2))*I-I*a*diff(b,`$`(LPDO__vars[1],2))-b^2*diff(b,LPDO__vars[2]) -2*diff(a,LPDO__vars[2])^2-2*diff(b,LPDO__vars[1])^2+a*diff(a,`$`(LPDO__vars[2],2)) +4*c^2+b*diff(a,`$`(LPDO__vars[2],2))*I-I*a*diff(b,`$`(LPDO__vars[2],2))+b*diff(b,`$`(LPDO__vars[1],2)) +a*diff(a,`$`(LPDO__vars[1],2))+b*diff(b,`$`(LPDO__vars[2],2))-a*b*diff(b,LPDO__vars[1]) +a^2*diff(b,LPDO__vars[1])*I-b*a*diff(a,LPDO__vars[2])-2*I*diff(b,LPDO__vars[1])*diff(b,LPDO__vars[2]) -I*b^2*diff(a,LPDO__vars[2])+2*I*diff(a,LPDO__vars[1])*diff(a,LPDO__vars[2])-2*I*a*diff(c,LPDO__vars[2]) -6*I*c*diff(b,LPDO__vars[1])+2*I*b*diff(c,LPDO__vars[1])-2*I*diff(a,LPDO__vars[1])*diff(b,LPDO__vars[1]) +6*I*c*diff(a,LPDO__vars[2])+2*I*diff(a,LPDO__vars[2])*diff(b,LPDO__vars[2])+a*b*diff(b,LPDO__vars[2])*I -I*b*a*diff(a,LPDO__vars[1])+2*b*diff(c,LPDO__vars[2])+4*diff(a,LPDO__vars[2])*diff(b,LPDO__vars[1]) +2*c*diff(b,LPDO__vars[2])+2*c*diff(a,LPDO__vars[1])-a^2*diff(a,LPDO__vars[1])+2*a*diff(c,LPDO__vars[1]) -c*b^2-c*a^2)/(-4*c+b^2+a^2-2*diff(a,LPDO__vars[1])-2*I*diff(a,LPDO__vars[2])+2*I*diff(b,LPDO__vars[1]) -2*diff(b,LPDO__vars[2])); S1:= LPDO__create2( 1, 0, 1, a1, b1, c1): RETURN(LPDO__simplify(S1)); fi; end: #---------------------------------------------------------------------------------------- # # for L=Dxy+... or L=Dxx + Dyy returns Ln for n=+/-1 # LPDO__Lapl_trans := proc(n::integer, L::array) if LPDO__nvars() <> 2 then ERROR(`LPDO__Lapl_trans : illegal number of vars`); fi; if LPDO__degree(L) <> 2 then ERROR(`LPDO__Lapl_trans : illegal degree of operator`); fi; if LPDO__value(L, [2,0])=0 and LPDO__value(L, [1,1])=1 and LPDO__value(L, [0,2])=0 then RETURN(LPDO__Laplace_trans_classic(n,L)); elif LPDO__value(L, [2,0])=1 and LPDO__value(L, [1,1])=0 and LPDO__value(L, [0,2])=1 then RETURN(LPDO__Laplace_trans_gen_form(n,L)); else ERROR(`LPDO__Lapl_trans: operator should have the symbol of the form XY or X^2 + Y^2`); fi; end: # # #---------------------------------------------------------------------------------------- # # (k,h) -> (k_1,h1) for L=Dxy+... # checked # LPDO__Lapl_chain_step1 := proc(k,h) local h1,k_1: if h=0 and k=0 then RETURN([0,0]); elif k=0 then RETURN( [0 , simplify(eval(2*h-k - ((diff(h, y, x))/h-(diff(h, x))*(diff(h, y))/h^2) ))] ); elif h=0 then RETURN( [simplify(eval(2*k-h - ((diff(k, y, x))/k-(diff(k, x))*(diff(k, y))/k^2) )) ,0] ); else h1:=eval(2*h-k - (diff(h,x,y)/h - diff(h,x)*diff(h,y)/h^2)): k_1:=eval(2*k-h - (diff(k,x,y)/k - diff(k,x)*diff(k,y)/k^2)): RETURN([k_1 ,h1]); fi; end: # #---------------------------------------------------------------------------------------- # checked # h_2 -> h_1 -> h -> h1 # || || || || # k_2 ->k_1 -> k -> k1 ->k2 # # L -> [[k,h],[k_1,h1],...,[k_n,hn]] # n steps # LPDO__Lapl_chain_classic := proc(n::integer, L::array) local a,b,c,h,k,i,res,h_new,k_new,hp,kp,x,y; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; if LPDO__degree(L) <> 2 then ERROR(`illegal degree of operator`); fi; if LPDO__value(L, [2,0]) <> 0 or LPDO__value(L, [1,1]) <> 1 or LPDO__value(L, [0,2]) <> 0 then ERROR(`Illegal symbol of operator`); fi; if n<0 then ERROR(`The first argument denotes the number of interations of Laplace transformation method and must be non-negative`); fi; h := simplify(eval(LPDO__Lapl_h(L))); k := simplify(eval(LPDO__Lapl_k(L))); x:=LPDO__vars[1]: y:=LPDO__vars[2]: if n=0 then RETURN([[k,h]]): fi: res:=[[k,h],map(el->simplify(eval(el)),LPDO__Lapl_chain_step1(k,h))]; if n>1 then for i from 3 to n+1 do # i is the number of the new sub-list in list res if res[i-1][2]=0 then h_new:=0: else hp:=res[i-1][2]: h_new:=2*hp-res[i-2][2]-(diff(hp,x,y)/hp - diff(hp,x)*diff(hp,y)/hp^2): fi: if res[i-1][1]=0 then k_new:=0: else kp:=res[i-1][1]: k_new:=2*kp-res[i-2][1]-(diff(kp,x,y)/kp - diff(kp,x)*diff(kp,y)/kp^2): fi: res:=[op(res),[simplify(eval(k_new)),simplify(eval(h_new))]]; od; fi: RETURN(res); end: # #---------------------------------------------------------------------------------------- # # finds the chain of invariants [[k,h],[k_1,h1],...,[k_n,hn]] for operators of the form # S=Dxx + Dyy + aDx + bDy + c # LPDO__Lapl_chain_gen_form := proc(n::integer, S::array) local i,S1,S_1,res; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; if LPDO__degree(S) <> 2 then ERROR(`illegal degree of operator`); fi; if LPDO__value(S, [2,0]) <> 1 then ERROR(`illegal symbol[2,0] of operator`); fi; if LPDO__value(S, [1,1]) <> 0 then ERROR(`illegal symbol[1,1] of operator`); fi; if LPDO__value(S, [0,2]) <> 1 then ERROR(`illegal symbol[0,2] of operator`); fi; res:=[[LPDO__Lapl_k(S),LPDO__Lapl_h(S)]]; S1:=S; S_1:=S; for i from 1 to n-1 do S1:=LPDO__Laplace_trans_gen_form(1,S1); S_1:=LPDO__Laplace_trans_gen_form(-1,S_1); res:=[op(res),[LPDO__Lapl_k(S_1),LPDO__Lapl_h(S1)]]; od; RETURN(res); end: # #---------------------------------------------------------------------------------------- # # h_2 -> h_1 -> h -> h1 # || || || || # k_2 ->k_1 -> k -> k1 ->k2 # # L -> [[k,h],[k_1,h1],...,[k_n,hn]] # n steps # # L=Dxy+... or L=Dxx + Dyy # LPDO__Lapl_chain := proc(n::integer, L::array) if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; if LPDO__degree(L) <> 2 then ERROR(`illegal degree of operator`); fi; if LPDO__value(L, [2,0])=0 and LPDO__value(L, [1,1])=1 and LPDO__value(L, [0,2])=0 then RETURN(LPDO__Lapl_chain_classic(n,L)); elif LPDO__value(L, [2,0])=1 and LPDO__value(L, [1,1])=0 and LPDO__value(L, [0,2])=1 then RETURN(LPDO__Lapl_chain_gen_form(n,L)); else ERROR(`operator should have the symbol of the form XY or X^2 + Y^2`); fi; end: #------------------------------------------------------------------------------------------------ # For the chain of Laplace invaraints # returns h_{i+1} for given hp=h_i and hpp=h_{i-1} # LPDO__Lapl_inv_h_next:= proc(hpp,hp) local x,y: if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; x:=LPDO__vars[1]: y:=LPDO__vars[2]: RETURN(simplify(eval(2*hp-hpp-(diff(hp,x,y)/hp - diff(hp,x)*diff(hp,y)/hp^2)))): end: #------------------------------------------------------------------------------------------------ # For the chain of Laplace invaraints # returns k_{-i-1} for given kp=k_{-i} and kpp=k_{-i+1} # LPDO__Lapl_inv_k_next:= proc(kpp,kp) local x,y: if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; x:=LPDO__vars[1]: y:=LPDO__vars[2]: RETURN(simplify(eval(2*kp-kpp-(diff(kp,x,y)/kp - diff(kp,x)*diff(kp,y)/kp^2)))): end: # #---------------------------------------------------------------------------------------- # # solver of PDEs u_xy+au_x + bu_y +cu=0 corresponding to factorizable L # Factorization_solver:=proc(L::array) local a,b,c,h,k,x,y: if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; if LPDO__degree(L) <> 2 then ERROR(`Laplace transformation method works for hyperbolic operators of second order only`); fi; h:=LPDO__Lapl_h(L): k:=LPDO__Lapl_k(L): a:=LPDO__value(L,[1,0]): b:=LPDO__value(L,[0,1]): c:=LPDO__value(L,[0,0]): x:=LPDO__vars[1]: y:=LPDO__vars[2]: if h=0 then RETURN(simplify(eval(exp(-int(a,y))*(_F1(x) + int(_F2(y)*exp(int(a,y)-int(b,x)),y))))): fi: if k=0 then RETURN(simplify(eval(exp(-int(b,x))*(_F1(y) + int(_F2(x)*exp(int(b,x)-int(a,y)),x))))): fi: end: #-------------------------------------------------------------------------------------- # n -number of interations # Solve by Laplace method by < or = n interations. # Laplace_solver_op:=proc(n::integer, L::array) local chain, i,j,a,b,c,x,y,u,ui,Li,L_i,kp,kpp,hp,hpp; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; if LPDO__degree(L) <> 2 then ERROR(`Laplace transformation method works for hyperbolic operators of second order only`); fi; chain:=[[LPDO__Lapl_k(L),LPDO__Lapl_h(L)]]: #print(`chain_is`,chain): a:=LPDO__value(L,[1,0]): b:=LPDO__value(L,[0,1]): c:=LPDO__value(L,[0,0]): x:=LPDO__vars[1]: y:=LPDO__vars[2]: Li:=L: L_i:=L: hpp:=chain[1][1]: #=k hp:=chain[1][2]: #=h kpp:=chain[1][2]: #=h kp:=chain[1][1]: #=k #LPDO__print(Li): LPDO__print(L_i): if chain[1][2]=0 or chain[1][1]=0 then RETURN(simplify(eval(Factorization_solver(L)))): fi: for i from 1 to n do # i - current iteration chain:=[op(chain),[LPDO__Lapl_inv_k_next(kpp,kp),LPDO__Lapl_inv_h_next(hpp,hp)]]: #print(`iteration,current_chain`,i,chain); if chain[i+1][2]=0 then # hi=0 #print(`h is zero`): Li:=LPDO__Lapl_trans(1,Li): #print(`current Li_ is`): #LPDO__print(Li): ui:=Factorization_solver(Li): u:=ui*exp(int(b,x)): #print(`take in brackets u=`,u): for j from i-1 by -1 to 0 do #print(`run jfrom`,i-1,`to_ 0`): #print(`current j=`,j): u:=1/chain[j+1][2]*diff(u,x): od: RETURN(simplify(eval(exp(-int(b,x))*u))): fi: if chain[i+1][1]=0 then # ki=0 print(`k is zero`): L_i:=LPDO__Lapl_trans(-1,L_i): ui:=Factorization_solver(L_i): u:=ui*exp(int(a,y)): for j from i-1 by -1 to 0 do u:=1/chain[j+1][1]*diff(u,y): od: RETURN(simplify(eval(exp(-int(a,y))*u))): fi: Li:=LPDO__Lapl_trans(1,Li): L_i:=LPDO__Lapl_trans(-1,L_i): hpp:=chain[i][2]: hp:=chain[i+1][2]: kpp:=chain[i][1]: kp:=chain[i+1][1]: od: RETURN(NULL): end: #-------------------------------------------------------------------------------------- # n -number of interations # Solve by Laplace method by < or = n interations. # L is either operator in the form D_xy + a D_x + b D_y + c or # L is the corresponding PDE u_xy + a u_x + b u_y + c u = 0 # one may need to spesify the third argument func, the dependable variable with respect to which is our PDE Laplace_solver:=proc(n::integer, L, func:=0) local els,i,u,a,b,c,opp; if type(L,array)=true then RETURN(Laplace_solver_op(n,L)): fi: if func<>0 then u:=func: else els:=indets(L); u:=0; for i from 1 to nops(els) while u=0 do if evalb(diff(els[i],x,y) in els)=true then u:=els[i]: fi: od: fi: c:=PDEtools[dsubs](eval(u)=1,L); a:=PDEtools[dsubs](eval(u)=x,L)-c*x; b:=PDEtools[dsubs](eval(u)=y,L)-c*y; opp:= Create_LPDO_from_coeff([c,a,b,0,1]): if simplify(eval(LPDO__apply(opp,u)-L))=0 then RETURN(Laplace_solver_op(n,opp)): else if func=0 then ERROR(`Either the form of the equation is not u_xy + a u_x + b u_y + c u = 0, or we cannot guess the dependent variable, use the third argument to specify it.`); else ERROR(`It seems that the form of the equation is not u_xy + a u_x + b u_y + c u = 0`); fi: fi: end: #---------------------------------------------------------------------------------------- # LPDO__factor111:= proc(r1,s1,r2,s2,r3,s3,a::array) local F1,F2,F3,k1,k2,k3,Appr,obst,sol,test; F1:=LPDO__create1(r1,s1,k1): F2:=LPDO__create1(r2,s2,k2): F3:=LPDO__create1(r3,s3,k3): Appr:=LPDO__mult(F1,LPDO__mult(F2,F3)): obst:=LPDO__minus(a,Appr); sol:=pdsolve({obst[3],obst[4],obst[5],obst[2],obst[1],obst[0]},{k1,k2,k3}); assign(sol); obst:=LPDO__simplify(obst): F1:=LPDO__simplify(F1): F2:=LPDO__simplify(F2): F3:=LPDO__simplify(F3): print(`The factors of the incomplete factorization in order. First factor:`); LPDO__print(F1); print(`Second factor:`); LPDO__print(F2); print(`Third factor:`); LPDO__print(F3); print(`The obstacle:`); LPDO__print(obst); RETURN([eval(F1),eval(F2),eval(F3),eval(obst)]); end: # #---------------------------------------------------------------------------------------- # # does for the operator "a" with the symbol XY(pX+qY) # output conditions of factorizations of type (pX+qY)(XY) in terms of the invariants of THE full system of invariants # LPDO__Fact_Cond_S_XY:= proc(p,q,I1,I2,I3,I4,I5) local Cond10, Cond00,Ir,Ip,Iq,I1x,I1y,I2x,I2y,px,py,qx,qy,co1,z00_; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; Cond10:=I3*q^3-diff(I1,y)*q^2*p+diff(q,y)*I1*q*p+2*diff(p,y)*I1*q^2 -I4*p^3+q*p^2*diff(I1,x)-2*diff(q,x)*p^2*I1-p*q*diff(p,x)*I1 -3*p*q*I2; Ir:= -p^2*q^2*diff(I4,y)+1/2*p*q^3*diff(I1,x,y)-p^3*q*diff(I4,x)-3/2*p*q^2*diff(q,x)*diff(I1,y)+1/p*q^3*I5 +p^2*q^2*diff(I1,`$`(x,2))+(q^2*diff(p,x)*diff(q,y)-p*q^2*diff(p,`$`(x,2))-1/p*q^3*diff(p,x)*diff(p,y) -3/2*p*q^2*diff(q,x,y)+2*q^2*diff(p,x)^2+4*p*q*diff(q,x)*diff(p,x)-2*p^2*q*diff(q,`$`(x,2)) +q^2*diff(q,x)*diff(p,y)+5*p*q*diff(q,x)*diff(q,y)+6*p^2*diff(q,x)^2)*I1+(3*p^3*diff(q,x)+3*p^2*q*diff(q,y))*I4 -q*I1*diff(I1,x)+p*I1*I4+(2*diff(q,x)+1/p*q*diff(p,x))*I1^2+(-4*p^2*q*diff(q,x)-2*p*q^2*diff(p,x) -3/2*p*q^2*diff(q,y))*diff(I1,x); Ip:= p: Iq := q: I1x := diff(I1,x): I1y := diff(I1,y): I2x := diff(I2,x): I2y := diff(I2,y): px:= diff(Ip,x): py := diff(Ip,y): qx := diff(Iq,x): qy := diff(Iq,y): co1 := + 4*Iq^2*px + Iq*I1/Ip +4*Ip*Iq*qx + 2*Iq^2*qy + 2*Iq^3*py/Ip; z00_ := -2*Ip*Iq^2*I2x - Iq^3*I2y + co1*I2; Cond00:= Ir + z00_; RETURN([simplify(eval(Cond10)),simplify(eval(Cond00))]); ERROR(`the procedure is implemented for first arguments p,q only`); end: #---------------------------------------------------------------------------------------- # # does for the operator "a" with the symbol XYS # output conditions of factorizations of type (S)(X)(Y) in terms of the invariants of THE full system of invariants # LPDO__Fact_Cond_S_X_Y:= proc(p,q,i1,i2,i3,i4,i5) local Cond, list1; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; Cond:= -i4*p^2+q*p*diff(i1,x)-2*diff(q,x)*p*i1-q*diff(p,x)*i1-q*i2; list1:=LPDO__Fact_Cond_S_XY(p,q,i1,i2,i3,i4,i5); RETURN([op(list1),simplify(eval(Cond))]); ERROR(`ERROR in LPDO__Fact_Cond_S_X_Y`); end: #---------------------------------------------------------------------------------------- # # does for the operator "a" with the symbol XYS # output conditions of factorizations of type (Y)(X)(S) in terms of the invariants of THE full system of invariants # LPDO__Fact_Cond_Y_X_S:= proc(p,q,i1,i2,i3,i4,i5) local i1t, i2t, i3t, i4t, i5t; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; i1t:=LPDO__Inv_Transposed(1,i1,i2,i3,i4,i5,0,p,q,0); i2t:=LPDO__Inv_Transposed(2,i1,i2,i3,i4,i5,0,p,q,0); i3t:=LPDO__Inv_Transposed(3,i1,i2,i3,i4,i5,0,p,q,0); i4t:=LPDO__Inv_Transposed(4,i1,i2,i3,i4,i5,0,p,q,0); i5t:=LPDO__Inv_Transposed(5,i1,i2,i3,i4,i5,0,p,q,0); RETURN(LPDO__Fact_Cond_S_X_Y(-p,-q,i1t,i2t,i3t,i4t,i5t)); ERROR(`ERROR in LPDO__Fact_Cond_Y_X_S`); end: #---------------------------------------------------------------------------------------- # # does for the operator "a" with the symbol XYS # output conditions of factorizations of type (S)(Y)(X) in terms of the invariants of THE full system of invariants # LPDO__Fact_Cond_S_Y_X:= proc(p,q,i1,i2,i3,i4,i5) local Cond, list1; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; Cond:=-i4*p^2+q*p*diff(i1,x)-2*diff(q,x)*p*i1-diff(p,x)*q*i1-2*q*i2; list1:=LPDO__Fact_Cond_S_XY(p,q,i1,i2,i3,i4,i5); RETURN([op(list1),simplify(eval(Cond))]); ERROR(`ERROR in LPDO__Fact_Cond_S_X_Y`); end: #---------------------------------------------------------------------------------------- # # does for the operator "a" with the symbol XYS # output conditions of factorizations of type (X)(Y)(S) in terms of the invariants of THE full system of invariants # LPDO__Fact_Cond_X_Y_S:= proc(p,q,i1,i2,i3,i4,i5) local i1t, i2t, i3t, i4t, i5t; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; i1t:=LPDO__Inv_Transposed(1,i1,i2,i3,i4,i5,0,p,q,0); i2t:=LPDO__Inv_Transposed(2,i1,i2,i3,i4,i5,0,p,q,0); i3t:=LPDO__Inv_Transposed(3,i1,i2,i3,i4,i5,0,p,q,0); i4t:=LPDO__Inv_Transposed(4,i1,i2,i3,i4,i5,0,p,q,0); i5t:=LPDO__Inv_Transposed(5,i1,i2,i3,i4,i5,0,p,q,0); RETURN(LPDO__Fact_Cond_S_Y_X(-p,-q,i1t,i2t,i3t,i4t,i5t)); ERROR(`ERROR in LPDO__Fact_Cond_X_Y_S`); end: #---------------------------------------------------------------------------------------- # # does for the operator "a" with the symbol XY(pX+qY) # output conditions of factorizations of type (X)(SY) in terms of the invariants of THE full system of invariants # LPDO__Fact_Cond_X_SY:= proc(p,q,i1,i2,i3,i4,i5) local Cond10, Cond00,Ir; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; Cond10:=-2/q*diff(q,x)*diff(p,x)-1/q*p*diff(q,`$`(x,2))+diff(p,`$`(x,2))+1/q^2*i4+2/q^2*p*diff(q,x)^2; Ir:= -3/2*diff(q,x)*q*p^2*diff(i1,y)-q^3*p*diff(i3,x)+i5*q^2+1/2*q^2*p^2*diff(i1,x,y)+(-1/2*diff(q,y)*p^2*q-diff(p,y)*p*q^2)*diff(i1,x)+(diff(q,x)*p*q^2+2*diff(p,x)*q^3)*i3+(diff(p,x)*diff(p,y)*q^2+2*diff(q,x)*diff(q,y)*p^2-1/2*diff(q,x,y)*p^2*q+3*diff(q,x)*diff(p,y)*p*q-diff(p,x,y)*p*q^2)*i1 ; Cond00:= Ir+ (-4*p^2*diff(q,x))*i2 +p^2*q*diff(i2,x); RETURN([simplify(eval(Cond10)),simplify(eval(Cond00))]); ERROR(`the procedure is implemented for first arguments p,q only`); end: # #---------------------------------------------------------------------------------------- # # does for the operator "a" with the symbol XYS # output conditions of factorizations of type (S)(Y)(X) in terms of the invariants of THE full system of invariants # LPDO__Fact_Cond_X_S_Y:= proc(p,q,i1,i2,i3,i4,i5) local Cond, list1; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; Cond:=1/p^2*i3-1/q/p*diff(i1,y)+1/q^2/p*diff(q,y)*i1+2/q/p^2*diff(p,y)*i1-2*i2/p/q^2; list1:=LPDO__Fact_Cond_X_SY(p,q,i1,i2,i3,i4,i5); RETURN([op(list1),simplify(eval(Cond))]); ERROR(`ERROR in LPDO__Fact_Cond_S_X_Y`); end: #---------------------------------------------------------------------------------------- # # does for the operator "a" with the symbol XYS # output conditions of factorizations of type (Y)(S)(X) in terms of the invariants of THE full system of invariants # LPDO__Fact_Cond_Y_S_X:= proc(p,q,i1,i2,i3,i4,i5) local i1t, i2t, i3t, i4t, i5t; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; i1t:=LPDO__Inv_Transposed(1,i1,i2,i3,i4,i5,0,p,q,0); i2t:=LPDO__Inv_Transposed(2,i1,i2,i3,i4,i5,0,p,q,0); i3t:=LPDO__Inv_Transposed(3,i1,i2,i3,i4,i5,0,p,q,0); i4t:=LPDO__Inv_Transposed(4,i1,i2,i3,i4,i5,0,p,q,0); i5t:=LPDO__Inv_Transposed(5,i1,i2,i3,i4,i5,0,p,q,0); RETURN(LPDO__Fact_Cond_X_S_Y(-p,-q,i1t,i2t,i3t,i4t,i5t)); ERROR(`ERROR in LPDO__Fact_Cond_Y_S_X`); end: #---------------------------------------------------------------------------------------- # # does for the operator "a" with the symbol XY(pX+qY) # output conditions of factorizations of type (Y)(SX) in terms of the invariants of THE full system of invariants # LPDO__Fact_Cond_Y_SX:= proc(p,q,i1,i2,i3,i4,i5) local Cond01, Cond00,Ir; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; Cond01:=1/p^2*i3-2/p*diff(p,y)*diff(q,y)-1/p*q*diff(p,`$`(y,2))+diff(q,`$`(y,2))+2/p^2*diff(p,y)^2*q; Ir:= -p^3*q*diff(i4,y)-3/2*diff(q,x)*diff(i1,y)*p^2*q+i5*q^2+1/2*diff(i1,x,y)*p^2*q^2+(p^2*q*diff(p,y)+2*diff(q,y)*p^3)*i4+(-diff(p,y)*p*q^2-1/2*diff(q,y)*p^2*q)*diff(i1,x)+(3*diff(q,x)*diff(p,y)*p*q+3*diff(q,x)*diff(q,y)*p^2-3/2*diff(q,x,y)*p^2*q)*i1; Cond00:= Ir+ (4*q^2*diff(p,y))*i2 -q^2*p*diff(i2,y); RETURN([simplify(eval(Cond01)),simplify(eval(Cond00))]); ERROR(`ERROR in LPDO__Fact_Cond_Y_SX`); end: # #---------------------------------------------------------------------------------------- # # does for the operator "a" with the symbol XY(pX+qY) # output conditions of factorizations of type (SY)(X) in terms of the invariants of THE full system of invariants # LPDO__Fact_Cond_SY_X:= proc(p,q,i1,i2,i3,i4,i5) local Cond10, Cond00; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; Cond10:=-(2*diff(q,x)*i1*p-diff(i1,x)*q*p+diff(p,x)*i1*q+i4*p^2+2*q*i2)/q^2/p^2; Cond00:= p^2*i4*diff(p,y)*q-p^3*diff(i4,y)*q+1/2*p^2*diff(i1,x,y)*q^2+i5*q^2-i2*i1 -p*diff(i2,y)*q^2+2*p^3*i4*diff(q,y)+4*i2*q^2*diff(p,y)+3*i1*diff(q,x)*diff(q,y)*p^2 -3/2*i1*diff(q,x,y)*p^2*q+3*i1*diff(q,x)*diff(p,y)*p*q-1/2*p^2*diff(i1,x)*q*diff(q,y) -p*diff(i1,x)*q^2*diff(p,y)-3/2*q*p^2*diff(i1,y)*diff(q,x); RETURN([simplify(eval(Cond10)),simplify(eval(Cond00))]); ERROR(`the procedure is implemented for first arguments p,q only`); end: # #---------------------------------------------------------------------------------------- # # does for the operator "a" with the symbol XY(X+Y) # output conditions of factorizations of type (XY)(S) in terms of the invariants of THE full system of invariants # LPDO__Fact_Cond_XY_S:= proc(p,q,i1,i2,i3,i4,i5) local Cond10, Cond00,Ir; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; Cond10:=-q*p*i2-2*diff(p,y)*diff(q,y)*q^3*p-q^3*p*diff(p,x)*diff(p,y)+diff(q,y)*q*p^3*diff(q,x)+2*p^3*diff(q,x)*q*diff(p,x)+diff(q,`$`(y,2))*q^3*p^2-q^4*p*diff(p,`$`(y,2))+q^3*p^2*diff(p,x,y)-q^2*p^3*diff(q,x,y)-diff(p,`$`(x,2))*q^2*p^3+q*p^4*diff(q,`$`(x,2))+2*q^4*diff(p,y)^2+i3*q^3-i4*p^3-2*p^4*diff(q,x)^2; Cond00:= -1/2*1/p*(2*p^5*q^2*diff(q,`$`(x,3))+10*p^4*diff(q,x)*i4-2*i5*q^3+2*q^3*p*diff(i1,x)*diff(p,y)-2*q^3*diff(p,y)*diff(p,x)*i1+3*q^2*diff(q,x)*diff(i1,y)*p^2+q^2*diff(i1,x)*diff(q,y)*p^2+4*q^2*p^3*diff(q,x)^2*diff(p,y)+8*p^2*q*diff(q,x)*i2+2*i4*p^2*diff(p,y)*q^2-2*i3*q^3*p*diff(q,x)-2*q*i4*p^3*diff(q,y)+8*q^2*diff(p,x)^2*p^3*diff(q,x)-28*p^4*diff(q,x)^2*q*diff(p,x)-4*p^4*q*diff(q,x)^2*diff(q,y)-4*diff(p,x)*q*p^3*i4+2*q^4*p*diff(i3,x)-4*p^3*diff(q,x)^2*i1-2*i1*i4*p^2+20*p^5*diff(q,x)^3+4*q^2*p^3*diff(q,x)*diff(q,y)*diff(p,x)-6*q^2*p*diff(q,x)*i1*diff(p,y)+4*p^2*q*diff(q,x)*i1*diff(p,x)-4*q*i1*diff(q,x)*diff(q,y)*p^2-2*p^4*q^3*diff(p,`$`(x,3))-2*p^2*q^2*diff(i2,x)-4*i3*q^4*diff(p,x)-4*p^2*q^3*diff(p,x)*diff(q,x)*diff(p,y)-q^3*diff(i1,x,y)*p^2-2*p^4*q*diff(i4,x)-2*q^3*p^3*diff(q,`$`(x,2))*diff(p,y)-2*q^3*diff(p,`$`(x,2))*p^3*diff(q,y)+2*q^3*p*diff(p,x,y)*i1+q^2*diff(q,x,y)*i1*p^2+2*q^2*diff(q,y)*p^4*diff(q,`$`(x,2))+2*p^2*q^4*diff(p,`$`(x,2))*diff(p,y)+10*p^4*q^2*diff(q,`$`(x,2))*diff(p,x)-2*p^2*q^2* diff(p,`$`(x,2))*i1+2*p^3*q*diff(q,`$`(x,2))*i1-16*p^5*q*diff(q,x)*diff(q,`$`(x,2))-4*p^3*q^3*diff(p,`$`(x,2))*diff(p,x)+10*diff(p,`$`(x,2))*q^2*p^4*diff(q,x)); RETURN([simplify(eval(Cond10)),simplify(eval(Cond00))]); ERROR(`the procedure is implemented for first arguments p,q only`); end: # #---------------------------------------------------------------------------------------- # # does for the operator "a" with the symbol XY(X+Y) # output conditions of factorizations of type (SX)(Y) in terms of the invariants of THE full system of invariants # LPDO__Fact_Cond_SX_Y:= proc(p,q,i1,i2,i3,i4,i5) local Cond1, Cond0,Ir; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; Cond1:=(-i3*q^2+q*p*diff(i1,y)-diff(q,y)*p*i1-2*q*diff(p,y)*i1+2*p*i2)/p^2/q^2; Cond0:=(q^2*p^2*(diff(i1, y, x))-2*q^2*p*i1*(diff(p, y, x))-2*q^2*p*(diff(p, x))*(diff(i1, y))-2*q^2*p*(diff(p, y))*(diff(i1, x))+6*q^2*i1*(diff(p, x))*(diff(p, y))-q*p^2*i1*(diff(q, y, x))-q*p^2*(diff(q, y))*(diff(i1, x))-q*p^2*(diff(q, x))*(diff(i1, y))+2*q*p*i1*(diff(p, x))*(diff(q, y))+2*q*p*i1*(diff(p, y))*(diff(q, x))+2*p^2*i1*(diff(q, y))*(diff(q, x))+2*q*p^2*(diff(i2, x))-4*q*p*(diff(p, x))*i2-4*p^2*i2*(diff(q, x))-2*q^2*i5)/(-2)/q^3/p^3; RETURN([simplify(eval(Cond1)),simplify(eval(Cond0))]); ERROR(`ERROR in LPDO__Fact_Cond_SX_Y`); end: # #--------------------------------------------------------------------------------------- # outputs a set of sets of conditions s.t. input operator is factorable if and # only if one of the sets of conditions holds LPDO__Is_Factorable:=proc(p,q,i1,i2,i3,i4,i5) local resu; resu:=[[`S_XY`, op(LPDO__Fact_Cond_S_XY(p,q,i1,i2,i3,i4,i5))], [`X_SY`, op(LPDO__Fact_Cond_X_SY(p,q,i1,i2,i3,i4,i5))], [`Y_SX`, op(LPDO__Fact_Cond_Y_SX(p,q,i1,i2,i3,i4,i5))], [`SY_X`, op(LPDO__Fact_Cond_SY_X(p,q,i1,i2,i3,i4,i5))], [`XY_S`, op(LPDO__Fact_Cond_XY_S(p,q,i1,i2,i3,i4,i5))], [`SX_Y`, op(LPDO__Fact_Cond_SX_Y(p,q,i1,i2,i3,i4,i5))], [`S_X_Y`,op(LPDO__Fact_Cond_S_X_Y(p,q,i1,i2,i3,i4,i5))], [`S_Y_X`,op(LPDO__Fact_Cond_S_Y_X(p,q,i1,i2,i3,i4,i5))], [`Y_X_S`,op(LPDO__Fact_Cond_Y_X_S(p,q,i1,i2,i3,i4,i5))], [`X_Y_S`,op(LPDO__Fact_Cond_X_Y_S(p,q,i1,i2,i3,i4,i5))], [`X_S_Y`,op(LPDO__Fact_Cond_X_S_Y(p,q,i1,i2,i3,i4,i5))], [`Y_S_X`,op(LPDO__Fact_Cond_Y_S_X(p,q,i1,i2,i3,i4,i5))] ]; RETURN(eval(resu)); end: #--------------------------------------------------------------------------------------- # outputs a set of sets of conditions s.t. input operator is factorable if and # only if one of the sets of conditions holds LPDO__Is_Factorable_via_Coeff:=proc(a) local resu, i1,i2,i3,i4,i5,p,q; if LPDO__nvars() <> 2 then ERROR(` LPDO__Is_Factorable_via_Coeff: illegal number of vars`); fi; if LPDO__degree(a) <> 3 then ERROR(` LPDO__Is_Factorable_via_Coeff: illegal degree of operator`); fi; if LPDO__value(a,[3,0])<>0 or LPDO__value(a,[0,3])<>0 then ERROR(` LPDO__Is_Factorable_via_Coeff: illegal symbol of the operator, works with (pX+qY)XY only`); fi; i1:=LPDO__Inv(1,a); i2:=LPDO__Inv(2,a); i3:=LPDO__Inv(3,a); i4:=LPDO__Inv(4,a); i5:=LPDO__Inv(5,a); p:=LPDO__value(a,[2,1]); q:=LPDO__value(a,[1,2]); resu:=[[S_XY, op(LPDO__Fact_Cond_S_XY(p,q,i1,i2,i3,i4,i5))], [X_SY, op(LPDO__Fact_Cond_X_SY(p,q,i1,i2,i3,i4,i5))], [Y_SX, op(LPDO__Fact_Cond_Y_SX(p,q,i1,i2,i3,i4,i5))], [SY_X, op(LPDO__Fact_Cond_SY_X(p,q,i1,i2,i3,i4,i5))], [XY_S, op(LPDO__Fact_Cond_XY_S(p,q,i1,i2,i3,i4,i5))], [SX_Y, op(LPDO__Fact_Cond_SX_Y(p,q,i1,i2,i3,i4,i5))]]; end: #----------------------------------------------------------- # # compute the invariant number "n" for an operator with the symbol pXXX, # n=1,2,3,4,5 if a02 neq 0 # n=1,..,4, if a02=0 # n=1,..,3, if a02=a11=0 # # LPDO__Inv_XXX :=proc(n::posint,a::array) local p, a20,a11,a02,a10,a01,a00; if n>5 then ERROR(` LPDO__Inv_XXX: illegal index of invariant `); fi; if LPDO__nvars() <> 2 then ERROR(` LPDO__Inv_XXX: illegal number of vars`); fi; if LPDO__degree(a) <> 3 then ERROR(` LPDO__Inv_XXX: illegal degree of operator`); fi; if (LPDO__value(a,[2,1])<>0) or (LPDO__value(a,[1,2])<>0) or (LPDO__value(a,[0,3])<>0) then ERROR(` LPDO__Inv: the second argument should be an operator with the symbol of the form XY(pX+qY)`); fi; p:=LPDO__value(a,[3,0]); a20:=LPDO__value(a,[2,0])/p; a11:=LPDO__value(a,[1,1])/p; a02:=LPDO__value(a,[0,2])/p; a10:=LPDO__value(a,[1,0])/p; a01:=LPDO__value(a,[0,1])/p; a00:=LPDO__value(a,[0,0])/p; # general case a02 <> 0 if a02<>0 then print(`Case a02<>0. Another cases are a02=0 and a11<>0 or a02=a11=0. In these cases invariants are different!`); if (n=1) then RETURN(eval(a02)); elif (n=2) then RETURN(eval(a11)); elif (n=3) then RETURN(simplify(eval((6*a10*a02-2*a20^2*a02-3*a11*a01+a11^2*a20-6*diff(a20,LPDO__vars[1])*a02)/a02/6))); elif (n=4) then RETURN(simplify(eval(-(-3*diff(a01,LPDO__vars[1])*a02+diff(a11,LPDO__vars[1])*a20*a02+a11*diff(a20,LPDO__vars[1])*a02+3*diff(a02,LPDO__vars[1])*a01-diff(a02,LPDO__vars[1])*a11*a20+2*a02^2*diff(a20,LPDO__vars[2]))/a02/3))); elif (n=5) then RETURN(simplify(eval(-(-108*a00*a02+36*a10*a20*a02+27*a01^2-18*a01*a11*a20-8*a20^3*a02+18*a11*diff(a20,LPDO__vars[2])*a02+3*a11^2*a20^2+54*diff(a01,LPDO__vars[2])*a02-18*diff(a11,LPDO__vars[2])*a20*a02-54*diff(a02,LPDO__vars[2])*a01+18*diff(a02,LPDO__vars[2])*a11*a20+36*diff(a20,`$`(LPDO__vars[1],2))*a02)/a02/108))); else ERROR(`the first argument should be 1,..,5`); fi; elif a02=0 and a11<>0 then print(`Case a02=0 and a11<>0. Another cases are a02<>0 and a02=a11=0. In these cases invariants are different!`); if (n=1) then RETURN(eval(a11)); elif (n=2) then RETURN(eval(a01-1/3*a11*a20)); elif (n=3) then RETURN(simplify(eval(diff(a10,LPDO__vars[1])-2/3*a20*diff(a20,LPDO__vars[1])-1/a11*diff(a11,LPDO__vars[1])*a10+1/3*1/a11*diff(a11,LPDO__vars[1])*a20^2+1/a11*diff(a11,LPDO__vars[1])*diff(a20,LPDO__vars[1])-1/3*a11*diff(a20,LPDO__vars[2])-diff(a20,`$`(LPDO__vars[1],2))))); elif (n=4) then RETURN(simplify(eval(a00-1/a11*a01*a10+1/3*1/a11*a01*a20^2+1/a11*a01*diff(a20,LPDO__vars[1])-1/27*a20^3-1/3*a11*diff(a20,LPDO__vars[2])-1/3*a20*diff(a20,LPDO__vars[1])-1/3*diff(a20,`$`(LPDO__vars[1],2))))); else ERROR(`the generating set consists of four invariants In, where n=1,..,4`); fi; elif a02=0 and a11=0 and a01<>0 then print(`Case a02=a11=0. In other cases invariants are different!!`); if (n=1) then RETURN(eval(a01)); elif (n=2) then RETURN(eval(a10-1/3*a20^2-diff(a20,LPDO__vars[1]))); elif (n=3) then RETURN(simplify(eval(diff(a00,LPDO__vars[1])-1/3*diff(a10,LPDO__vars[1])*a20-1/3*a10*diff(a20,LPDO__vars[1]) -1/a01*diff(a01,LPDO__vars[1])*a00+1/3*1/a01*diff(a01,LPDO__vars[1])*a10*a20 -2/27*1/a01*diff(a01,LPDO__vars[1])*a20^3+1/3*1/a01*diff(a01,LPDO__vars[1])*diff(a20,`$`(LPDO__vars[1],2)) -1/3*a01*diff(a20,LPDO__vars[2])+2/9*a20^2*diff(a20,LPDO__vars[1])-1/3*diff(a20,`$`(LPDO__vars[1],3))))); else ERROR(`the generating set consists of three invariants In, where n=1,..,3`); fi; else if (n=1) then RETURN( expand(eval(a10-1/3*a20^2-diff(a20,LPDO__vars[1])))); else RETURN(expand(eval( a00-1/3*a10*a20+2/27*a20^3-1/3*diff(a20,LPDO__vars[1],LPDO__vars[1]) ))); fi; fi; end: #----------------------------------------------------------- # # compute the invariant number "n" for an operator with the symbol pXXY, # n=1,2,3,4,5 # # LPDO__Inv_XXY := proc(n::posint,a::array) local p,a20,a11,a02,a10,a01,a00; if n>5 then ERROR(` LPDO__Inv_XXX: illegal index of invariant `); fi; if LPDO__nvars() <> 2 then ERROR(` LPDO__Inv_XXX: illegal number of vars`); fi; if LPDO__degree(a) <> 3 then ERROR(` LPDO__Inv_XXX: illegal degree of operator`); fi; if (LPDO__value(a,[3,0])<>0) or (LPDO__value(a,[1,2])<>0) or (LPDO__value(a,[0,3])<>0) then ERROR(` LPDO__Inv: the second argument should be an operator with the symbol XXY`); fi; p:=LPDO__value(a,[2,1]); a20:=LPDO__value(a,[2,0])/p; a11:=LPDO__value(a,[1,1])/p; a02:=LPDO__value(a,[0,2])/p; a10:=LPDO__value(a,[1,0])/p; a01:=LPDO__value(a,[0,1])/p; a00:=LPDO__value(a,[0,0])/p; if (n=1) then RETURN(eval(a02)); elif (n=2) then RETURN(eval(diff(a11,LPDO__vars[2])-2*diff(a20,LPDO__vars[1]))); elif (n=3) then RETURN(simplify(eval(a10-a20*a11-2*diff(a20,LPDO__vars[1])))); elif (n=4) then RETURN(simplify(eval(a01-1/4*a11^2-2*a02*a20-1/2*diff(a11,LPDO__vars[1])))); elif (n=5) then RETURN(simplify(eval(a00-1/2*a10*a11-a01*a20+1/2*a20*a11^2-a02*diff(a20,LPDO__vars[2])+a02*a20^2-diff(a20,`$`(LPDO__vars[1],2))))); else ERROR(`the first argument should be 1,..,5`); fi; end: #----------------------------------------------------------------- # # compute the invariant number "n" for an operator with the symbols LPDO__XY(X+qY) and pXXX and pXXY # LPDO__Inv := proc(n::posint,a::array) local res; if n>7 then ERROR(` LPDO__Inv: illegal index of invariant `); fi; if LPDO__nvars() <> 2 then ERROR(` LPDO__Inv: illegal number of vars`); fi; if LPDO__degree(a) <> 3 then ERROR(` LPDO__Inv: illegal degree of operator`); fi; ## CASE of the SYMBOL XXX: if (LPDO__value(a,[3,0])<>0) and (LPDO__value(a,[2,1])=0) and (LPDO__value(a,[1,2])=0) and (LPDO__value(a,[0,3])=0) then RETURN(LPDO__Inv_XXX(n,a)); fi; ## CASE of the SYMBOL XXY: if (LPDO__value(a,[2,1])<>0) and (LPDO__value(a,[3,0])=0) and (LPDO__value(a,[1,2])=0) and (LPDO__value(a,[0,3])=0) then RETURN(LPDO__Inv_XXY(n,a)); fi; ## CASE of the SYMBOL XY(pX+qY): if (LPDO__value(a,[3,0])=0) and (LPDO__value(a,[0,3])=0) then RETURN(LPDO__Inv_XYS(n,a)); fi; ERROR(` LPDO__Inv: the symbol of the operator should be XXX,XXY or XY(pX+qY)`); end: #----------------------------------------------------------------- # # Necessary and sufficient invariant condition of factorization with the factorization type (X)(XY) # # LPDO__Fact_Cond_X_XY:= proc(i1,i2,i3,i4,i5) local x,y,u; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; x:=LPDO__vars[1]; y:=LPDO__vars[2]; if LPDO__warnings=true then print(`case i3=0 in LPDO__Fact_Cond_X_XY requires that`, 0=i4 + u(x,y)^2+ diff(u(x,y),x),`has solution w.r.t. u`): fi: if i3=0 then RETURN([i1,i3,i5]); else RETURN([[case_i3_isnot_0, [i1,simplify(eval(i4-3/i3^2*i5*diff(i3,x)+1/i3*diff(i5,x) -1/i3*diff(i3,`$`(x,2))+1/i3^2*i5^2+2/i3^2*diff(i3,x)^2))]], [case_i3=0, [i1,i3,i5]]]); fi; end: #----------------------------------------------------------------- # # Necessary and sufficient invariant condition of factorization with the factorization type (X)(XY) # checked # Assuming r1[x] = -i4-(1/4)*a11^2+(1/2)*a11[x]+r1*a11-r1^2 can be solved for r1. LPDO__Fact_Cond_X_X_Y:= proc(i1,i2,i3,i4,i5) local x,y; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; x:=LPDO__vars[1]; y:=LPDO__vars[2]; RETURN([simplify(eval(i1)),simplify(eval(i3)),simplify(eval(i5))]); end: #----------------------------------------------------------------- # # Necessary and sufficient invariant condition of factorization with the factorization type (X)(XY) # checked # # LPDO__Fact_Cond_XY_X:= proc(i1,i2,i3,i4,i5) local x,y,u; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; x:=LPDO__vars[1]; y:=LPDO__vars[2]; if LPDO__warnings=true then if i2=i3 then print(`LPDO__Fact_Cond_XY_X requires also that`, 0=i4-2*diff(i1,y)+u(x,y)^2 + diff(u(x,y),x),`has solution w.r.t. u`): else print(`case i2=i3 in LPDO__Fact_Cond_XY_X requires also that`, 0=i4-2*diff(i1,y)+u(x,y)^2 + diff(u(x,y),x),`has solution w.r.t. u`): fi: fi: if i2=i3 then RETURN([i1,simplify(eval(diff(i2,x)-2*i5+2*diff(i4,y)))]); # provided r[x] = i4+(1/4)*a11^2+(1/2)*a11[x]-r*a11+r^2 can be solved for r else RETURN([[`case i2=i3`,[i1,simplify(eval(diff(i2,x)-2*i5+2*diff(i4,y)))]],[`case i2<>i3`,[i1,simplify(eval(1/4*(4*diff(i1,`$`(y,2))*diff(i3,x)-8*diff(i1,`$`(y,2))*diff(i2,x)+4*diff(i4,x,y)*i3-4*diff(i4,x,y)*i2 -4*diff(i1,x,`$`(y,2))*i3+4*diff(i1,x,`$`(y,2))*i2+2*diff(i2,`$`(x,2))*i3-2*diff(i2,`$`(x,2))*i2-8*diff(i4,y)*diff(i1,`$`(y,2)) +8*diff(i1,`$`(y,2))*i5+4*diff(i1,`$`(y,2))^2+16*diff(i1,y)*i3*i2-8*diff(i1,y)*i3^2+4*i4*i2^2+4*i4*i3^2-8*diff(i1,y)*i2^2 -2*diff(i3,x)*diff(i2,x)-4*diff(i4,y)*diff(i3,x)+8*diff(i4,y)*diff(i2,x)+4*i5*diff(i3,x)-8*i5*diff(i2,x)-4*diff(i5,x)*i3+4*diff(i5,x)*i2 -8*diff(i4,y)*i5+3*diff(i2,x)^2+4*diff(i4,y)^2+4*i5^2-8*i4*i3*i2)/(i3-i2)^2))]]]); fi: end: #----------------------------------------------------------------- # # Necessary and sufficient invariant condition of factorization with the factorization type (Y)(XX) # # LPDO__Fact_Cond_Y_XX:= proc(i1,i2,i3,i4,i5) local x,y; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; x:=LPDO__vars[1]; y:=LPDO__vars[2]; RETURN([simplify(eval(-i5+diff(i4,y)-diff(i1,`$`(y,2))+1/2*diff(i2,x)))]); end: #----------------------------------------------------------------- # # Necessary and sufficient invariant condition of factorization with the factorization type (Y)(XX) # # LPDO__Fact_Cond_XX_Y:= proc(i1,i2,i3,i4,i5) local x,y; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; x:=LPDO__vars[1]; y:=LPDO__vars[2]; RETURN([simplify(eval(i3)), simplify(eval(i5))]); end: #----------------------------------------------------------------- # # Necessary and sufficient invariant condition of factorization with the factorization type (X)(XX) # # LPDO__Fact_Cond_X_XX_a02_0:= proc(i1,i2,i3,i4) local x,y; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; x:=LPDO__vars[1]; y:=LPDO__vars[2]; RETURN([simplify(eval((-i2^3+i4*i1^3+6*diff(i1,x)^3-3*i2*i1*diff(i2,x)+4*i2*i1*diff(i1,`$`(x,2))+5*diff(i1,x)*i1*diff(i2,x)-6*diff(i1,x)*i1 *diff(i1,`$`(x,2))+6*diff(i1,x)*i2^2-11*i2*diff(i1,x)^2-i1^2*diff(i2,`$`(x,2))+i1^2*diff(i1,`$`(x,3)))/i1^3-i3))]); end: #----------------------------------------------------------------- # # Necessary and sufficient invariant conditions of # LPDOs of the eq. class defined by the values I1, ... ,I5 of the # generating set of invaraints to be factorable with the factorization type (Var2)(Var1^2) # # LPDO__Fact_Cond_Var2_Var1Var1:= proc(I1,I2,I3,I4,I5) if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; RETURN([simplify(eval(I3-I2)),simplify(eval(-I5+diff(I4,LPDO__vars[2])-diff(I1,`$`(LPDO__vars[2],2))+1/2*diff(I2,LPDO__vars[1])))]); end: #------------------------------------------------------------------ # # Invariants for operators of the form L=Dxx + aDx + b Dy + c , b <>0 # LPDO__Inv_strict_parabolic_ord2:=proc(n,L) local a,b,c,x,y; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; if b=0 then ERROR(`LPDO__Inv_strict_parabolic: the coefficient at Dy must be non zero`); fi; if LPDO__degree(L) <> 2 then ERROR(` LPDO__Inv_strict_parabolic: illegal degree of operator`); fi; if n>3 then ERROR(`LPDO__Inv_strict_parabolic: illegal index of invariant `); fi; if (LPDO__value(L,[0,2])<>0) or (LPDO__value(L,[1,1])<>0) then ERROR(` LPDO__Inv_strict_parabolic: the second argument should be an operator with the symbol XXY`); fi; x:=LPDO__vars[1]; y:=LPDO__vars[2]; a:=LPDO__value(L,[1,0]); b:=LPDO__value(L,[0,1]); c:=LPDO__value(L,[0,0]); if (n=1) then RETURN(eval(b)); else RETURN(eval(diff(c,x)-1/2*a*diff(a,x)-1/b*diff(b,x)*c+1/4*1/b*diff(b,x)*a^2+1/2*1/b*diff(b,x)*diff(a,x)-1/2*b*diff(a,y)-1/2*diff(a,x,x))); fi; end: #----------------------------------------------------------- # computes the symbol of an operator # LPDO__Symbol:=proc(L::array) local d, i, S; d:=LPDO__degree(L); S:=0; for i from 0 to d do S:=S+LPDO__value(L,[d-i,i])*X^(d-i)*Y^i; end do: RETURN(eval(S)); end: # #-------------------------------------------------------------------------------- # Darboux transformations theory: #-------------------------------------------------------------------------------- # Darboux transformation. Given L,z1,z2, find M,M1,L1 # Let L - laplacian operator and z1,z2 - linear independent solutions # Find the genralized Laplace transformation: # # M*L = L1*M1 #where (with coefficient-multiplicator) # | z z_x z_y | #M1(z) = | z1 z1_x z1_y | # | z2 z2_x z2_y | # # returns [M,L1,M1] Darboux:= proc(L,z1,z2) local a,b,c,z1x,z1y,z2x,z2y,p0,q,qx,r1,r, a1,b1,c1,M1,M,L1; a := LPDO__value(L,[1,0]); b := LPDO__value(L,[0,1]); c := LPDO__value(L,[0,0]); z1x := diff(eval(z1),x): z1y := diff(eval(z1),y): z2x := diff(eval(z2),x): z2y := diff(eval(z2),y): p0 := -z1 *z2y + z2 *z1y; if p0=0 then ERROR(`Determinant formulae imply a Laplace Y-transformation, use procedure DarbouxX`); fi; q := simplify((z1 *z2x - z2 *z1x)/p0); qx := diff(q,x); r1 := (z1x*z2y - z1y*z2x)/p0; r := simplify(r1 - diff(q,x)/q + diff(q,y)); a1 := a; b1 := simplify(b - diff(q,x)/q); c1 := simplify(c + diff(b,y) + (diff(b,x)-a*qx-diff(r,x))/q + (qx*(r-b)-diff(qx,x))/q^2 + 2*qx^2/q^3); M := LPDO__create1(1,q,r): M1 := LPDO__create1(1,q,r1): L1 := LPDO__create2(0,1,0,a1,b1,c1): RETURN([eval(M),eval(L1), eval(M1)]); end: #--------------------------------------------------------------------------------- # Darboux X-transformation. Given L,z1, find M,M1,L1 # Let L - laplacian operator and z1 - solution of L # Find the genralized Laplace transformation: # M*L = L1*M1 # returns [M,L1,M1] # DarbouxX := proc(L,z1) local a,b,c,z1x,z1y,z1yy,ay, a1, c1, F, M1,M,L1; a := LPDO__value(L,[1,0]); b := LPDO__value(L,[0,1]); c := LPDO__value(L,[0,0]); z1x := diff(z1,x): z1y := diff(z1,y): z1yy:= diff(z1y,y): ay := diff(a,y): M1 := LPDO__create1(0,1, -z1y/z1): M := LPDO__create1(0,1, -(z1yy+a*z1y+ay*z1)/(z1y+a*z1) ): L1 := LPDO__create2(0,1,0,a1(x,y), b,c1(x,y)): F := LPDO__minus(LPDO__mult(M,L),LPDO__mult(L1,M1)): #LPDO__print(F); a1(x,y) := factor(solve(LPDO__value(F,[1,0]), a1(x,y))); #F := LPDO__simplify(F): LPDO__print(F); #print(`a1=`, a1(x,y)); c1(x,y) := factor(solve(LPDO__value(F,[0,0]), c1(x,y))); #print(`c1=`, c1(x,y)); #F := LPDO__simplify(F): LPDO__print(F); L1 := LPDO__simplify(L1): M1 := LPDO__simplify(M1): M := LPDO__simplify(M): RETURN([eval(M),eval(L1), eval(M1)]); end: #---------------------------------------------------------------------------------- # Darboux Y-transformation. Given L,z1, find M,M1,L1 # Let L - laplacian operator and z1 - solution of L # Find the genralized Laplace transformation: # M*L = L1*M1 # returns [M,L1,M1] # DarbouxY := proc(L,z1) local a,b,c,z1x,z1y,z1xx,bx, c1, F, b1, M1,M,L1,sol; a := LPDO__value(L,[1,0]); b := LPDO__value(L,[0,1]); c := LPDO__value(L,[0,0]); z1x := diff(z1,x): z1y := diff(z1,y): z1xx:= diff(z1x,x): bx := diff(b,x): M1 := LPDO__create1(1,0, -z1x/z1): M := LPDO__create1(1,0, b1(x,y)-b-z1x/z1 ): L1 := LPDO__create2(0,1,0,a, b1(x,y),c1(x,y)): F := LPDO__minus(LPDO__mult(M,L),LPDO__mult(L1,M1)): sol:=solve({LPDO__value(F,[1,0]),LPDO__value(F,[0,0])},[c1(x,y),b1(x,y)]); assign(sol); # F := LPDO__simplify(F): L1 := LPDO__simplify(L1): M1 := LPDO__simplify(M1): M := LPDO__simplify(M): RETURN([eval(M),eval(L1), eval(M1)]); end: #----------------------------------------------------------- # given list of coefficients, creates an LPDO Create_LPDO_from_coeff:= proc(a::list) local res,i,Coeff_in_list_form; res:=LPDO__create(); for i from 0 to nops(a)-1 do Coeff_in_list_form:=inttovec(i,LPDO__nvars()); LPDO__add_value(res, a[i+1] , Coeff_in_list_form): end do; eval(res); end: #----------------------------------------------------------- # DIFF(f,[i1, ..., i_n])= D_{x1}^{i1} ... D_{xn}^{in} (f) # DIFF:= proc(f,lis::list) local i,j,res; if LPDO__nvars()<>nops(lis) then ERROR(`the length of the second argument is`,nops(lis),`while it must be equal to the number of variables`, LPDO__nvars()); fi; res:=f; for i from 1 to nops(lis) do for j from 1 to lis[i] do res:=diff(res,LPDO__vars[i]); end do: end do: RETURN(res); end: #------------------------------------------------------------ # returns substitutions needed for computation of generating sets of invariants # my means of moving frames #f Subs_list:=proc(G) local res,x,y,z; if LPDO__nvars() > 3 then ERROR(`the procedure works for 2 and 3 varaibles only`); fi; if LPDO__nvars() = 1 then ERROR(`the procedure works for 2 and 3 varaibles only`); fi; if LPDO__nvars() = 3 then x:=LPDO__vars[1]; y:=LPDO__vars[2]; z:=LPDO__vars[3]; res:=[[diff(G,x)=Gx, diff(G,y)=Gy, diff(G,z)=Gz], [diff(G,x,y)=Gxy,diff(G,x,z)=Gxz,diff(G,y,z)=Gyz,diff(G,x,x)=Gxx,diff(G,y,y)=Gyy,diff(G,z,z)=Gzz], [diff(G,x,x,x)=Gxxx,diff(G,x,x,y)=Gxxy,diff(G,x,x,z)=Gxxz, diff(G,x,y,y)=Gxyy,diff(G,x,y,z)=Gxyz,diff(G,x,z,z)=Gxzz, diff(G,y,y,y)=Gyyy,diff(G,y,y,z)=Gyyz,diff(G,z,z,z)=Gzzz]]; else x:=LPDO__vars[1]; y:=LPDO__vars[2]; res:=[[diff(G,x)=Gx, diff(G,y)=Gy], [diff(G,x,y)=Gxy,diff(G,x,x)=Gxx,diff(G,y,y)=Gyy], [diff(G,x,x,x)=Gxxx,diff(G,x,x,y)=Gxxy, diff(G,x,y,y)=Gxyy, diff(G,y,y,y)=Gyyy]]; fi: RETURN(eval(res)); end: #---------------------------------------------------------- # Projection from the left # Returns M-A*L without DxDy LPDO__pl := proc(L, M) local deg, ix, iy, M1, M2, co; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); end if; if LPDO__degree(L) <> 2 then ERROR(`illegal degree of operator`); end if; if LPDO__value(L, [2,0])<>0 or LPDO__value(L, [1,1])<>1 or LPDO__value(L, [0,2])<>0 then ERROR(`operator should have the symbol of the form XY or X^2 + Y^2`); end if; M1 := LPDO__simplify(M); deg := LPDO__degree(M1); if deg<2 then RETURN(M); end if; for ix from deg by -1 to 1 do for iy from deg-ix by -1 to 1 do co := LPDO__value(M1, [ix, iy]); if co<> 0 then #print(778, ix, iy); M2 := LPDO__create(); LPDO__add_value(M2, -co, [ix-1,iy-1]);# LPDO__print(M2); M1 := LPDO__add(M1, LPDO__mult(M2, L));# LPDO__print(M1); end if; end do; end do; eval(M1); end: #------------------------------------------- # Projection from the right # Returns M-L*A without DxDy LPDO__pr := proc(L, M) local deg, ix, iy, M1, M2, co; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); end if; if LPDO__degree(L) <> 2 then ERROR(`illegal degree of operator`); end if; if LPDO__value(L, [2,0])<>0 or LPDO__value(L, [1,1])<>1 or LPDO__value(L, [0,2])<>0 then ERROR(`operator should have the symbol of the form XY or X^2 + Y^2`); end if; M1 := LPDO__simplify(M); deg := LPDO__degree(M1); if deg<2 then RETURN(M); end if; for ix from deg by -1 to 1 do for iy from deg-ix by -1 to 1 do co := LPDO__value(M1, [ix, iy]); if co<> 0 then #print(778, ix, iy); M2 := LPDO__create(); LPDO__add_value(M2, -co, [ix-1,iy-1]);# LPDO__print(M2); M1 := LPDO__add(M1, LPDO__mult(L, M2));# LPDO__print(M1); end if; end do; end do; eval(M1); end: #--------------------------------------------- # Returns the gauged evolution invariants for pair (L,M), where L=D_x D_x D_y + .. and M=D_x+m or M=D_y+m LPDO__pairs_invariants3:=proc(L::array,M::array) local a00,a10,a01,a20,a11,a02,I0,I1,I2,I3,I4,m; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars, currently implemented for two variables only`); end if; if LPDO__degree(L) <> 3 then ERROR(`illegal degree of first operator, this procedure works for the first operator in the form D_x D_x D_y + ..`); end if; if LPDO__value(L, [3,0])<>0 or LPDO__value(L, [1,2])<>0 or LPDO__value(L, [0,3])<>0 or LPDO__value(L, [2,1])<>1 then ERROR(`the first operator should be in the form in the form D_x D_x D_y + ...`); end if; if LPDO__degree(M) <> 1 then ERROR(`this procedure works for the second operator in the form D_x + m or D_y + m only`); end if; if LPDO__value(M, [1,0])=1 and LPDO__value(M, [0,1])=1 then ERROR(`this procedure works for the second operator in the form D_x + m or D_y + m only`); end if; m:=LPDO__value(M, [0,0]); a00:=LPDO__value(L, [0,0]): a10:=LPDO__value(L, [1,0]): a01:=LPDO__value(L, [0,1]): a20:=LPDO__value(L, [2,0]): a11:=LPDO__value(L, [1,1]): a02:=LPDO__value(L, [0,2]): if LPDO__value(M, [1,0])=1 then I0:=a02; I1:=a11-2*m; I2:=diff(m,y)-diff(a20,x); I3:=a01-a11*m-2*a02*a20-diff(m,x)+m^2; I4:=(-a10+a11*a20+2*(diff(m, y)))*m+a00+a02*a20^2-a01*a20-(diff(m, y))*a11-a02*(diff(a20, y))-(diff(m, y, x)): else I0:=a02; I1:=a20-m; I2:=diff(a11, y)-2*(diff(m, x)); I3:=a10-a20*a11-2*(diff(m, x)); I4:=a02*m^2+((1/4)*a11^2-a01+(1/2)*(diff(a11, x)))*m+a00-(1/2)*a10*a11-a02*(diff(m, y))-(1/2)*a20*(diff(a11, x))+(1/4)*a20*a11^2-diff(m, x, x); fi: RETURN([simplify(I0),simplify(I1),simplify(I2),simplify(I3),simplify(I4)]); end: #-------------------------------------------- # finds pairs invariants w.r.t. gauge transformations of the pairs: m,h,q, R and w.r.t. to gauged evolution: I1, I2, and I3. LPDO__pairs_invariants:=proc(L,M) local a,b,q,r,R,h,m,I2,I3; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars, currently implemented for two variables only`); end if; if LPDO__degree(L) <> 2 then ERROR(`illegal degree of first operator, currently implemented for oeprator of order two only`); end if; if LPDO__value(L, [2,0])<>0 or LPDO__value(L, [1,1])<>1 or LPDO__value(L, [0,2])<>0 then ERROR(`operator should have the symbol of the form XY`); end if; if LPDO__degree(M) <> 1 then ERROR(`illegal degree of second operator, currently implemented for oeprator of order 1 only`); end if; a:=LPDO__value(L, [1,0]): b:=LPDO__value(L, [0,1]): if LPDO__value(M, [1,0])<>1 then ERROR(`the coefficient at D_x of the second operator must be 1`); end if; q:=LPDO__value(M, [0,1]): r:=LPDO__value(M, [0,0]): R:=eval(r-b-q*a); h:=eval(LPDO__Lapl_h(L)); m:=eval(h-LPDO__Lapl_k(L)); I2:=simplify(2*m+diff(R/q,x) - diff(R,y)); I3:=simplify(2*h-R^2/2/q+diff(R/q,x)); RETURN([[m,h,q,R],[q,I2,I3]]); end: #-------------------------------------------- # returns conditions for the existence of a Darboux transformation in terms of given invariants I1,I2,I3 of pair M, L w.r.t. gauged evolution. LPDO__DT_conditions_in_gauged_evolution_invariants:=proc(q,I2,I3) local F1,F2; if LPDO__nvars() <> 2 then ERROR(`this procedure works only for two independent variables`); end if; F2:=simplify(eval(diff(I3, x) -(-(-3*(diff(q, x))*(diff(q, y, x))*q+q^4*(diff(I3, y))+(diff(q, y))*q^3*I3-(diff(q, x))*I3*q^2+(diff(q, y, x, x))*q^2+3*(diff(q, y))*(diff(q, x))^2-q*(diff(q, y))*(diff(q, x, x)))/q^3))): F1:=simplify(eval(I2-(-(-(diff(q, y))*(diff(q, x))+(diff(q, y, x))*q)/q^2))); RETURN([F1,F2]); end: #-------------------------------------------- # finds pairs invariants w.r.t. to gauged evolution: I1, I2, and I3 in terms of # invariants of pairs w.r.t. gauge transformations of the pairs: m,h,q, R and LPDO__pairs_invariants_evolution_via_gauged:=proc(m,h,q,R) local I2,I3; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars, currently implemented for two variables only`); end if; if q=0 then ERROR(`third invariant, invariant q cannot be zero`); end if; if LPDO__degree(M) <> 1 then ERROR(`illegal degree of second operator, currently implemented for oeprator of order 1 only`); end if; if LPDO__value(M, [1,0])<>1 then ERROR(`the coefficient at D_x of the second operator must be 1`); end if; I2:=simplify(2*m+diff(R/q,x) - diff(R,y)); I3:=simplify(2*h-R^2/2/q+diff(R/q,x)); RETURN([q,I2,I3]); end: #------------------------------------------- # returns joint gauge invariants for operators L=D_xy + aD_x +bD_y +c and M=m_d D_x^d + ... m_{-d} D_y^d # M must be without any mixed derivatives # returns [m,h,[R_1,R_2, ... R_d],[R_{-1}, .., R_{-d}]] LPDO__joint_gauge_invariants:=proc(L::array,M::array) local d, coeff_x,coeff_y,coeff_0,i,R_x,R_y,R_0,R_current,R_current_y,w,j,a,b,forBell,forBell_y,m,h; if LPDO__nvars() <> 2 then ERROR(`this procedure works only for two independent variables`); end if; if LPDO__value(L, [2,0])<>0 or LPDO__value(L, [1,1])<>1 or LPDO__value(L, [0,2])<>0 then ERROR(`operator should have the symbol of the form XY`); end if; d:=LPDO__degree(M); coeff_x:=[]: coeff_y:=[]: coeff_0:=LPDO__value(M,[0,0]): for i from 1 to d do coeff_x:=[op(coeff_x),LPDO__value(M,[i,0])]; coeff_y:=[op(coeff_y),LPDO__value(M,[0,i])]; od: R_x:=[]; R_y:=[]; a:=LPDO__value(L,[1,0]): b:=LPDO__value(L,[0,1]): for j from 1 to d do R_current:=0; R_current_y:=0; for w from j to d do forBell:=[-b]; forBell_y:=[-a]; for i from 1 to w-j-1 do forBell:=[op(forBell),diff(-b, [x$i])]: forBell_y:=[op(forBell_y),diff(-a, [y$i])]: od: R_current:=R_current+coeff_x[w]*binomial(w,j)*CompleteBellB(w-j, op(forBell)): R_current_y:=R_current_y+coeff_y[w]*binomial(w,j)*CompleteBellB(w-j, op(forBell_y)): od: R_x:=[op(R_x),R_current]; R_y:=[op(R_y),R_current_y]; od: R_0:=coeff_0; for w from 1 to d do forBell:=[-b]; forBell_y:=[-a]; for i from 1 to w-1 do forBell:=[op(forBell),diff(-b, [x$i])]: forBell_y:=[op(forBell_y),diff(-a, [y$i])]: od: R_0:=R_0+coeff_x[w]*CompleteBellB(w, op(forBell))+coeff_y[w]*CompleteBellB(w, op(forBell_y)): od: h:=LPDO__Lapl_h(L); m:=LPDO__Lapl_m(L); RETURN([m,h,R_x,R_y,R_0]); end: #################################################################################################################### # # User VERY Friendly Functions #################################################################################################################### # # ginv(L) returns a set of generating differential invariants for L with respect to the gauge transformations # ginv(L,M) returns a set of generating joint differential invariants for L and M with respect to the gauge transformations # ginv:=proc(L::array:=LPDO__create0(0), M::array:=LPDO__create0(0)) local i,j,counter,res,degL,degM; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars, currently implemented for two variables only`); end if: degL:=LPDO__degree(L); degM:=LPDO__degree(M); if degL=0 then ERROR(`your first argument must be an LPDO of order larger than 0`); fi: if degM=0 then # we are asked to find generating invariants for an operator if degL> 3 or degL=1 then ERROR(`not implemented to search invariants of the individual operators of orders larger than 3 or of order 1`); fi: res:=[]: #print(`individual operator case`); if degL=3 then ## CASE of the SYMBOL XXX: if (LPDO__value(L,[3,0])<>0) and (LPDO__value(L,[2,1])=0) and (LPDO__value(L,[1,2])=0) and (LPDO__value(L,[0,3])=0) then for i from 1 to 5 do res:=[op(res),LPDO__Inv_XXX(i,L)]: od: fi; ## CASE of the SYMBOL XXY: if (LPDO__value(L,[2,1])<>0) and (LPDO__value(L,[3,0])=0) and (LPDO__value(L,[1,2])=0) and (LPDO__value(L,[0,3])=0) then for i from 1 to 5 do res:=[op(res),LPDO__Inv_XXY(i,L)]: od: fi: ## CASE of the SYMBOL XY(pX+qY): if (LPDO__value(L,[3,0])=0) and (LPDO__value(L,[0,3])=0) and (LPDO__value(L,[2,1])<>0) and (LPDO__value(L,[2,1])<>0) then for i from 1 to 7 do res:=[op(res),LPDO__Inv_XYS(i,L)]: od: fi; RETURN(map(el->simplify(eval(el)),res)); else ## case of operators of order 2 if LPDO__value(L, [2,0])=0 and LPDO__value(L, [1,1])=1 and LPDO__value(L, [0,2])=0 or LPDO__value(L, [2,0])=1 and LPDO__value(L, [1,1])=0 and LPDO__value(L, [0,2])=1 then RETURN([eval(LPDO__Lapl_h(L)),eval(LPDO__Lapl_k(L))]): else if LPDO__value(L,[1,0])=0 then RETURN([]): else RETURN([`these invariants exist if b<>0`, [eval(LPDO__Inv_strict_parabolic_ord2(1,L)),eval(LPDO__Inv_strict_parabolic_ord2(2,L))]]): fi: fi: fi; else # we are asked to find joint generating invariants for a pair of operators # case L=Dxy + aDx + bDy + c if degL=2 and LPDO__value(L, [2,0])=0 and LPDO__value(L, [1,1])=1 and LPDO__value(L, [0,2])=0 then if degM=1 and LPDO__value(M, [1,0])=1 and LPDO__value(M, [0,1])<>0 then RETURN(map(el->simplify(eval(el)),LPDO__pairs_invariants(L,M)[1])): else counter:=0: for i from 1 to degM do for j from 1 to degM do counter:=LPDO__value(M,[i,j])+counter: od: od: if counter<>0 then ERROR(`the second argument, an LPDO cannot have mixed derivatives`); fi: RETURN(map(el->simplify(eval(el)),LPDO__joint_gauge_invariants(L,M))): fi: fi: # case L=Dxxy + ... M=D_x+m or M=D_y+m if degL=3 and degM=1 and LPDO__value(L, [3,0])=0 and LPDO__value(L, [2,1])=1 and LPDO__value(L, [1,2])=0 and LPDO__value(L, [0,3])=0 then if LPDO__value(M, [1,0])=0 or LPDO__value(M, [0,1])=0 then RETURN(map(el->simplify(eval(el)),LPDO__pairs_invariants3(L,M))): fi: fi: fi: ERROR(`"ginv" is not implemented for given operators`,L,`and`,M); end: #################################################################################################################################### ######################################################################### # # performs *INVERTIBLE* DT for L with symbol XXY with M=Dy+m on invariants # LPDO__XXY_DT_with_Dy:=proc(b,I2,I3,I4,r) local J1,J2,J3,J4,J5,x,y; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; x:=LPDO__vars[1]; y:=LPDO__vars[2]; if I3<>0 then ERROR(`Invertible (!!) DT with dy exists only if I3=0`); fi; J1:=b; J2:=I2-(2*((diff(r, x))*(diff(r, y))-(diff(r, y, x))*r))/r^2; J3:=J2; J4:=I4+((diff(r, y))*b+(diff(b, y))*r)/r; J5:=r+diff(I2,x)/2+diff(I4,y)+(diff(r, y, y))*b/r-(diff(r, y))^2*b/r^2+(diff(r, y))*(diff(b, y))/r-(diff(r, x, x))*(diff(r, y))/r^2+(diff(r, y, x, x))/r+2*(diff(r, x))^2*(diff(r, y))/r^3-2*(diff(r, x))*(diff(r, y, x))/r^2; RETURN([J1,J2,J3,J4,J5]); end: ########################################################################## # # performs *INVERTIBLE* DT for L with symbol XXY with M=Dx+m on invariants, # here L=CM+r and I2 <> I3 # LPDO__XXY_DT_with_Dx:=proc(i1,i2,i3,i4,i5,r) local J1,J2,J3,J4,J5,x,y; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; x:=LPDO__vars[1]; y:=LPDO__vars[2]; if i1<>0 then ERROR(`Invertible (!!) DT with dx exists only if I1=0`); fi; if i2=i3 then ERROR(`This procedure is not working for case I2=I3. In this case you can apply a Darboux transformation with principal symbol p_x^2 using procedure LPDO__XXY_DT_with_Dxx`); fi; J1:=0; J2:=i2-(diff(r, y, x))/r+(diff(r, y))*(diff(r, x))/r^2; J3:=(2*i3^3-5*i3^2*i2+4*i3*i2^2-i2^3-(diff(i3, y))*(diff(i2, x))-2*(diff(i3, y))*(diff(i4, y))+2*(diff(i3, y))*i5-2*(diff(i3, y))*r+(diff(i2, y))*(diff(i2, x))+2*(diff(i2, y))*(diff(i4, y))-2*(diff(i2, y))*i5+2*(diff(i2, y))*r+(diff(i2, y, x))*i3-(diff(i2, y, x))*i2+2*(diff(i4, y, y))*i3-2*(diff(i4, y, y))*i2-2*(diff(i5, y))*i3+2*(diff(i5, y))*i2+2*(diff(r, y))*i3-2*(diff(r, y))*i2)/2/(i3-i2)^2; J4:=(8*(diff(i2, x))*r^3+8*i5^2*r^2-16*i5*r^3+3*(diff(r, x))^2*i3^2+3*(diff(r, x))^2*i2^2+8*(diff(i4, y))^2*r^2+16*(diff(i4, y))*r^3+2*(diff(i2, x))^2*r^2+8*r^4-8*i4*r^2*i3*i2+4*(diff(i4, y))*r*(diff(r, x))*i3-4*(diff(i4, y))*r*(diff(r, x))*i2+2*(diff(i2, x))*r*(diff(r, x))*i3-2*(diff(i2, x))*r*(diff(r, x))*i2-4*i5*r*(diff(r, x))*i3+4*i5*r*(diff(r, x))*i2+4*i4*r^2*i3^2+4*i4*r^2*i2^2+8*(diff(i4, y))*(diff(i2, x))*r^2-16*(diff(i4, y))*i5*r^2-8*(diff(i2, x))*i5*r^2+4*r^2*(diff(r, x))*i3-4*r^2*(diff(r, x))*i2-6*(diff(r, x))^2*i3*i2-2*(diff(r, x, x))*r*i3^2-2*(diff(r, x, x))*r*i2^2+4*(diff(r, x, x))*r*i3*i2)/(-4)/(i3-i2)^2/r^2; J5:=(-12*(diff(i4, y, y))*i5*r*i2+4*(diff(i4, y, y))*(diff(r, x))*i3*i2-6*(diff(i4, y, y))*(diff(i2, x))*r*i3+12*(diff(i5, y))*(diff(i4, y))*r*i3-12*(diff(i5, y))*(diff(i4, y))*r*i2-12*(diff(i5, y))*i5*r*i3+12*(diff(i5, y))*i5*r*i2-4*(diff(i5, y))*(diff(r, x))*i3*i2-6*(diff(i2, x))*r*(diff(r, y))*i3+6*(diff(i2, x))*r*(diff(r, y))*i2-12*(diff(i4, y))*r*(diff(r, y))*i3+12*(diff(i4, y))*r*(diff(r, y))*i2+12*i5*r*(diff(r, y))*i3-12*i5*r*(diff(r, y))*i2+4*(diff(r, y))*(diff(r, x))*i3*i2+11*(diff(i2, x))*r*i3^2*i2-10*(diff(i2, x))*r*i3*i2^2+22*(diff(i4, y))*r*i3^2*i2-20*(diff(i4, y))*r*i3*i2^2-12*(diff(i3, x))*r*i3^2*i2+2*(diff(i2, y, x))*(diff(r, x))*i3*i2+12*(diff(i3, x))*r*i3*i2^2-10*i5*r*i3^2*i2+8*i5*r*i3*i2^2+12*(diff(i3, y))*(diff(i2, x))*(diff(i4, y))*r-12*(diff(i3, y))*(diff(i2, x))*i5*r+(diff(i3, y))*(diff(i2, x))*(diff(r, x))*i3-(diff(i3, y))*(diff(i2, x))*(diff(r, x))*i2-24*(diff(i3, y))*(diff(i4, y))*i5*r+2*(diff(i3, y))*(diff(i4, y))*(diff(r, x))*i3-2*(diff(i3, y))*(diff(i4, y))*(diff(r, x))*i2-2*( diff(i3, y))*i5*(diff(r, x))*i3+2*(diff(i3, y))*i5*(diff(r, x))*i2+2*(diff(i3, y))*r*(diff(r, x))*i3-2*(diff(i3, y))*r*(diff(r, x))*i2-12*(diff(i2, y))*(diff(i2, x))*(diff(i4, y))*r+12*(diff(i2, y))*(diff(i2, x))*i5*r-(diff(i2, y))*(diff(i2, x))*(diff(r, x))*i3+(diff(i2, y))*(diff(i2, x))*(diff(r, x))*i2+24*(diff(i2, y))*(diff(i4, y))*i5*r-2*(diff(i2, y))*(diff(i4, y))*(diff(r, x))*i3+2*(diff(i2, y))*(diff(i4, y))*(diff(r, x))*i2+2*(diff(i2, y))*i5*(diff(r, x))*i3-2*(diff(i2, y))*i5*(diff(r, x))*i2-2*(diff(i2, y))*r*(diff(r, x))*i3+2*(diff(i2, y))*r*(diff(r, x))*i2+6*(diff(i5, y))*(diff(i2, x))*r*i3-6*(diff(i5, y))*(diff(i2, x))*r*i2-4*(diff(i2, x))*r*i3^3+3*(diff(i2, x))*r*i2^3-8*(diff(i4, y))*r*i3^3+6*(diff(i4, y))*r*i2^3+4*(diff(i3, x))*r*i3^3-4*(diff(i3, x))*r*i2^3+4*i5*r*i3^3-2*i5*r*i2^3-2*r^2*i3^2*i2+4*r^2*i3*i2^2+7*(diff(r, x))*i3^3*i2-9*(diff(r, x))*i3^2*i2^2+5*(diff(r, x))*i3*i2^3+3*(diff(i3, y))*(diff(i2, x))^2*r+12*(diff(i3, y))*(diff(i2, x))*r^2+12*(diff(i3, y))*(diff(i4, y))^2*r+24*(diff(i3, y))* (diff(i4, y))*r^2+12*(diff(i3, y))*i5^2*r-24*(diff(i3, y))*i5*r^2-3*(diff(i2, y))*(diff(i2, x))^2*r-12*(diff(i2, y))*(diff(i2, x))*r^2-12*(diff(i2, y))*(diff(i4, y))^2*r-24*(diff(i2, y))*(diff(i4, y))*r^2-12*(diff(i2, y))*i5^2*r+24*(diff(i2, y))*i5*r^2+12*(diff(i5, y))*r^2*i3-12*(diff(i5, y))*r^2*i2+2*(diff(i5, y))*(diff(r, x))*i3^2+2*(diff(i5, y))*(diff(r, x))*i2^2-12*r^2*(diff(r, y))*i3+12*r^2*(diff(r, y))*i2-2*(diff(r, y))*(diff(r, x))*i3^2-2*(diff(r, y))*(diff(r, x))*i2^2-6*(diff(i2, y, x))*i5*r*i2-2*r^2*i2^3-2*(diff(r, x))*i3^4-(diff(r, x))*i2^4+12*(diff(i3, y))*r^3-12*(diff(i2, y))*r^3+6*(diff(i2, y, x))*i5*r*i3+6*(diff(i2, y, x))*(diff(i4, y))*r*i2-6*(diff(i2, y, x))*(diff(i4, y))*r*i3-2*(diff(i4, y, y))*(diff(r, x))*i2^2-6*(diff(i2, y, x))*r^2*i3+12*(diff(i4, y, y))*r^2*i2-2*(diff(i4, y, y))*(diff(r, x))*i3^2-12*(diff(i4, y, y))*r^2*i3-(diff(i2, y, x))*(diff(r, x))*i2^2+6*(diff(i2, y, x))*r^2*i2-(diff(i2, y, x))*(diff(r, x))*i3^2+6*(diff(i4, y, y))*(diff(i2, x))*r*i2-12*(diff(i4, y, y))*(diff(i4, y))* r*i3+12*(diff(i4, y, y))*(diff(i4, y))*r*i2+12*(diff(i4, y, y))*i5*r*i3-3*(diff(i2, y, x))*(diff(i2, x))*r*i3+3*(diff(i2, y, x))*(diff(i2, x))*r*i2)/4/r/(i3-i2)^3; RETURN([J1,J2,J3,J4,J5]); end: ######################################################################### # # # performs *INVERTIBLE* DT for L with symbol XY(X+Y) with M=Dx+m on invariants # LPDO__XYS_DT_with_Dx:=proc(i1,i2,i3,i4,i5,op::integer:=0) local J1,J2,J3,J4,J5,x,y,d; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; x:=LPDO__vars[1]; y:=LPDO__vars[2]; if simplify(i4+2*i2-diff(i1,x))<>0 then print(`Invertible DT with dx for such operator exists only if i4+2*i2-diff(i1,x)=0, i.e.`,factor(i4+2*i2-diff(i1,x))=0); if op=0 then RETURN([]); fi: fi; d:=simplify(eval(2*(-i5+i1*i2+diff(i4,y)/2))); if d=0 then print(`The operator has a factorization, where the right factor has principal symbol p_x. It does not have DT with p_x of type I.`); RETURN([]); fi; J1:=i1-(diff(d, x))/d; J2:=-(diff(d, y))*(diff(d, x))/d^2+(diff(d, y, x))/d+i2; J3:=i2+(diff(d, y, x))/d+i3-(diff(d, y))*(diff(d, x))/d^2; J4:=0; J5:=(diff(d, y))*(diff(d, x))^2/d^3+diff(i3, x)-(1/2)*(diff(i4, y))-(diff(i2, y))+(1/2)*(diff(d, y, x, x))/d-(1/2)*d-(1/2)*(diff(d, y))*(diff(d, x, x))/d^2-(diff(d, x))*(diff(d, y, x))/d^2; #####print(`new tail is`,factor(J5-diff(J4,y)/2-J2*J1)); RETURN([J1,J2,J3,J4,J5]); end: ######################################################################### # # # performs *INVERTIBLE* DT for L with symbol XY(X+Y) with M=Dy+m on invariants # LPDO__XYS_DT_with_Dy:=proc(i1,i2,i3,i4,i5) local J1,J2,J3,J4,J5,x,y,d; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; x:=LPDO__vars[1]; y:=LPDO__vars[2]; if simplify(factor(diff(i1,y)+2*i2-i3))<>0 then print(`Invertible DT with dy for such operator exists only if diff(i1,y)+2*i2-i3=0, i.e.,`,factor(diff(i1,y)+2*i2-i3)=0); RETURN([]); fi; d:=simplify(eval(i5-(1/2)*(diff(i3, x)))); if d=0 then print(`The operator has a factorization, where the right factor has principal symbol p_y. It does not have DT with p_y of type I.`); RETURN([]); fi; J1:=-(diff(d, y))/d+i1; J2:=(diff(d, y))*(diff(d, x))/d^2+i2-(diff(d, y, x))/d; J3:=0; J4:=-i2+(diff(d, y, x))/d+i4-(diff(d, y))*(diff(d, x))/d^2; J5:=diff(i4, y)+diff(i2, x)-(1/2)*(diff(i3, x))+i2*i1+(1/2)*(diff(d, y, y, x))/d-(1/2)*(diff(d, x))*(diff(d, y, y))/d^2-(diff(d, y, x))*i1/d-(diff(d, y))*i2/d+(diff(d, y))*(diff(d, x))*i1/d^2+d; RETURN([J1,J2,J3,J4,J5]); end: ######################################################################### # # # performs *INVERTIBLE* DT for L with symbol XY(X+Y) with M=Dx+Dy+m on invariants # LPDO__XYS_DT_with_Dx_Dy:=proc(i1,i2,i3,i4,i5) local J1,J2,J3,J4,J5,x,y,f; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; x:=LPDO__vars[1]; y:=LPDO__vars[2]; if factor(i3-i4-i2)<>0 then print(`Invertible DT with dx+dy for such operator exists only if i3-i4-i2=0, i.e.`,simplify(i3-i4-i2)=0); RETURN([]); fi; f:=simplify(eval(i5+i1*i4+diff(i1,x,y)/2)); if f=0 then print(`The operator has a factorization, where the right factor has principal symbol p_x+p_y. It does not have DT with p_x+p_y of type I.`); RETURN([]); fi; J1:=i1+(diff(f, x))/f+(diff(f, y))/f; J2:=i3-i4; J3:=-(diff(f, y))^2/f^2+2*i3+diff(i1, y)-i4+(diff(f, y, y))/f-(diff(f, y))*(diff(f, x))/f^2+(diff(f, y, x))/f; J4:=-i3+diff(i1, x)+2*i4+(diff(f, x, x))/f-(diff(f, x))^2/f^2-(diff(f, y))*(diff(f, x))/f^2+(diff(f, y, x))/f; J5:=-(1/2)*(diff(f, y))*(diff(f, x, x))/f^2-(diff(f, y))*(diff(f, y, x))/f^2-(diff(f, x))*(diff(f, y, x))/f^2-(1/2)*(diff(f, x))*(diff(f, y, y))/f^2+f+(1/2)*(diff(f, y, y, x))/f+(diff(f, y))*(diff(f, x))^2/f^3+(diff(f, y))^2*(diff(f, x))/f^3+(1/2)*(diff(f, y, x, x))/f-(diff(f, y))*i4/f+diff(i3, x)+diff(i4, y)-i1*i4-(diff(f, x))*i4/f+(1/2)*(diff(i1, y, x)); RETURN([J1,J2,J3,J4,J5]); end: ##################################################### ######################################################################### ########################################################################## # # performs a step of Equivalence-Shift procedure for L of arbitrarily form # LPDO__a_step_of_Eq_and_Shift:=proc(L,M,N,L1,Minv,Ninv,A,G,S,T) local Ln, Mn, Nn, L1n, Minvn, Ninvn, An,Gn; Mn:=LPDO__add(M,LPDO__mult(S,L)): Nn:=LPDO__add(N,LPDO__mult(L1,S)): Ln:=LPDO__add(L,LPDO__mult(T,Mn)): L1n:=LPDO__add(L1,LPDO__mult(Nn,T)): An:=LPDO__minus(A,LPDO__mult(Minv,S)): Gn:=LPDO__minus(G,LPDO__mult(S,Ninv)): Minvn:=LPDO__minus(Minv,LPDO__mult(An,T)): Ninvn:=LPDO__minus(Ninv,LPDO__mult(T,Gn)): RETURN([Ln,Mn,Nn,L1n,Minvn,Ninvn,An,Gn]); end: ########################################################################## # # performs a step of Equivalence-Shift procedure for L of arbitrarily form # Does not compute the auxilary operators, just the change in L and M, L1, N # # M -> M+S*L # N -> N+L1*S # LPDO__a_step_of_Eq_and_Shift_basic:=proc(L,M,N,L1,S,T) local Ln, Mn, Nn, L1n, Minvn, Ninvn, An,Gn; Mn:=LPDO__add(M,LPDO__mult(S,L)): Nn:=LPDO__add(N,LPDO__mult(L1,S)): Ln:=LPDO__add(L,LPDO__mult(T,Mn)): L1n:=LPDO__add(L1,LPDO__mult(Nn,T)): RETURN([Ln,Mn,Nn,L1n]); end: ########################################################################## ########################################################################## # # returns full data for invertible DT for L=CM+f with symbol XXY # for type I # # C,M -- operators, M is the same one as in NL=L1M # f -- a function LPDO__data_for_invertible_DT_of_Type_I:=proc(C,M,f) local F,L,N,L1,Minv,Ninv,A,G,tmp; if f=0 then ERROR(`The third argument, f cannot be zero, since then an operator L is divisible by M, which means that the transformation cannot be invertible`); fi; F:=Create_LPDO_from_coeff([f]): L:=LPDO__add(LPDO__mult(C,M),F): #LPDO__print(L); N:=LPDO__conj(M,1/f): L1:=LPDO__add(LPDO__mult(N,C),F): A:=Create_LPDO_from_coeff([1/f]): G:=Create_LPDO_from_coeff([1/f]): Minv:=LPDO__mult1(-1/f,C): tmp:=Create_LPDO_from_coeff([-1/f]): Ninv:=LPDO__mult(C,tmp): RETURN([L,M,N,L1,Minv,Ninv,A,G]); end: ########################################################################## # # performs a step of Equivalence-Shift procedure for L with symbol XXY # for type I so for L=CM+f # # C,M -- operators, M is the same one as in NL=L1M # f -- a function LPDO__a_step_of_Eq_and_Shift_for_Type_I:=proc(C,M,f,S,T) local res,res2; res:=LPDO__data_for_invertible_DT_of_Type_I(C,M,f,S,T); res2:=LPDO__a_step_of_Eq_and_Shift(op(res),S,T): RETURN(res2); end: ########################################################################## # # performs a step of invertible Darboux transformation of type I # RETURNs [N,L1,Minvn,Ninvn,A,G], where # LPDO__DT_of_type_I:=proc(C,M,f) local N, L1, F, T, Minv, Ninv,A,G; N:=LPDO__conj(M,1/f); F:=Create_LPDO_from_coeff([f]); L1:=LPDO__add(LPDO__mult(N,C),F); A:=Create_LPDO_from_coeff([1/f]); G:=Create_LPDO_from_coeff([1/f]); ### G:=Create_LPDO_from_coeff([1/f]); T:=Create_LPDO_from_coeff([-1/f]); Minv:=LPDO__mult(T,C); Ninv:=LPDO__mult(C,T); RETURN([N,L1,Minv,Ninv,A,G]); end: ############################################################################################# ############################################################################################# # # the chain can be input as an array of numbers 1,2, and 3, e.g. [1,2,1,3,3,1], # which says that first one should apply a DT with X, then with Y, then with X, then with X+Y symbols. # LPDO__XYS_chain_of_first_order_invertible_DTs:=proc(chain,i1,i2,i3,i4,i5) local i, currentInvs,Invs; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; Invs:=[[i1,i2,i3,i4,i5]]; currentInvs:=[i1,i2,i3,i4,i5]; ###print(currentInvs); for i from 1 to nops(chain) do if nops(currentInvs)=0 then RETURN(currentInvs); fi: if chain[i]=1 then currentInvs:=LPDO__XYS_DT_with_Dx(op(currentInvs)); elif chain[i]=2 then currentInvs:=LPDO__XYS_DT_with_Dy(op(currentInvs)); elif chain[i]=3 then currentInvs:=LPDO__XYS_DT_with_Dx_Dy(op(currentInvs)); else ERROR(`the first argument should be a list of numbers 1,2,3 in some order`); fi: currentInvs:=map(el->factor(simplify(el)),currentInvs); #print(currentInvs); Invs:=[op(Invs),currentInvs]; if nops(currentInvs)=0 then RETURN(Invs); fi; od: RETURN(Invs); end: ################################################################################################ ################################################################################################ ################################################################################################ # # a -- the point of the chart from which you can test existence of Darboux transformations with p_x, p_y and p_x+p_y # p1,p2,p3 are the three new points # q1,q2,q3 are the values of the five invariants at these points # # LPDO__XYS_Chart_Plot_step_of_chain_of_first_order_invertible_DTs:=proc(a,i1,i2,i3,i4,i5) local p1,p2,p3,q1,q2,q3,p,q; if LPDO__nvars() <> 2 then ERROR(`illegal number of vars`); fi; p1:=a: p2:=a: p3:=a: q1:=[]: q2:=[]: q3:=[]: p:=[]: q:=[]: if simplify(eval(i4+2*i2-diff(i1,x)))=0 then p1:=[a[1]+1, a[2]]: q1:=LPDO__XYS_DT_with_Dx(i1,i2,i3,i4,i5); if nops(q1)<>0 and q1<>[i1,i2,i3,i4,i5] then p:=[p1,op(p)]; q:=[q1,op(q)]; fi: #else print(`Invertible DT with dx for such operator exists only if i4+2*i2-diff(i1,x)=0, i.e.`,factor(i4+2*i2-diff(i1,x))=0); fi: if simplify(eval(diff(i1,y)+2*i2-i3))=0 then p2:=[a[1], a[2]+1]: q2:=LPDO__XYS_DT_with_Dy(i1,i2,i3,i4,i5); if nops(q2)<>0 and q2<>[i1,i2,i3,i4,i5] then p:=[p2,op(p)]; q:=[q2,op(q)]; fi: #else print(`Invertible DT with dy for such operator exists only if diff(i1,y)+2*i2-i3=0, i.e.,`,factor(diff(i1,y)+2*i2-i3)=0); fi: if simplify(eval(i3-i4-i2))=0 then p3:=[a[1]-1, a[2]-1]: q3:=LPDO__XYS_DT_with_Dx_Dy(i1,i2,i3,i4,i5); if nops(q3)<>0 and q3<>[i1,i2,i3,i4,i5] then p:=[p3,op(p)]; q:=[q3,op(q)]; fi: #else print(`Invertible DT with dx+dy for such operator exists only if i3-i4-i2=0, i.e.`,simplify(i3-i4-i2)=0); fi: RETURN([p,q]): end: ################################################################################################ # # computes (including plotting) the next level of DTs chart # # output: # pl=[vector, vector, vector] # A:=[A1,A2,A3] # A1=[[p1,p2,p3],[q1,q2,q3]] # p1,p2,p3 are the three new points # q1,q2,q3 are the values of the five invariants at these points # LPDO__Charts_next_level:=proc(res) local i,A, el,pl,pts,invs,col,list_of_colors,j; pl:=[]; pts:=[]; invs:=[]; list_of_colors:=[[127, 255, 212],[255, 69, 0],[255, 255, 0],[154, 205, 50],[255, 0, 255],[138, 43, 226],[255, 215, 0]]: j:=rand(1..nops(list_of_colors))(); for i from 1 to nops(res[1]) do A:=LPDO__XYS_Chart_Plot_step_of_chain_of_first_order_invertible_DTs(res[1][i],op(res[2][i])); for el in A[1] do pl:=[plottools:-arrow(res[1][i],el,0.00000000001,0.08,0.06,color=ColorTools:-Color( list_of_colors[j])),op(pl)]: pts:=[op(A[1]),op(pts)]: od: for el in A[2] do invs:=[op(A[2]),op(invs)]: od: od: RETURN([pl,pts,invs]); end: ############################################################################################### # # LPDO__clean_chart:=proc(l) local i,pts,invs; if nops(l[2])=0 then RETURN(l); fi: pts:=[l[2][1]]; invs:=[l[3][1]]; for i from 1 to nops(l[2]) do if member(l[2][i],pts)=false then pts:=[op(pts),l[2][i]]; invs:=[op(invs),l[3][i]]; fi: od: RETURN([l[1],pts,invs]); end: ############################################################################################### # # computes (including plotting) b levels of DTs chart # Unfortunately it does not trace cycles. # LPDO__XYZ_Chart:=proc(b,list_of_invs) local el,l,i,full_list; l:=LPDO__Charts_next_level([[[0,0]],[list_of_invs]]): if nops(l[1])=0 then print(`First empty level is 1`); RETURN([[],[],[]]); fi: l:=LPDO__clean_chart(l); l[3]:=map(el->simplify(el),l[3]): full_list:=l; for i from 1 to b-1 do l:=LPDO__Charts_next_level([l[2],l[3]],op): l:=LPDO__clean_chart(l); l[3]:=map(el->simplify(el),l[3]): if nops(l[1])=0 then print(`First empty level is`,i+1); RETURN(full_list); fi: full_list[1]:=[op(l[1]),op(full_list[1])]: full_list[2]:=[op(l[2]),op(full_list[2])]: full_list[3]:=[op(l[3]),op(full_list[3])]: od: RETURN(full_list); end: ################################################################################################