Say we have a discretised a coupled nonlinear system of two PDEs to give a system of ODEs which approximates the original system, $$ \frac{\partial u}{\partial t} = F_1(t,\boldsymbol{u,v}) \\ \frac{\partial v}{\partial t} = F_2(t,\boldsymbol{u,v}) \\ $$ where $\boldsymbol{u} = [u_1 \cdots u_N]^{T}$ and $\boldsymbol{v} = [v_1 \cdots v_N]^{T}$ are a vector of solution variables with the same number of elements as discretised points in the system. $F$ is a vector which depends on how the PDEs has been discretised.

When PDEs are left in this semi-discrete form (i.e. the transient term is not discretised) they can be solved using conventional ODE solvers, this technique is known as the Method of Lines (MOL).

Normally the ODE solver can be used to step the system forward in time until a desired time point or until the system reaches steady-state. However, can the MOL technique be used to solve directly for the steady-state value?

Internally the ODE solver will be computing a Jacobian and using Newton's method to solve the system, these are the prerequisites for solving the system directly for steady state. So it would seem that it has all the necessary information to do so. Moreover, to solve the above system for the steady state one would apply Newton's method but would set the transient term to zero,

$$ 0 = F_1(t,\boldsymbol{u,v}) \\ 0 = F_2(t,\boldsymbol{u,v}) \\ $$

Is there a way to do this with MOL solvers (e.g. in MATLAB or Python)?

  • 1
    $\begingroup$ Do you mean $d\boldsymbol{v}/dt$ in the second equation? $\endgroup$ – Bill Barth Sep 25 '13 at 13:12

A timestepping method can be considered always as an iterative solver for the steady state problem if the solution is guaranteed to run into a steady state. There are several aspects here:

  1. If the timestepping scheme is implicit and sufficiently stable, you can increase the time steps when you approach the limit. This is often referred to as pseudo timestepping and can be considered a stabilization of the Newton method you use to solve the nonlinear system. This might be necessary if several steady states exist and you want to stay close to the initial value.

  2. If the timestepping scheme is explicit, it amounts to a defect correction scheme of the form \begin{gather} x^{(k+1)} = x^{(k)} - \omega \mathcal F(x^{(k)}), \end{gather} where the damping parameter $\omega$ coresponds to the timestep. This might work nicely, if the problem is not stiff. For stiff problems, the time step will be very small and you will need excruciatingly many iterations.

In both cases, there are more efficient methods, depending on the actual system you solve. In particular, if a Newton method with globalization converges to the desired solution, it will be faster by orders of magnitude.

| cite | improve this answer | |
  • 1
    $\begingroup$ I agree with you 99% of the time about Newton's method being faster, but I've found with some really degenerate PDE (e.g. the Stefan problem) that Newton fails and pseudo-timestepping works. $\endgroup$ – Daniel Shapero Sep 25 '13 at 17:44
  • $\begingroup$ The typical answer in such cases is to run a pseudo-timestepping procedure until you get close enough to the steady state for Newton's method to start working. Then switch to Newton's method and get to the solution in just a few iterations. $\endgroup$ – Wolfgang Bangerth Sep 26 '13 at 1:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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