Formal Power Series Solutions of AODEs

This program is accompanying the paper “Formal Power Series Solutions of Algebraic Ordinary Differential Equations” written by Sebastian Falkensteiner, Thieu Vo and Yi Zhang.
author of this notebook: Sebastian Falkensteiner

Problem

In general we want to test wheter a differential polynomial is m-th non-vanishing (see the definition in the paper) or not.

In the last section we also do a statistical evaluation of Kamke’s collection. In the paper, we present a more detailed result as one can obtain by running the code implemented here.

ATTENTION: Running the last section “Statistical Evaluation” can be time consuming.

Functions

ClearAll[y,x,c,t];

DiffOrder[F]
Input: F... differential polynomial
Output: order of F as natural number

ClearAll[DiffOrder];
DiffOrder[F_]:=Max[Cases[F,Derivative[a_][y][x]->a,{0,Infinity}]];

TotalDerivative[F,k]
Input: F... differential polynomial, k... natural number
Output: k-th derivative of F with respect to x

ClearAll[TotalDerivative];
TotalDerivative[F_,0]:=F;
TotalDerivative[F_,k_Integer]/;k>0:=D[F,{x,k}];

PartialDerivative[F,k]
Input: F... differential polynomial, k... natural number
Output: partial derivative of F with respect to y^{(k)} (or 0 if k is negative)

ClearAll[PartialDerivative];
PartialDerivative[F_,k_Integer]/;k<0:=0;
PartialDerivative[F_,k_Integer]:=D[F,D[y[x],{x,k}]];

Help functions for evaluation:

ClearAll[InitialCoeffRules,EvaluateCoeffRules];
InitialCoeffRules[k_]:=Table[D[y[x],{x,i}]->c[i],{i,0,k}];
InitialCoeffRules[k_,coeffs_]:=InitialCoeffRules[k]/.Table[c[i]->coeffs[[i+1]],{i,0,k}];
EvaluateCoeffRules[coeffs_]:=MapIndexed[c[First[#2]-1]->#1&,coeffs];

ClearAll[Evaluation];
Evaluation[exp_,coeffrules_]:=(exp/.coeffrules)/.{x->0};

SeparantMatrix[F,m]
Input: F... differential polynomial, m... natural number
Output: the non-evaluated separant matrix of F as it is defined in the paper

ClearAll[SeparantMatrix];
SeparantMatrix[F_,m_]:=Module[{n=DiffOrder[F],matrix=ConstantArray[0,{m+1,m+1}]},Clear[c];Table[matrix[[j,i]]=TotalDerivative[PartialDerivative[F,n+j-i],j-1],{i,1,m+1},{j,1,m+1}];Evaluation[matrix,InitialCoeffRules[n+Max[n,m]]]];

SeparantVector[F,m (,c)]
Input: F... differential polynomial, m... natural number, optional: coeffs... list of initial values
Output: the separant matrix of F multiplied with the coefficient vector comparible to Theorem 4.6. Either non-evaluated or evaluated at c if they are given as input

ClearAll[SeparantVector];
SeparantVector[F_,m_]:=Table[Binomial[2m+1+t,i],{i,0,m}].SeparantMatrix[F,m];
SeparantVector[F_,m_,coeffs_List]:=Module[{},Clear[c];SeparantVector[F,m]/.EvaluateCoeffRules[coeffs]];

JetIdeal[F,k]
Input: F... differential polynomial, k... natural number
Output: k-th Jet Ideal of F as a list

ClearAll[JetIdeal];
JetIdeal[F_,k_]:=Module[{n=DiffOrder[F]},Clear[c];Evaluation[Table[TotalDerivative[F,i],{i,0,k}],InitialCoeffRules[n+k]]];

JetCoeff[F,k]
Input: F... differential polynomial, k... natural number
Output: the Ideal generated by the entries of the k-th separant matrix of F

ClearAll[JetCoeff];
JetCoeff[F_,k_]:=Module[{n=DiffOrder[F]},GroebnerBasis[Flatten[SeparantMatrix[F,k]],Reverse[Table[c[i],{i,0,n+k}]],CoefficientDomain->RationalFunctions]];

TestIdeal[F,k]
Input: F... differential polynomial, k... natural number
Output: the simplified ideal consisting the (2k)-th Jet Ideal and the ideal generated by the entries of the k-th separant matrix of F

ClearAll[TestIdeal];
TestIdeal[F_,k_]:=Module[{n=DiffOrder[F]},Simplify[GroebnerBasis[Union[JetIdeal[F,2*k],JetCoeff[F,k]],Reverse[Table[c[i],{i,0,n+k}]],CoefficientDomain->RationalFunctions]]];

TestNonVanishing[F,k]
Input: F... differential polynomial, N... natural number
Output: a natural number k such that F is k-th non-vanishing if the upper bound N has not been reached first (Then the output is “-1”).

ClearAll[TestNonVanishing];
TestNonVanishing[F_,N_]:=Module[{ideal=TestIdeal[F,0],k=0},While[k≤N&&!ideal==={1},ideal=TestIdeal[F,k];k++];If[k==0,0,If[k==N+1,-1,k-1]]];

The following functions are used for the statistical evaluation:

ClearDenominator[F]
Input: F... differential polynomial
Output: F multiplied with its common denominator

ClearAll[ClearDenominator];
ClearDenominator[F_]:=Simplify[F*Denominator[Together[F]]];

TestPolynomial[F]
Input: F... differential polynomial
Output: “True” if F is polynomial in x,y[x],...,y^{(n)}[x] and “False” otherwise

ClearAll[TestPolynomial];
TestPolynomial[F_]:=Module[{n=DiffOrder[F],var},var=Union[{x},Table[D[y[x],{x,i}],{i,0,n}]];PolynomialQ[F,var]];

Examples

A first non-vanishing AODE which separant is always vanishing.

F1:=x*y'[x]+y[x]^2-y[x]-x^2;

TestNonVanishing[F1,5]

MatrixForm[SeparantMatrix[F1,1]]

A vanishing differential polynomial (first-order, autonomous).

F2:=y'[x]^2+y[x]^3;
TestNonVanishing[F2,10]

TestIdeal[F2,5]

We slightly modify F2. For a generic value of a the AODE is non-vanishing, only for a=0 we obtain the previous vanishing differential polynomial.

TestNonVanishing[a*x+y'[x]^2+y[x]^3,10]

For every natural number m we obtain a differential polynomial with non-vanishing number equal to m.

F3:=(y'[x]^2+y[x]^2)/2+x^(2*m);

m=4;
TestNonVanishing[F3,10]

m=5;
TestNonVanishing[F3,10]

F4:=x+y[x]-x* y[x]+y[x]^2-y[x]^3+(x-(-1+x)*y[x]+(1+x)*y[x]^2)*y'[x];

TestNonVanishing[F4,10]

SeparantVector[F4,2,{0,-1,-4,-39}]

Compute the separant vector at the problematic point of the following AODE:

F5:=y'[x]^2+y'[x]-2*y[x]-x;
TestNonVanishing[F5,10]

SeparantVector[F5,1,{-1/8,-1/2,0}]

F6:=y'[x]*y''[x]^2+2*x^3*y[x]+x^2;
TestNonVanishing[F6,10]

MatrixForm[SeparantMatrix[F6,2]]

SeparantVector[F6,2,{c0,-1/c3^2,0,c3,c4}]

We cannot find a bound for the roots of the first non-zero entry of the separant vector of the following AODE:

F7:=x*y''[x]+y'[x]*y[x];
TestNonVanishing[F7,10]

SeparantVector[F7,1,{c0,0,0,c3,c4}]

A 1-st non-vanishing differential polynomial which corresponding curve is not smooth.

F8:=x*y'[x]+y[x]^2;
TestNonVanishing[F7,10]

Statistical Evaluation

Import Kamke’s collection of single ODEs.
IMPORTANT: Change the import adress according to the file location of the file “Kamke_ODE.m” .

ClearAll[ode];
ode=ConstantArray[0,1855];
string = << Kamke_ODE.m;

Select the irreducible AODEs and clear denominators.

ode=ode//.{a_==b_:>a-b};

ode=Map[ClearDenominator[#]&,ode];

ode=ode/.{0->Nothing};

ode=Select[ode,TestPolynomial[#]&];

ode=Quiet[Select[ode,IrreduciblePolynomialQ[Factor[#]]&]];

Test up to non-vanishing number 4 and insert vanishing number into list.
vanishing number = -1 iff F has a vanishing number of 4 or higher or is vanishing.

ode=Map[{#,{DiffOrder[#],TestNonVanishing[#,4]}}&,ode];

Return the statistic evaluation in the following shape:
{n,m}->x where x is the number of m-th non-vanishing AODEs of order n (note: n= -1 if the AODE is possibly vanishing)

Counts[Map[#[[2]]&,ode]]

count=Counts[Map[#[[2]][[2]]&,ode]]

Percentage of possibly vanishing AODEs in Kamke’s collection.

N[count[-1]/(Total[count]-count[-1])]