search

sympy / sympy

A computer algebra system written in pure Python
https://sympy.org/
Other
12.85k stars 4.41k forks source link

Numerically solving system of PDAEs using polynomial power series #15015

Open Foadsf opened 6 years ago

Foadsf commented 6 years ago

I requested this feature firstly here in the Mailing list and this Reddit post, and then followed by this post on cs.stackexchange, which was kindly responded by @asmeurer.

Basically I'm looking for a method similar to Mathematica's AsymptoticDSolveValue which takes an ODE and initial conditions, plus a positive integer d, and returns a power series (i.e. polynomial of the degree d) as approximation of the solution (more info here). I think it is feasible to implement a function to solve system of partial differential equations. A function:

AsymptoticPdeSolve(eqns,fs,vars,D)

Where eqns is the set of symbolic partial differential and algebraic expressions plus boundary conditions, fs are the set of functions we want to solve, vars are the variables and D is the dimension of our multivariate polynomials.

P.S.

  1. Maple's dsolve also has a series option.
Foadsf commented 6 years ago

@celliern was so kind to implement one version of the idea here. I haven't checked the code yet, but I have some other points to add:

1. The PDAEs are not necessarily always linear. So maybe it is a good idea to first expand them using Taylor/Maclaurin power series. But before that we might need to limit the PDAEs to polynomial form, where each monomial are in the form of y1*...*yn or y1 * y1' *y1" ... 2. In some cases implementation of boundary conditions results in zero coefficients, or they contradict the PDEs. in these cases I think it is better minimize the integration of the error in the boundary. For example if ∂𝛺1 is part of the boundary where BC1 applies, and P is the symbolic polynomial compatible with the PDEs, then ∮_{∂𝛺1} (BC1(P))^2 ds should be minimized. ( [BFGS](https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm) might be useful here )
certik commented 6 years ago

It would be also useful to first implement an equivalent AsymptoticDSolveValue into SymPy, that would be valuable on its own.

Regarding AsymptoticPdeSolve the best is to take some example and code a prototype, see if it works, then improve upon it.

Foadsf commented 6 years ago

It would be also useful to first implement an equivalent AsymptoticDSolveValue into SymPy, that would be valuable on its own.

Great point. Totally agree. I just think PDEs are a completely different beasts, and AsymptoticPdeSolve won't be an easy extension to AsymptoticDSolveValue.

Regarding AsymptoticPdeSolve the best is to take some example and code a prototype, see if it works, then improve upon it.

We now have one implementation by @celliern here and some ideas by me and @asmeurer over here.

Foadsf commented 6 years ago

I'm here discussing with the FriCAS team about the same idea, apparently there is a seriesSolve in Axiom-FriCAS. I will fins the source code and link in here.

P.S.1. The source code of FriCAS seriesSolve can be found here and here

P.S.2. Surprisingly enough there are some traces of the seriesSolve function by @manoj-kumar-s in Sympy mailing list over here and also on his blogpost back in 2013!

Foadsf commented 6 years ago

Whitepapers:

Here I will list all the relevant papers, books, posts and keywords I can find around:

  1. Aroca, Fuensanta, J. Cano, and Françoise Jung. "Power series solutions for non-linear PDE's." Proceedings of the 2003 international symposium on Symbolic and algebraic computation. ACM, 2003.
  2. the wikipedia page on this matter is really insufficient!
  3. Cauchy–Kowalevski theorem has been pointed to me both in Sympy mailing list and here on Reddit.
  4. Frobenius method
  5. Bender, Carl M., and Steven A. Orszag. Advanced mathematical methods for scientists and engineers I: Asymptotic methods and perturbation theory. Springer Science & Business Media, 2013.
  6. Chebyshev polynomials
  7. Spectral method --> Chebfun and DEDALUS
  8. Formal power series
  9. Runge's phenomenon
  10. Asymptotic expansion
  11. Cauchy–Hadamard theorem
  12. WKB approximation
Foadsf commented 6 years ago

Considering the fact that multiplication of long multivariate polynomials can be inefficient in sympy, I'm developing a multivariate polynomial class using numpy. please take a look at it over here and let me know how do you think.

P.S. There are also other libraries: POLYNOMIAL, LibPoly and berlin seidler and matthewhr

jksuom commented 6 years ago

There are two polynomial implementations in SymPy, dense and sparse. The former can be inefficient for many variables but not the latter one.

Foadsf commented 6 years ago

@jksuom If the later is more efficient than transforming the polynomial to a numpy class and then numpy ndarray multiplication and then changing back to the sympy, then it is ok. How can I use the later?

I'm also considering doing the differentiation over there. maybe I can implement the linear PDEs also there. Then it might be faster I think. not sure though.

jksuom commented 6 years ago

The implementation is in polys.rings. There are several public interfaces, ring, xring, vring and sring. The last one is probably most convenient for constructing both rings and polynomials.

Foadsf commented 6 years ago

I'm surprised nobody pointed out that SymPy already has the ability to solve ODEs using formal power series! dsolve(eq, hint='1st_power_series') examples here.

certik commented 6 years ago

Good. So you are halfway done already. ;)

Foadsf commented 6 years ago

I'm surprised again to find out there has been already attempts to implement sequences and formal power series in SymPy, and nobody has pointed it so far!

One can generate sequenced symbolic formal power series using Sequence and PowerSeries of that branch. If I knew this exist I wouldn't have to do it again (here, here and here)