Differential equations with free variables

[vsht]

Hi,

as to my understanding (and according to Issue #11), the Python version of Fuchsia can deal with systems that depend on a free variable, with the tacit assumption that the variable doesn't take values that lead to singularities.

I have a small example (3x3 matrix), where this works nicely and the epsilon form is found straightforwardly

{{(4 - 6*eps - 3*z + 2*eps*z)/(2*(-1 + z)*z), (-1 + 2*eps)/(2*(-1 + z)), 0}, 
 {-(-1 + 2*eps)/(2*(-1 + z)), -(-2*eps + z + 2*eps*z)/(2*(-1 + z)*z), 0}, 
 {-(((-1 + 2*eps)*z)/((-1 + z)*(1 + 2*c*z))), 
  -(((-1 + 2*eps)*(-2 + z))/((-1 + z)*(1 + 2*c*z))), 
  -((-eps - 2*c*z + 2*eps*z + 2*c*eps*z + 2*c*z^2)/
    ((-1 + z)*z*(1 + 2*c*z)))}}

However, with another similar matrix

{{-((-2 + 3*eps + z - 2*c*z + 4*c*eps*z)/((-1 + z)*z)), 
  ((-1 + 2*eps)*(1 + 4*c + 4*c^2*z))/(2*(-1 + z)*z), 0}, 
 {(-2*(-1 + 2*eps))/(-1 + z), (eps - z - 2*c*z + 4*c*eps*z)/((-1 + z)*z), 
  0}, {(-4*(-1 + 2*eps)*z)/((-1 + z)*(1 + 2*c*z)), 
  (2*(-1 + 2*eps))/(-1 + z), 
  -((-eps - 2*c*z + 2*eps*z + 2*c*eps*z + 2*c*z^2)/
    ((-1 + z)*z*(1 + 2*c*z)))}}

Fuchsia claims that it has found the canonical form, but the resulting matrix

{{-2*(2*c*eps-c-2)/(z-1)-eps/z,4*(2*c*eps-c-2)/z+1/2*((16*c^2-8*c+1)*eps^2+4*c^2-4*(4*c^2+7*c-2)*eps+16*c+16)/(eps*(z-1)),0},{-2*eps/(z-1),2*((2*c-1)*eps-c-2)/(z-1)+3*eps/z,0},{-4*eps/(c*(z-1)),2*((4*c-1)*eps-2*c-4)/(c*(z-1))+4*eps/(c*z),4*c*eps/(2*c*z+1)-eps/(z-1)-eps/z}}

is obviously not of the form eps*B. I wonder why that happens and if one can regard this as a bug: I wouldn't complain if Fuchsia would simply fail, but finishing successfully and returning something that is not correct, is somewhat odd.

If the free variable c is replaced by some prime number, Fuchsia will return the proper epsilon form, so that the matrix itself doesn't seem to be problematic.

Cheers, Vladyslav

[magv]

Your understanding is correct, and this example looks like a giant bug...

I'll look into it.

In the mean while, you can use the new fuchsia. It produces this epsilon form matrix:

{{-eps*z^(-1)-2*eps*(-1+z)^(-1)*(-2+c),-8*z^(-1)*(2*eps-eps*c)+(25*eps+4*eps*c^2-20*eps*c)*(-1+z)^(-1),0},
{-eps*(-1+z)^(-1),-2*(3*eps-eps*c)*(-1+z)^(-1)+3*eps*z^(-1),0},
{eps*(-1+z)^(-1)*c^(-1),-2*eps*z^(-1)*c^(-1)+(5*eps-2*eps*c)*(-1+z)^(-1)*c^(-1),4*eps*(2*z+c^(-1))^(-1)-eps*z^(-1)-eps*(-1+z)^(-1)}}

... using this transformation:

{{2*eps*z^(-1),-2*(1+2*z*c-3*eps*z-4*eps-2*eps*z*c)*z^(-2),0},
{0,4*(-1+2*eps)*z^(-1),0},
{0,0,-8*(-1+2*eps)*(1+2*z*c)^(-1)*c}}

[vsht]

Thanks for the prompt reply! Indeed the new fuchsia does the job pretty well here :+1:

I would be glad to use it more often, but unfortunately it seems to have some limitations when working with matrices that have mixed singularities in the free variable and eps. Like in this 2x2 matrix

{{(-2*(-1 + 7*eps - 10*eps^2 + y^2 + eps*y^2 - 10*eps^2*y^2 + 2*y^4 - 
       14*eps*y^4 + 20*eps^2*y^4))/((-1 + y)*
     y*(1 + y)*(-1 + eps - y^2 + 5*eps*y^2)), (-4*eps*
     y)/((-1 + y)*(1 + y)*(-1 + eps - y^2 + 
       5*eps*y^2))}, {(-4*(-1 + 3*eps)^2*(-1 + 2*eps + 
       4*eps*y^2))/(y^3*(-1 + eps - y^2 + 5*eps*y^2)), (2*(-1 + 
       4*eps - 7*eps^2 - 2*eps*y^2 + 2*eps^2*y^2 - y^4 + 4*eps*y^4 + 
       5*eps^2*y^4))/((-1 + y)*
     y*(1 + y)*(-1 + eps - y^2 + 5*eps*y^2))}}

which yields

terminate called after throwing an instance of 'std::runtime_error'
  what():  pfmatrix(): can't solve equations of 2nd degree or higher

It can work in simple cases, where the mixed polynomial is only linear in y' andepsbut higher order polynomials seem to be a show stopper. CANONICA can deal with such systems since Meyer's algorithm is based on expanding the matrix ineps`. For example

RecursivelyTransformSectors[{mat}, {y}, {{1, 2}}, {1, 1}, 
 TDeltaDenominatorDegree -> 3, TDeltaNumeratorDegree -> 3]

returns

{{{2*eps*(-1 + y^2), (-1 + y^2)*(-y^2 + eps*(-4 + 5*y^2))}, {((1 - 3*eps)^2*(-1 + y^4))/y^2, 
   (-2*(1 - 3*eps)^2*(-1 + y^2))/y^2}}, 
 {{{(2*eps*(-3 + y^2))/(y*(-1 + y^2)), (4*eps*(5 - 4*y^2))/(y*(-1 + y^2))}, 
   {(4*eps)/(y - y^3), (4*eps*(3 - 2*y^2))/(y*(-1 + y^2))}}}}

Unfortunately, CANONICA has many other limitations, so I'd be more than happy to switch to the CPP version of fuchsia if my cases would be covered. However, I perfectly understand the difference between a dev preview and a final version, so I'm not complaining. Just would like to point out that this would be a nice feature for the stable release of the new fuchsia.

[magv]

Yes, that is the main limitation of the new version. I hope this limitation to be temporary, but I can't say when will it be fully lifted (it is a coding problem, not a problem in principle).

The somewhat good news is that the 1/(-1 + eps - y^2 + 5*eps*y^2) pole in your matrix is spurious. In fact, it can't physical, just look at it! Spurious poles can be completely removed, and do not appear in the epsilon form. In that sense it's the easy case. I'll try to make it work with the new Fuchsia.

I've previously removed such poles manually by 1) changing the master selection (with guesswork), and/or 2) making sure that my IBP gives the least possible number of masters (for example FIRE is known to give too many at times; I've had better luck with Kira). But there is an algorithmic way of removing these poles too.

[vsht]

I fully agree with you that such poles can be usually removed by choosing a better basis. But since this procedure needs to be done by hand, it may take quite some time on larger systems or when there are many systems to deal with. So essentially it becomes a question of saving time: Can I still get away with a "bad basis" if my tool for finding the canonical form can handle that, or do I absolutely have to get my hands dirty and fix the basis.

I surely understand what you mean with a coding problem, as we all have only a finite amount of time to spend on developing tools. It is already good enough that Fuchsia is still being developed and maintained, unlike CANONICA.