Wolfram Technology Conference 2017

A package for symbolic PDE analysis

Davit Shahnazaryan, KAUST
Joint work with Diogo Gomes and Mher Safaryan

Outline

Motivation

Introduction (conservation/dissipation laws)

Previous work

Key functions

Variational derivative

Monomial generation

Find group of discrete symmetries and generate polynomials

Integration by parts

Cleaning Null-Lagrangians

Simplifying integral expressions

Undetermined Coefficients Matrix

Finding conditions for positivity

Finding conservation laws

Finding dissipation laws

Motivation

In physics conservation/dissipation laws describe the conservation/dissipation of fundamental
physical quantities such as mass, energy, linear momentum and entropy.

Computations are often tedious and not obvious, when you do it by hand.

For nonlinear PDEs, the existence of a sufficiently large (in principal infinite) number of conservation laws assures complete integrability.

The problem

The initial problem that we started to solve is to find conservation (dissipation) laws for a given nonlinear evolution PDE or system of PDEs. For simplicity, consider a single evolution equation of the following form:

(Subscript[u, t])(x,t)=F(x,u(x,t), (Subscript[u, x])(x,t),(Subscript[u, xx])(x,t), (Subscript[u, xxx])(x,t),. . . ), x∈Ω

where x = (x,y,z,...)is high dimensional variable in the periodical domain   Ω (∂Ω=∅),

and (Subscript[u, x])(x,t) = ((Subscript[u, x])(x,t), (Subscript[u, y])(x,t), (Subscript[u, z])(x,t), . . .) represents the set of all partial derivatives with respect to spatial variables.

We are looking for expressions/functionals of the form

E(u)=E(x,u(x,t), (Subscript[u, x])(x,t),(Subscript[u, xx])(x,t), (Subscript[u, xxx])(x,t),. . . )

for which

(d/(dt))(∫_Ω)E(u)dx=0 - conservation law

(∫_Ω)E(u)dx≥constant, (d/(dt))(∫_Ω)E(u)dx≤0 - dissipation law

Simple example: heat equation

Consider the heat equation in one spatial dimension:

(Subscript[u, t])(x,t)=(a^2)(Subscript[u, xx])(x,t)

Conservation of energy

(d/(dt))(∫_Ω)u(x,t)dx =(a^2)(∫_Ω)(((∂^2)u(x,t))/(∂(x^2)))dx = (a^2)(∫_(∂Ω))((∂u(x,t))/(∂x))dx =0

Dissipation of entropy

(d/(dt))(∫_Ω)u(x,t)log u(x,t)dx=(Underscript[-∫, Ω])(((((∂u(x,t))/(∂x))))^2)((dx)/(u(x,t)))<0, if u(x,t)>0

Previous work

L.D. Poole and W. Hereman, ConservationLawsMD.m  : A Mathematica package for the symbolic computation of conservation laws of systems of nonlinear partial differential equations in multiple space dimensions (2009, Mathematica 7).

Ü. Göktas and W. Hereman, condens2009.m: A Mathematica program for the symbolic computation of conserved densities for systems of nonlinear evolution equations (1996, updated: 2009).

P.J. Adams and W. Hereman, TransPDEDensityFlux.m: A Mathematica program for the symbolic computation of conserved densities and fluxes for systems of partial differential equations with transcendental nonlinearities (2002, Mathematica 5).

Dimitri Tsoubeli and Stylianos Dimas, Sym.m: Partial differential equations, algebraic computing and nonlinear systems(2008, Mathematica 6)

Tasks

Generate Monomials/Polynomials

Integrate by parts

Clean Null-Lagrangians

Simplifying integral identities

Applying Euler-Lagrange operator

Find Undetermined Coefficients Matrix

And solve in respect to coefficients to find conserved and dissipated quantities

Euler-Lagrange operator or Variational Derivative

If the expression E(u) has the form

E(u)=E(x,u(x,t), (Subscript[u, x])(x,t),(Subscript[u, xx])(x,t), (Subscript[u, xxx])(x,t),. . . )
E(u)≡0, if u(x,t)≡0

Then

(∫_Ω)E(u)dx ≡0

is equivalent to

(Subscript[ℒ, u])E≡0 (null-lagrangian)

where variational derivative VarD or Euler-Lagrange operator is defined as

(Subscript[ℒ, u])= (∂/(∂u))- (d/(dx))(∂/(∂(Subscript[u, x])))+((d^2)/(d(x^2)))(∂/(∂(Subscript[u, xx])))-((d^3)/(d(x^3)))(∂/(∂(Subscript[u, xxx])))+ . . . , for u=u(x)

(Subscript[ℒ, u])= (∂/(∂u))-( (∂/(∂x))(∂/(∂(Subscript[u, x])))+ (∂/(∂y))(∂/(∂(Subscript[u, y]))))+(((∂^2)/(∂(x^2)))(∂/(∂(Subscript[u, xx])))+((∂^2)/(∂x∂y))(∂/(∂(Subscript[u, xy])))+((∂^2)/(∂(y^2)))(∂/(∂(Subscript[u, yy]))))-. . .,for u=u(x,y)

Time Derivative

Using the Euler-Lagrange operator we can define time derivative for the expression E(u) as follows:

(d/(dt))(∫_Ω)E(u)dx=(Underscript[∫ , Ω])(d/(dt))E(u)dx=(Underscript[∫ , Ω])(((∂E)/(∂u))ut+((∂E)/(∂(Subscript[u, x])))uxt+((∂E)/(∂(Subscript[u, xx])))(Subscript[u, xxt])+. . .)dx=(Underscript[∫ , Ω])(Subscript[ℒ, u])E(u) ·(Subscript[u, t])dx=(Underscript[∫ , Ω])(Subscript[ℒ, u])E(u) ·F(u)dx

Then using the above statement we can state that

E(u) is conseved if  (Subscript[ℒ, u])[(Subscript[ℒ, u])E(u) ·F(u)]≡0

E(u) is dissipated if  (Underscript[∫ , Ω])E(u)dx≥constant && (Underscript[∫ , Ω])(Subscript[ℒ, u])E(u) ·F(u)dx≤0

Variational Derivative

VarD[(1/2)((u'[x])^2), u, x]

VarD[(Sqrt[1+((Derivative[1,0][u][x,y])^2)]),u,{x,y}]//Simplify

VarD[(((m[x])^2)/2)+(((u'[x])^2)/(2 m[x])), {u,m}, {x}]

Monomial generation

SetOptions[Manipulator,Appearance->"Open"];
Manipulate[TraditionalForm@Monomials[vars, args, der, deg, Type->type],
Style["Monomial Generation",14,Italic],Delimiter,
{{vars,{u},"variables"},{u,v},ControlType->TogglerBar},
{{args,{x},"arguments"},{x,y,z},ControlType->TogglerBar},
{{deg,1,Style["maxDegree"]},1,5,1},
{{der,0,Style["maxDerivative"]},0,5,1},
{type,{"All", "Even","Odd"}},ContentSize -> {400,400},ControlPlacement->Left, LabelStyle->{14}]

Monomial generation with generators

Monomials[{u},{x},2,2, Generators->{(E^(u[x]))}]

Monomials[{u,v},{x},2,2, Generators->{(E^(u[x])), (1/(v[x]))}]

Monomials[{u,v},{x},2,2, Generators->{(E^(u[x])),((v[x])/(u[x]))}]

Symmetry group and polynomial generation

symmetryGroup=FindSymmetryGroup[{(e^(u[x,y])) ,((D[v[x,y],{x,3},{y,2}])^3)}, {u, v}, {x, y}];

symmetryGroup//TableForm

Manipulate[value,Style["Symmetric Polynomials",14,Italic],Delimiter,
{{deg,1,Style["maxDegree"]},1,5,1},
{{der,0,Style["maxDerivative"]},0,5,1},Style["Options",12,Italic],
{{value,{},""},
Row[{Button["Evaluate",value=TraditionalForm@StrictlySymmetricPols[{u,v}, {x,y}, der, deg, symmetryGroup[[All,2]], symmetryGroup[[All,3]]] ]}]&},
ContentSize -> {600,400},ControlPlacement->Left, LabelStyle->14]

Integration By Parts

Since we are considering periodical domain Ω (∂Ω=∅), then integration by parts will look like

(∫_Ω)(Subscript[u, x])(x,t)v(x,t)dx =-(∫_Ω)u(x,t)(Subscript[v, x])(x,t)dx

Examples

IntegrateByParts[x][a u'[x]^2 + b u[x]u''[x]]//Simplify

IntegrateByParts[x,y][D[u[x,y],x,y,x,x] D[v[x,y],y,x,y,y]]

IntegrateByParts[x,y][p E^u[x]u''[x]-q E^u[x]u'[x]^2]//Simplify

Enhanced Null Lagrangian Cleaner

ENLCleaner[{u''[t]v[t],u'[t]v'[t],u[t]v''[t]},{u,v},{t}]

ENLCleaner[{u'[x],v'[x],u'[x]v[x]+u[x]v'[x],E^u[x]u'[x]},{u},{x}]

ENLCleaner[{E^u[x]u''[x],E^u[x]u'[x]^2},{u},{x}]

Cancel Null Lagrangians or Simplifying integral identities

CancelNL[((u[x])^2)+5((u'[x])^2)+7u[x]u''[x], {u}, {x}]

CancelNL[p (E^(u[x]))((u'[x])^2)+(E^(u[x]))u''[x],{u},{x},{p}]

CancelNL[D[Derivative[1,3][v][x,y] Derivative[3,1][u][x,y],y],{u},{x,y}]

Undetermined Coefficients Matrix

ClearAll[u,v,x,a,b,c,d,e,p,q]
coeffs = {a,b,c,d,e};
pars ={p,q};
generators ={q u[x], p E^u[x],2^p/v[x],Sqrt[v[x]]/q};
monomials=Monomials[{u,v},{x},2,2, Generators->generators];
expr = Total[Total/@RandomChoice[monomials[[;;9]],{Length@coeffs,3}].coeffs];
Print[expr//TraditionalForm];
MatrixForm/@UCM[expr, coeffs,pars]//TraditionalForm

Finding conservation laws

Key steps of algorithm

Generate appropriate initial expressions (monomials/polynomials) - (Subscript[e, 1]),(Subscript[e, 2]),. . . ,(Subscript[e, n])

Monomial, SymmetricPolynomials

ENLCleaner, IntegrateByParts

Construct the general form of density and compute time derivative - (d/(dt))((Subscript[C, 1])(Subscript[e, 1])+(Subscript[C, 2])(Subscript[e, 2]),. . . +(Subscript[C, m])(Subscript[e, m]))

VarD, CancelNL, IntegrateByParts

Determine the unknown constant coefficients (Subscript[C, 1]),(Subscript[C, 2]),. . . ,(Subscript[C, m])

UCM, Solve

Examples

KdV equation: (Subscript[u, t]) + q(Subscript[uu, x]) + (Subscript[u, xxx]) = 0

monomials=Monomials[{u},{x},2,4]
monomials//Length

symmetryGroup= FindSymmetryGroup[{q u[x]u'[x]+ u'''[x]},{u},{x}]

symMonomials = StrictlySymmetricPols[{u},{x},2,4,symmetryGroup[[All,2]],symmetryGroup[[All,3]]]

monomials=ENLCleaner[monomials,{u},{x},{}]
monomials//Length

coefficients=Table[(Subscript[a, k]),{k,1,Length@monomials}]
generalForm=coefficients.monomials

VarD[VarD[generalForm,{u},{x}].{-q u[x] u'[x]-Derivative[3][u][x]},{u},{x}]//Simplify

matrix =UCM[%,coefficients,{q},Factors->False];
matrix//MatrixForm

Quiet@Solve[matrix.coefficients==0,coefficients]

generalForm/.%

ConservedDensity[{-q u[x] u'[x]-Derivative[3][u][x]}, {u},{x}, Parameters->{q}, Deg->{4}, Der->{2}]

ConservedDensity[{-q u[x] u'[x]-Derivative[3][u][x]}, {u},{x}, Parameters->{q}, Deg->{4}, Der->{2},Testing->True]

Heat equation: (Subscript[u, t])=(a^2)Δu

ConservedDensity[{(a^2)u''[x]},{u},{x}, Der->{2}, Deg->{3}, Parameters->{a}]

ConservedDensity[{(a^2)u''[x]},{u},{x}, Generators->{u[x],u'[x],x}, Parameters->{a}, Deg->{3}]

ConservedDensity[{(a^2)(Derivative[2,0,0][u][x,y,z]+Derivative[0,2,0][u][x,y,z]+Derivative[0,0,2][u][x,y,z])},{u},{x,y,z}, 
Generators->{u[x,y,z]}, Deg->{2}, Parameters->{a}]

Mean Field Games

(Subscript[v, t])=(((Subscript[v, x])^2)/2)+Log m                           Hamilton-Jacobi equation
(Subscript[m, t])=(Subscript[(m v), x])≡v[x] m'[x]+m[x] v'[x]             Transport equation  (Fokker-Planck)

ConservedDensity[{(((D[v[x],x])^2)/2)+Log[m[x]], D[m[x] v[x], x]},{v, m},{x}, Deg->{2}, Der->{1}]

(Subscript[v, t])=(Subscript[((1/2)(m^(-1+c))(v^2)), x])
(Subscript[m, t])=(Subscript[(((m )^c)v), x])

ConservedDensity[{D[v[x]^2/(2 m[x]^(1-c)), x], D[ m[x]^(c) v[x], x]}, {v, m}, {x},
Generators->{v[x], m[x],  m[x]^(c)}, Parameters->{c}, Deg->{2}]

Wave equation: (Subscript[u,  t t])=(a^2)Δu

(Subscript[u, t])=v
(Subscript[v, t])=(a^2)Δu

ConservedDensity[{v[x,y], (a^2) (Derivative[2,0][u][x,y]+Derivative[0,2][u][x,y])}, {u,v}, {x,y}, Parameters->{a}, Der->{1}, Deg->{2}]

Finding dissipation laws: conditions for positivity

Short way

PositivityBasic[expr_, vars_, facts_:True]:=Resolve[ForAll[vars,facts,expr≥0], Reals];

PositivityBasic[a ((u[x])^3) + b u'[x], {u[x],u'[x]}, u[x]≥0]

PositivityBasic[a + b u[x] + c ((u[x])^2), {u[x]}]

Issues

non - polynomial expressions

PositivityBasic[b (E^(u[x])) + a ((u[x])^2), {u[x]}]

PositivityBasic[a + (b/((u[x])^2)), {u[x]}]

high order polynomial expressions

PositivityBasic[a + b u[x] + c ((u[x])^2)u'[x] + d ((u[x])^3) + e ((u[x])^4), {u[x],u'[x]}]

Long way

Polynomial cases: General positivity - necessary and sufficient conditions

Positivity[a ((u[x])^3) + b u'[x], {u}, Facts->u[x]≥0, Type->"General"]

Positivity[a + b u[x] + c ((u[x])^2), {u}, Facts->u[x]>= 0, Type->"General"]

Non-polynomial cases: General positivity - sufficient conditions

Positivity[b (E^(u[x])) + a ((u[x])^2), {u}, Facts->u[x]≥0, Type->"General"]

Positivity[b Y[x] + a ((u[x])^2), {u,Y}, Facts->u[x]≥0&&Y[x]≥0, Type->"General"]

However some cases they can be also necessary:

Positivity[a + (b/((u[x])^2)), {u}, Type->"General"]

Polynomial cases: Assymptotic positivity - necessary conditions

Positivity[a + b u[x] + c ((u[x])^2)u'[x] + d ((u[x])^3) + e ((u[x])^4), {u}, Type->"Assymptotic"]

Positivity[a + b u[x] + c ((u[x])^2) + d ((u[x])^3) + e ((u[x])^7), {u}, Type->"Assymptotic"]

Simplify[a + b u[x] + c ((u[x])^2) + d ((u[x])^3) + e ((u[x])^7), %]

PositivityBasic[%, {u[x]}]//Simplify

Finding dissipation laws: examples

Heat equation: (Subscript[u, t])=(a^2)Δu

monomials=Monomials[{u},{x,y},0,2,Generators->{u[x,y],Derivative[1,0][u][x,y],Derivative[0,1][u][x,y],Log[u[x,y]]}]
monomials//Length

monomials=ENLCleaner[monomials,{u},{x,y},{}]
monomials//Length

coefficients=Table[(Subscript[a, k]),{k,Length@monomials}];
generalForm = coefficients.monomials

timeDer = VarD[generalForm ,u,{x,y}]( Derivative[2,0][u][x,y]+ Derivative[0,2][u][x,y])

newTimeDer=CancelNL[timeDer,{u},{x,y}];
newTimeDer=IntegrateByParts[x,y][newTimeDer]

newTimeDer=((u[x,y])^2)newTimeDer//Simplify

constraints=Positivity[-newTimeDer,{u},Type->"General"]

constraints=(constraints//Simplify//BooleanConvert)/.Or->List

constraints//TableForm

Simplify[generalForm,constraints[[4]]]

Map[(%/.#)&,Thread[coefficients->#]&/@IdentityMatrix[Length@coefficients]]//TableForm

Conclusions and further research

Conserved quantities

Dissipated quantities

Overdetermined systems solver

Thank you