restart;
with(DifferentialGeometry):
with(DifferentialGeometry:-Tools):
with(JetCalculus):
Preferences("JetNotation", "JetNotation2"):
pEqs := eval([
G*r(t, a)*diff(h(a), a) + r(t, a)*diff(u(t, a), a)*u(t, a) + r(t, a)*diff(u(t, a), t) + diff(p(t, a), a),
diff(r(t, a), a)*u(t, a) + r(t, a)*diff(u(t, a), a) + diff(r(t, a), t),
-k*diff(T(t, a), a, a) + r(t, a)*T(t, a)*(diff(s(t, a), t) + u(t, a)*diff(s(t, a), a))
],
{
h(a) = lambda*a^1,
p(t,a) = P(r(t, a), T(t, a)),
s(t,a) = S(r(t, a), T(t, a))
});
DGsetup([t, a], [u, r, T], C, 5);
Eqs := convert(pEqs, DGjet);
Eqs2 := [TotalDiff(Eqs[1], t), TotalDiff(Eqs[1], a), TotalDiff(Eqs[2], t), TotalDiff(Eqs[2], a)]
EulerRestr := solve( {op(Eqs), op(Eqs2)}, {T[0,2], r[1,1], r[2,0], u[0,2], u[1,1], u[1,0], r[1,0]} ):
Invs := [
# x y
r[0,0], T[0,0],
u[0,1], # K
r[0,1], # L
T[0,1], # M
(T[1,0] + u[0,0]*T[0,1]) # N
];
onInvs := F->local k; eval(F, solve( {seq( Invs[k] = J||k, k = 1..nops(Invs) )}, {r[0,0], T[0,0], u[0,1], r[0,1], T[0,1], T[1,0]} ));
Tresse := (F, J1, J2) -> eval([A,B], solve(
{
A*TotalDiff(J1, [1,0]) + B*TotalDiff(J2, [1,0]) - TotalDiff(F, [1,0]),
A*TotalDiff(J1, [0,1]) + B*TotalDiff(J2, [0,1]) - TotalDiff(F, [0,1])
},
{A, B}
)):
E2 := simplify(onInvs(eval(
[
op(Tresse(Invs[3], Invs[1], Invs[2])),
op(Tresse(Invs[4], Invs[1], Invs[2])),
op(Tresse(Invs[5], Invs[1], Invs[2])),
op(Tresse(Invs[6], Invs[1], Invs[2]))
],
EulerRestr))):
nops(E2);
ex1:=eliminate({E2[1]-D31, E2[2]-D32, E2[3]-D41, E2[4]-D42, E2[5]-D51, E2[6]-D52, E2[7]-D61, E2[8]-D62}, {u[2,0], r[0,2], T[1,1], T[2,0]});
PDETools:-declare((P, S, K, L, M, N)(x, y));
SysQ0 := subs({
J1=x, J2=y,
J3=K(x,y),
J4=L(x,y),
J5=M(x,y),
J6=N(x,y),
D31=diff(K(x,y),x), D32=diff(K(x,y),y),
D41=diff(L(x,y),x), D42=diff(L(x,y),y),
D51=diff(M(x,y),x), D52=diff(M(x,y),y),
D61=diff(N(x,y),x), D62=diff(N(x,y),y)
}, convert(ex1[2],diff));
nops(SysQ0);
S1 := simplify(SysQ0[1],size);
S2 := simplify(SysQ0[2] - S1*N(x,y))/k/M(x,y);
AA := simplify(SysQ0[3] + S1*x*diff(P(x,y),y)*L(x,y))/k:
BB := simplify(SysQ0[4] - S1*x*M(x, y)*diff(P(x, y), y))/k:
S3 := simplify(eval(BB, solve({AA}, {diff(M(x, y), y)}))*L(x, y)/x/(x*K(x, y)*M(x, y) + L(x, y)*N(x, y)));
S4 := simplify(AA/x^2, size);
onPlanck := F -> eval(F,
{
P(x,y)=-R*x^2*y*diff(phi(x,y),x),
S(x,y)=R*(phi(x,y)+y*diff(phi(x,y),y))
});
Z := simplify( [seq(onPlanck(S||k), k=1..4)] );
virial := k ->local i; phi(x,y) = n/2*ln(y) - ln(x) - add(x^i/i*A||i(y), i=1..k);
ser := (Q, n) -> local k;
Q(x,y) = expand(
x^((Q||d)) * add( Q||k(y)*x^k,k=0..n)
);
ser(M, 2);
degsol := {Kd = 1, Ld = 2, Md = 1, Nd = 1};
exs := n->eval( {virial(n), ser(K, n), ser(L, n), ser(M, n), ser(N, n)}, degsol );
W1 := collect(expand(eval(Z[1], exs(0))), x) / x^2;
W2 := collect(expand(eval(Z[2], exs(0))), x) / x^2;
W3 := collect(expand(eval(Z[3], exs(0))), x) / x^3;
W4 := collect(expand(eval(Z[4], exs(0))), x) / x^4;
simplify(eval( simplify((collect(expand(eval(Z[1], exs(1)))/x^2, x)-W1)/ x), x=0));
simplify(eval( simplify((collect(expand(eval(Z[2], exs(1)))/x^2, x)-W2)/ x), x=0));
simplify(eval( simplify((collect(expand(eval(Z[3], exs(1)))/x^3, x)-W3)/ x), x=0));
V4:=simplify(eval( simplify((collect(expand(eval(Z[4], exs(1)))/x^4, x)-W4)/ x), x=0),size);
collect(V4, {diff(K1(y), y), diff(L1(y), y), diff(M1(y), y), diff(N1(y), y)});
Ws := simplify( { W||(1..4) } );
solve(
Ws,
{
diff(K0(y), y),
diff(L0(y), y),
diff(M0(y), y),
diff(N0(y), y)
});
#####################################
idealG := {
P(x,y) = R*x*y,
S(x,y) = R*ln(y^(n/2)*x)
};
vdW := {
P(x,y) = 8*y/(3/x-1)-3*x^2,
S(x,y) = ln(y^(4*n/3)*(3/x-1)^(8/3))
};
pdsolve(onPlanck(vdW));
DGsetup([x,y], [K,L,M,N], CJ, 2);
SJ := convert({S||(1..4)}, DGjet);
onSJ := F->simplify(eval(F, solve({op(SJ)},{K[0,1], M[1,0], L[0,1], N[1,0]})));
Sym0 := DGzip(
eval([
A1(x),
A2(y),
-(diff(A2(y), y)*x - x*diff(A1(x), x) + A1(x))*K[0, 0]/x + 0*L[0, 0],
L[0,0]*(A4(x)-diff(A2(y),y)),
M[0,0]*((A4(x)-diff(A2(y),y)) - diff(A1(x), x) + diff(A2(y), y)),
0
], {A2(y)=a21*y+a22, A4(x)= a11 + a41, A1(x)=a11*x}
), DGinfo(CJ, "FrameJetVectors")[1..6]);
Sym1 := Prolong(DGsimplify(eval(Sym0, {a11=1, a21=0, a22=0, a41=0})), 1);
Sym2 := Prolong(DGsimplify(eval(Sym0, {a11=0, a21=1, a22=0, a41=0})), 1);
Sym3 := Prolong(DGsimplify(eval(Sym0, {a11=0, a21=0, a22=1, a41=0})), 1);
Sym4 := Prolong(DGsimplify(eval(Sym0, {a11=0, a21=0, a22=0, a41=1})), 1);
Sym :=evalDG( alpha1*Sym1+alpha2*Sym2+alpha3*Sym3+1/2*(alpha1-alpha4)*Sym4 );
dissolve := p->{coeffs(collect(expand(p), DGinfo("FrameJetVariables")[7..14], 'distributed'), DGinfo("FrameJetVariables")[7..14])};
ssym:=foldl(`union`, op(map(dissolve@numer, onSJ( map( X->LieDerivative(Sym, X), SJ ) )))):
filterout:=S->foldl(`*`, 1, op(map(p->p[1], select(p->p[2]>0 and has(p, {A||(1..6)}), factors(S)[2]))));
ssym1:=map(filterout, ssym):
nops(ssym);
ssym[1];
eq1:=simplify( op(2,ssym) );
eq2:=simplify( op(3,ssym) );
eq3:=simplify( op(4,ssym) );
eq4:=simplify( op(5,ssym) );
psq1:=({coeffs(collect(eq1, [K[0,0],L[0,0],M[0,0],N[0,0]], 'distributed'), [K[0,0],L[0,0],M[0,0],N[0,0]])});
psq2:=({coeffs(collect(eq2, [K[0,0],L[0,0],M[0,0],N[0,0]], 'distributed'), [K[0,0],L[0,0],M[0,0],N[0,0]])});
psq3:=({coeffs(collect(eq3, [K[0,0],L[0,0],M[0,0],N[0,0]], 'distributed'), [K[0,0],L[0,0],M[0,0],N[0,0]])});
psq4:=({coeffs(collect(eq4, [K[0,0],L[0,0],M[0,0],N[0,0]], 'distributed'), [K[0,0],L[0,0],M[0,0],N[0,0]])});
nops(psq4);
pdsolve(eval( onPlanck(psq1 union psq2 union psq3 union psq4), {} ));
pdsolve(eval( onPlanck(psq1 union psq2 union psq3 union psq4), {alpha2=0} ));
pdsolve(eval( onPlanck(psq1 union psq2 union psq3 union psq4), {alpha1=0} ));
expand(pdsolve(eval( onPlanck(psq1 union psq2 union psq3 union psq4), {alpha4=0} )));
eval(Sym,{alpha1=1,alpha3=0});
combine(pdsolve(eval( onPlanck(psq1 union psq2 union psq3 union psq4), {alpha1=1,alpha3=0} )));
SvdW := simplify( map(numer, eval({S||(1..4)}, vdW)) );
pdsolve( simplify(eval({op(SvdW[1..4])}, M(x,y)=0)));
pdsolve( simplify(eval({op(SvdW[1..4])}, L(x,y)=0)));
# pdsolve( simplify(eval({op(SvdW[1..4])}, N(x,y)=0)));
pdsolve( simplify(eval({S||(1..4)}, K(x,y)=0)));
pdsolve( simplify(eval({S||(1..4)}, M(x,y)=0)));