# Write a symbolic partial differential equation system with index

I'm trying to find conditions for an endomorphism $L:TM\rightarrow TM$ verfying the condition

$$\oint_{j,k,l} \left(\mathcal{T}_L\right)_{jm}^i \left(\mathcal{T}_L\right)_{kl}^m = \left(\mathcal{T}_L\right)_{jm}^i \left(\mathcal{T}_L\right)_{kl}^m + \left(\mathcal{T}_L\right)_{km}^i \left(\mathcal{T}_L\right)_{lj}^m + \left(\mathcal{T}_L\right)_{lm}^i \left(\mathcal{T}_L\right)_{jk}^m= 0,$$

where $\mathcal{T}_L$ denotes the Nijenhuis torsion of $L$ and a sum over $m$ is understood.

Assuming $\mbox{dim}\,(M)=3$,the system should have $n^3(n-1)/2 =_{n=3} 27$ (independent) equations. So, there is a compact way (maybe the Table command mixed with DSolve) to write the 27 equations quickly in Mathematica? And does the program help me and will it solve symbolically the system giving the endomorphism $L$ in function of some real parameters?

Thanks :)

Pd: The Nijenhuis torsion of $L$ is written in terms of $L_j^i$ and $\frac{\partial L_j^i}{\partial x^k}$, such as show eq. 5 https://arxiv.org/pdf/1405.5118.pdf

EDIT:

Following the @bills' advice, I have decided try with an easier example. Let $f(x,y)=\sin(x+y)$. Then, is very easy check the identity

$$f(x,y) + \frac{\partial^2f}{\partial x\partial y} =0 .$$

I have introduced this partial differential equation in Mathematica using DSolve:

DSolve[f[x, y] + D[D[f[x, y], y], x] == 0, f, {x, y}]


and hoping for a condition about the function $f$, but the program returns the same code.

At this point, what (or where) is the problem? Is a syntax problem? Or maybe is because the function DSolve can not be used to get conditions about a function from a differential equation? Then, it is there some command in Mathematica which could do it?

• Maybe if you start in 2 dimensions, you can figure it out and gain experience towards the larger problem. Probably in 2D you could write out the equations and then see if DSolve (or other commands) can carry out the solution. – bill s Sep 13 '17 at 23:46
• @bills I'll try – Dog_69 Sep 14 '17 at 15:40
• @bills Please see my edit. Thanks. – Dog_69 Sep 14 '17 at 19:47
• @bills No no, don't be confused. My code in Mathematica is only the line DSolve[f[x, y] + D[D[f[x, y], y], x] == 0, f, {x, y}]. The programm doesn't know $f(x,y) = sin(x,y)$. – Dog_69 Sep 18 '17 at 14:59

## Summary:

There are infinitely many functions which satisfy $\left(1+\partial_x\partial_y\right)f(x,y)=0$, that is why Mathematica does not immediately give you the specific answer $sin(x+y)$. For example both $$-J_0\left(2 \sqrt{x} \sqrt{y}\right)\quad,\quad-\frac{2 \sqrt{y} J_1\left(2 \sqrt{x} \sqrt{y}\right)}{\sqrt{x}}$$ are functions that satisfy given differential equation.

## Code to solve homogeneous differential equations symbolically

Following code finds a series expansion for any homogeneous differential equation:

     ClearAll[symbolicDSolve];
symbolicDSolve = Function[{equation, function},
Module[{variables, exponent, exponents, sum, sum2, c, n, a,
equation2, recursion, result, boundary, finalequations,
finalresult},
variables = List @@ function;
exponents = Table[Unique["n"], Length[variables]];
exponent =
If[Length[exponents] == 1, Last[exponents], exponents];
sum[a_ + b_] := sum[a] + sum[b];
sum[x_^(n_ + a_)] :=
sum[x^n] /; And[MemberQ[variables, x], MemberQ[exponents, n]];
sum[b_ x_^(n_ + a_)] :=
sum[(b /. n -> n - a) x^n] /;
And[MemberQ[variables, x], MemberQ[exponents, n]];
simplify = {sum2[a_ ] + sum2[b_ ] :> sum2[a + b],
a_ sum2[b_] :> sum2[a b]};

equation2 = (equation //.
function :>
sum[(c @@ exponents) (Times @@ (variables^exponents))]) //.
sum[(D[(c @@
exponents) (Times @@ (variables^exponents)), ##] & @@
Transpose[{variables, {a}}])];
recursion =
Simplify[
Coefficient[
Last[equation2 /. sum -> sum2 //.
simplify], (Times @@ (variables^exponents))]];
result =
Evaluate[(Inactive[
Sum][(c @@
exponents) (Times @@ (variables^exponents)), ##] & @@
Table[{exponents[[i]], 0, Infinity}, {i,
Length[variables]}])] /.

Flatten[{RSolve[recursion == 0, (c @@ exponents), exponent]}]

]
];


## Analysis of the original question:

   In:= symbolicDSolve[
D[f[x, y], x, y] + f[x, y], f[x, y]]

Out= Inactive[Sum][((-1)^(-1 + n13) x^n13 y^n14 C[n13 - n14])/(
Pochhammer[2, -1 + n13] Pochhammer[2 - n13 + n14, -1 + n13]), {n13,
0, \[Infinity]}, {n14, 0, \[Infinity]}]


which is simply

$$\underset{\text{n13}=0}{\overset{\inf }{\sum }}\underset{\text{n14}=0}{\overset{\inf}{\sum }}\frac{(-1)^{\text{n13}-1} x^{\text{n13}}y^{\text{n14}}c_1(\text{n13}-\text{n14})}{(2)_{\text{n13}-1}(-\text{n13}+\text{n14}+2)_{\text{n13}-1}}$$

for undetermined coefficients $c(n13-n14)$. We see that differential equation as a constraint reduced the two dimensional unknown $c(n13,n14)$ to one dimensional unknown $c(n13-14)$ ; still, there are infinitely many solutions! Let us assume that $c(0)=1$ and all other $c(i)=0$. That means, out solution is $$f(x,y)=\underset{\text{a}=0}{\overset{\infty }{\sum }}\frac{(-1)^{a+1} c_1(0) x^a y^a}{\Gamma (a+1)^2}$$ We can let mathematica calculate this for us, we immediately get the Bessel function:

In:=Sum[((-1)^(1 + a)x^ay^aGammaC)/(Gamma[1+a]Gamma[1+a]),
{a,0,Infinity}]Out= -BesselJ[0, 2 Sqrt[x] Sqrt[y]] C


Of course we can directly check that the result satisfies original differential equation:

In:= (Derivative[1, 1][f][x, y] + f[x, y]) /.
f -> Function[{x, y}, -((2 Sqrt[y] BesselJ[1, 2 Sqrt[x] Sqrt[y]])/
Sqrt[x])] // FullSimplify

Out= 0


## Other examples of the code:

First order, one variable: $a f'(x)+b f(x)=0$:

     In:= symbolicDSolve[a D[f[x], x] + b f[x], f[x]] // Activate
Out= -((a E^(-((b x)/a)) C)/b)


$$f(x)=-\frac{a c_1 e^{-\frac{b x}{a}}}{b}$$

Second order, one variable: $f''(x)+a x f'(x)+ b f(x)=0$:

    In:=symbolicDSolve[f''[x]+axf'[x]+bf[x],f[x]]
//Activate//FullSimplify
Out= (E^(-(1/2) x (a x +
Sqrt[-4 b + a^2 x^2])) (Sqrt[-4 b +
a^2 x^2] (C - E^(x Sqrt[-4 b + a^2 x^2])
C)-   a x (C + E^(x Sqrt[-4 b + a^2 x^2]) C)))


$$f(x)=\frac{e^{-\frac{1}{2} x \left(\sqrt{a^2 x^2-4 b}+a x\right)} \left(\sqrt{a^2 x^2-4 b} \left(c_1-c_2 e^{x \sqrt{a^2 x^2-4 b}}\right)-a x \left(c_2 e^{x \sqrt{a^2 x^2-4 b}}+c_1\right)\right)}{b}$$

First order, three variables: $(\partial_x+z \partial_z)f(x,y,z)=0$:

  In:= symbolicDSolve[
D[f[x, y, z], {x, 1}] + D[f[x, y, z], {y, 1}],
f[x, y, z]] // FullSimplify

Out=
Inactive[Sum][(
x^n252 y^n253 z^
n254 Pochhammer[1 - n252 - n253, -1 + n252] C[n254][
n252 + n253])/
Gamma[1 + n252], {n252, 0, \[Infinity]}, {n253,
0, \[Infinity]}, {n254, 0, \[Infinity]}]


$$f(x,y,z)=\underset{\text{n252}=0}{\overset{\infty }{\sum}}\underset{\text{n253}=0}{\overset{\infty }{\sum}}\underset{\text{n254}=0}{\overset{\infty }{\sum}}\frac{x^{\text{n252}} y^{\text{n253}} z^{\text{n254}} (-\text{n252}-\text{n253}+1)_{\text{n252}-1} c_1(\text{n254})(\text{n252}+\text{n253})}{\Gamma (\text{n252}+1)}$$

Among doubly infinite solution space, let us choose the following one: $$c_1(n)(k)=\frac{1}{n!}e^{-k}$$. Then, we get

In:=
Sum[Sum[(x^n249 y^(k - n249) z^n251 Pochhammer[1 - k, -1 + n249])/
Gamma[1 + n249] Exp[-k]/(n251 !), {n251, 0, Infinity}, {k, 0,
Infinity}], {n249, a, Infinity}] // FullSimplify

Out= E^z x^a y^-a Gamma[a] Hypergeometric2F1Regularized[1, a,
1 + a, x/y]


which reads $$f_a(x,y,z)=e^z x^a y^{-a} \Gamma (a) \, _2\tilde{F}_1\left(1,a;a+1;\frac{x}{y}\right)$$

We summed over one of the parameters not from 0 but from $a$ to $\infty$ as the summation is not convergent at the origin around which we have made our series expansion within the code.

If we check the action of our differential operator on our result, we see that

$$(\partial_x+z \partial_z)f_a(x,y,z)=e^{z}x^{-1+a}y^{-a}\sim f_{a-1}(x,y,z)$$ in terms of its divergence around $x\sim 0$. As it stands, it satisfies the differential equation in leading order in regions away from the origin. We can also improve our code such that series expansion point can be chosen manually to avoid singular points.

• I don't know if your answer will solve my problem completely. Or even if I can understand it really well. But I think it is brilliant. I thank to your patience writing such a detailed answer. I have voted up of course. – Dog_69 Mar 20 '18 at 0:20

I don't see how to get the closed form solution. But it's possible to do it numerically. Imposing some arbitrary boundary conditions:

sol = NDSolve[{f[x, y] == -D[D[f[x, y], y], x], f[0, y] == 0,
f[x, 0] == Sin[x], f[x, 5] == 0}, f, {x, 0, 10}, {y, 0, 5}]

Plot3D[Evaluate[f[x, y] /. sol], {x, 0, 10}, {y, 0, 5},
PlotRange -> All] 