pyomo
pyomo copied to clipboard
Imaginary numbers in collocation points
Summary
When trying to discretize a PDE with collocation using a large amount of points (20), imaginary numbers occur in the returned points. This can throw an error further down the road.
Steps to reproduce the issue
from pyomo.dae.plugins.colloc import calc_cp
print(calc_cp(1, 0, 19))
[(0.003623538921546554+0j), (0.01889025162098419+0j), (0.04617986689257604+0j), (0.0843179433375663+0j), (0.13322082232455726+0j), (0.19039692429875266+0j), (0.2564281921284097+0j), (0.32602096766337535+0j), (0.40613716185242277+0j), (0.47264007102384903+0j), (0.5908680381819745-0.01669040549711183j), (0.5908680381819745+0.01669040549711183j), (0.7378135199125857-0.04910048600435542j), (0.7378135199125857+0.04910048600435542j), (0.8685143860954003-0.0520287893947382j), (0.8685143860954003+0.0520287893947382j), (0.9624876846319255-0.030883704390630654j), (0.9624876846319255+0.030883704390630654j), (0.9991872587024497+0j)]
Stepping through the debugger, I get the impression that this might not be the best-conditioned calculation:

Since root finding of high-order polynomials is also a badly-conditioned process, having bad results occur isn't surprising, but I wanted to at least put this issue on record in case anybody else encounters it.
@dallan-keylogic thanks for the bug report. I think this hasn't been encountered before because when we increase the discretization resolution for the collocation approach we typically increase the number of finite elements rather than the number of collocation points. I don't know of any work using the collocation approach that uses more than 3-6 collocation points.
If you think it's useful to support much larger numbers of collocation points I can rework the calculation you pointed to above. Otherwise, it would be easy to set an upper bound on the number of collocation points supported by the discretization transformation.
I was intrigued by the equivalence of collocation and spectral methods and the possibility of "spectral accuracy" with a large number of elements (although perhaps what I should have thought is because of "spectral accuracy" I don't even need 20 points).
Capping the number of points is acceptable.