I am researching heat diffusion in an optical element irradiated by laser. This problem is described by the PDE which I wrote down in this question. I am using an implicit numerical scheme to model the heat distribution evolving in time.

After some time has passed the system comes to a stationary distribution (corresponding with thermal equilibrium). One of the most interesting features of the model is the nonuniqueness of stationary state: if we'll start from another initial distribution, then we may come to another stationary state.

So, one of the task I need to solve is to find the set of stationary states for a given set of parameters.

I am solving it with following algorithm: I run the simulation until the heat distribution stops changing (difference between current and past distribution is lesser than given epsilon). So, after one stationary distribution is found, I increase and decrease laser power to find the next stationary state. By varying initial distributions for given set of parameters I compute a set of stationary distributions.

But, there are two problems. First, to compute the needed set with given accuracy I need to use small steps in space and time, small value of epsilon used to find stationary states and many steps of laser power variating. This is very computationally intensive.

The second problem is that it's easy to miss a distribution or count one as two, and I am not sure if there are other distributions which my algorithm cannot found.

So, what algorithms should I consider to use instead of just using implicit scheme and waiting until equilibrium is estabilished? There are two significant cases:

  • there is one possible stationary state, and the task is just to find it with given accuracy and minimal amount of computation

  • there is a set of possible stationary distributions and the task is to find them all.


1 Answer 1


If there is one steady state, then a standard method is to use pseudotransient continuation, which is essentially linearly implicit Euler with residual-based step size adaptivity. This approach has relatively strong convergence guarantees for a quadratically convergent method. Coffey, Kelley, and Keyes (2003)

Methods to provably compute all solutions are rarely available, though Baron can sometimes be used.

If the steady state has manifold structure parametrized by a low-dimensional space, then you can explore it using continuation techniques, see the book by Allgower and Georg or by Seydel. Some software is available for this task, e.g. Multifario.

  • 2
    $\begingroup$ Can't you sometimes phrase this as an eigenproblem, which would make getting all solutions somewhat tractable? $\endgroup$ Dec 6, 2011 at 17:18
  • 1
    $\begingroup$ Yes, but this problem is nonlinear and I got the impression that the presence of multiple solutions is due to nonlinearity instead of linear degeneracy. In the linear case, the space of steady state solutions is the space spanned by all eigenvectors associated with eigenvalues equal to zero. $\endgroup$
    – Jed Brown
    Dec 6, 2011 at 17:33
  • $\begingroup$ A general principle for finding stationary solutions is that the accuracy of the time stepping is not that important; you don't care about how you get to the stationary solution, only that you get there (exception: if there are multiple stationary solutions and you want to find the right one). That's why it is okay to use a first-order method like Euler and fairly big time steps. $\endgroup$ Dec 7, 2011 at 11:29

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.