Solving a system of nonlinear PDEs by minimization

Question

I have two coupled nonlinear partial differential equations of the form:

$ \begin{align} \dot{u} -f(u,u',u'',v,v',v'')=0 \\ \dot{v} -g(u,u',u'',v,v',v'')=0 \end{align} $

The boundary conditions are such that I know $u$ at time $t=0$ and $v$ at time $t=-\infty$. I also know the values of both functions at the two boundaries of coordinate interval. I need the solution for $x \in [0,1]$ and $t \in (-\infty,0]$ For a specific parameter value I also know an analytical solution and this allows me to tune this parameter little by little and then solve the equations by using as an initial guess for solution my solution for previous parameter value.

I'm solving this by writing $ S \equiv (\dot{u} -f(u,u',u'',v,v',v''))^2 + (\dot{v} -g(u,u',u'',v,v',v''))^2 $ and then finding the minimum of S and verifying that it is approximately zero. Currently I'm doing the minimization by discretizing $u$ and $v$ on my time and coordinate grid and then using conjugate gradient algorithm from SciPy optimization library. This gives me good results, but is terribly slow for larger grids. For a grid with $10^4$ grid points my solver is already not feasible.

My questions:

1) Which minimization algorithm is best for this sort of task? Conjugation gradient seems to require a lot of steps to converge to the solution.

2) What external libraries there are for solving this kind of problem? I'd prefer ones with a Python API. I looked at some finite element libraries and they seem like a bit of an overkill for this problem.

Answer 1

1) If you're just looking to solve the PDEs without any other optimization, then my answer would be "none of them". Algorithms that discretize partial differential equations and then solve them as algebraic equations are massively parallelizable. It is possible to solve a PDE over a billion point mesh. Algorithms for nonlinear programming have made great advances recently. It is now possible to use preconditioned iterative algorithms somewhat effectively to solve large-scale nonlinear programs with interior point methods, but I don't know of any nonlinear program that has been solved with a billion degrees of freedom, once we exclude special structure (typically, stochastic programs).

There is an approach to PDE-constrained optimization called simultaneous analysis and design (SAND, also called "all-at-once") that treats the PDEs as constraints to be solved within an optimization formulation. Solution of large-scale problem instances typically hinges on formulating a good preconditioner for the KKT system. This approach is different from your least-squares approach, and without having an optimization problem using your PDEs as constraints, I do not recommend this approach.

2) If you're dead set on doing this with optimization, I'd try pyipopt, which is a Python interface to the optimization library IPOPT, which is currently the best open-source nonlinear programming solver.

Really, though, I'd recommend using traditional discretization approaches for solving PDEs, for which there are many good Python libraries available (depending on what you want to do, any of PyClaw, FEniCS, FiPy, sfepy, firedrake, petsc4py, PyTrilinos, and so on...).