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.

  • $\begingroup$ Could you elaborate on the exact form of $f$ and $g$? That would determine the character of the PDE (hyperbolic, parabolic, totally bizarre) and in turn which solution method is best. $\endgroup$ May 21, 2014 at 15:56
  • $\begingroup$ It is basically a fourth order polynomial plus a "global" part which is proportional to the average value of $v$ over the coordinate axis times $v$ itself. Because of the global term, the problem maybe goes into the category "totally bizarre". $\endgroup$
    – Echows
    May 22, 2014 at 8:30

1 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...).

  • $\begingroup$ I solve a simplified version of this problem by just discretizing the problem, starting from an initial condition and propagating the solution in time. I didn't figure out how to do this for this problem for two reasons: 1) I have both initial AND final conditions for the solution. 2) I have the "global" term in the equations (see my comment above), i.e. the solution for one point doesn't depend only on the few neighboring points, but the solution everywhere in the coordinate axis. I will, however, look into the libraries you recommended, thanks for that. $\endgroup$
    – Echows
    May 22, 2014 at 12:56
  • 1
    $\begingroup$ You have a boundary value problem. One way to solve such a problem is to treat the time coordinate as if it were a spatial coordinate, mesh with a structured grid, and solve. You could also use shooting-type methods (use multiple shooting, not single shooting). $\endgroup$ May 22, 2014 at 18:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.