######################################################################## GeneralizedLaplace := table(): mypackagename:="GeneralizedLaplace"; ######################################################################## # # (C) Pieter Eendebak, eendebak@math.uu.nl # # See http://www.math.uu.nl/people/eendebak # ######################################################################## print(""||mypackagename||`, by Pieter Eendebak`); print(`This packages contains software to calculate generalized Laplace sequences.`); ######################################################################## # # init_equation initializes all structures # ######################################################################## GeneralizedLaplace[init_equation]:=proc(eq) local srho, ssigma, stau, rho, sigma, tau; global ivar,dvar,varrho, Fr, Ft, Fs, Fz, Fp, Fq, eqsimp, eqsimp_rho, eqsimp_sigma, eqsimp_tau, char_eq, mx, my, nx, ny, kappa, delta, sol, P, Q, X, Y; # initialize variables ivar:=[x,y]; dvar:=[u]; # for simplicity we take varrho (introduced on page 26 of Juras' thesis) equal to 1 varrho:=1; # figure out for which of the second order variables we can solve rho:=solve(eq,u[x,x]); srho:=nops([rho]); if(srho=0) then print("Could not solve for rho") else eqsimp_rho:=expr->simplify(jsubs(u[x,x]=[rho][1],expr, ivar, dvar), symbolic): print("Simplification method: rho = ",rho); fi; sigma:=solve(eq,u[x,y]); ssigma:=nops([sigma]); if(ssigma=0) then print("Could not solve for sigma") else eqsimp_sigma:=expr->simplify(jsubs(u[x,y]=[sigma][1],expr, ivar, dvar), symbolic): print("Simplification method: sigma = ", sigma); fi; tau:=solve(eq,u[y,y]); stau:=nops([tau]); if(stau=0) then print("Could not solve for tau") else eqsimp_tau:=expr->simplify(jsubs(u[x,y]=[tau][1],expr, ivar, dvar), symbolic): print("Simplification method: tau = ", tau); fi; # report results if(srho>1 or ssigma>1 or stau>1) then WARNING("Multiple solve methods possible") fi; if(srho+ssigma+stau=0) then error("Could not solve for any highest order variable") fi; if(srho>0) then eqsimp:=eqsimp_rho elif(ssigma>0) then eqsimp:=eqsimp_sigma; else eqsimp:=eqsimp_tau; fi; Fr:=partder(eq,u[x,x],ivar,dvar): Fs:=partder(eq,u[x,y],ivar,dvar): Ft:=partder(eq,u[y,y],ivar,dvar): Fz:=partder(eq,u,ivar,dvar): Fp:=partder(eq,u[x],ivar,dvar): Fq:=partder(eq,u[y],ivar,dvar): # the characteristic equation char_eq:= Fr*v^2 - Fs*v*w + Ft*w^2; sol:=solveHomQuadraticKappa(char_eq,v,w); mx:=eqsimp(sol[1][2][2]); my:=eqsimp(sol[1][2][1]); nx:=eqsimp(sol[1][1][2]); ny:=eqsimp(sol[1][1][1]); kappa:=eqsimp(sol[2]); delta:=mx*ny-my*nx; # the characteristic total vector fields, defined as operators here X:=jet->mx*totalder(jet,x,ivar,dvar)+my*totalder(jet,y,ivar,dvar): Y:=jet->nx*totalder(jet,x,ivar,dvar)+ny*totalder(jet,y,ivar,dvar): # formula 4.4 and 4.5 P:=delta^(-1)*( ny*X(nx)-nx*X(ny)+nx*Y(my)-ny*Y(mx) ); Q:=delta^(-1)*( mx*X(ny)-my*X(nx)+my*Y(mx)-mx*Y(my) ); print("Equation initialized"); end proc: ######################################################################## # # laplace_invariants # # The following function calculated the first order Laplace invariants # and associated functions from the global data generated by init_equation. # We define the generalized coefficients A, B, C and the invariants H0, K0 # using the formulas from Anderson and Juras . # ######################################################################## GeneralizedLaplace[laplace_invariants]:=proc() global A, B, C, AR, BR, CR, H0, K0, f0, b0; local A0, B0, C0; # formulas 4.12 in Juras thesis A0:=eqsimp(delta^(-1)*((kappa*Fp-X(nx))*ny-(kappa*Fq-X(ny))*nx)); B0:=eqsimp(delta^(-1)*(-(kappa*Fp-X(nx))*my+(kappa*Fq-X(ny))*mx)); C0:=kappa*Fz; # formulas 4.10a,b A:=A0-Y(varrho)/varrho; B:=B0-X(varrho)/varrho; C:=C0-X(Y(varrho))/varrho - A0*X(varrho)/varrho - B0*Y(varrho)/varrho+2*X(varrho)*Y(varrho)/varrho^2; # from formula 4.11, the conjugate functions AR:=eqsimp(A+P); BR:=eqsimp(B+Q); CR:=C; # equation 4.14 and 4.23 H0:=eqsimp(X(A)+A*B-C); K0:=eqsimp(Y(BR)+BR*AR-C); # the elements f0 and g0 contain all the data needed to calculate higher order invariants f0:=[A,B,C,H0]; b0:=[AR,BR,CR,K0]; [H0,K0]; #print("Invariants generated"); end proc: ######################################################################## # # forwardLaplace, backwardLaplace # # From the equation we can calculate the functions A, B and C. These functions # completely determined the Laplace invariants and the transformed equations. # Then inductively we can calculate all higher order invariants. # ######################################################################## GeneralizedLaplace[forwardLaplace]:=proc(data) local A, B, C, H, fA,fB,fC,fH,fK; global P,Q,X,Y; A:=data[1]; B:=data[2]; C:=data[3]; H:=data[4]; if(H=0) then #error "Transformation not possible for zero invariant" return [0,0,0,0] fi; # the prolongation formulas 4.31 fA:=A-Y(log(H))-P; fB:=B-Q; fC:=C-X(A)-B*Y(log(H))+Y(B); fH:=X(fA)+fA*fB-fC; [fA,fB,fC,fH]; end proc: GeneralizedLaplace[backwardLaplace]:=proc(data) local AR, BR, CR, K, bA,bB,bC,bK; global P,Q,X,Y; AR:=data[1]; BR:=data[2]; CR:=data[3]; K:=data[4]; if(K=0) then #error "Transformation not possible for zero invariant" return [0,0,0,0] fi; # formula 4.32 bA:=AR+P; bB:=BR-X(log(K))+Q; bC:=CR-Y(BR)-AR*X(log(K))+X(AR); bK:=Y(bB)+bA*bB-bC; [bA,bB,bC,bK]; end proc: ######################################################################## # # monge_invariants # # We calculate the Monge invariants using explicit equations. # ######################################################################## GeneralizedLaplace[monge_invariants]:=proc() local mx_r, my_r, mx_s, my_s, mx_t, my_t, nx_r, ny_r, nx_s, ny_s, nx_t, ny_t, t1, t2, t3, s1, s2, s3; global Ms, Mt; mx_r:=partder(mx, u[x,x], ivar, dvar): my_r:=partder(my, u[x,x], ivar, dvar): mx_s:=partder(mx, u[x,y], ivar, dvar): my_s:=partder(my, u[x,y], ivar, dvar): mx_t:=partder(mx, u[y,y], ivar, dvar): my_t:=partder(my, u[y,y], ivar, dvar): nx_r:=partder(nx, u[x,x], ivar, dvar): ny_r:=partder(ny, u[x,x], ivar, dvar): nx_s:=partder(nx, u[x,y], ivar, dvar): ny_s:=partder(ny, u[x,y], ivar, dvar): nx_t:=partder(nx, u[y,y], ivar, dvar): ny_t:=partder(ny, u[y,y], ivar, dvar): t1:=(ny_r*nx-nx_r*ny)*ny^2; t2:=(ny_s*nx-nx_s*ny)*nx*ny; t3:=(ny_t*nx-nx_t*ny)*nx^2; s1:=(my_r*mx-mx_r*my)*my^2; s2:=(my_s*mx-mx_s*my)*mx*my; s3:=(my_t*mx-mx_t*my)*mx^2; # from the formula's for M_sigma and M_tau on page 71 Mt:=(-1/(delta^3*varrho))*(t1-t2+t3): Ms:=(-1/(delta^3*varrho))*(s1-s2+s3): [Ms,Mt]; end proc: ######################################################################## # old #save (GeneralizedLaplace, `GeneralizedLaplace.m`): # set variables mypackagenamelib:=cat(mypackagename, ".lib"): cdir:=currentdir(): mypackagenameFull:=""||cdir||"/"||mypackagenamelib; # remove old lib system("rm "||mypackagenamelib); # create new lib print("Creating archive"); march('create', mypackagenamelib, 100); savelibname:=mypackagenameFull; savelib(convert(mypackagename,symbol)); marchlist:=march('list', mypackagenameFull); map(march, ['gc','reindex','pack'], mypackagenameFull); march('setattribute', mypackagenameFull, mode="READONLY"); print(marchlist); quit; ########################################################################