Can the method of lines technique be used to solve a ODE directly for the steady-state value?

Question

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)?

Answer 1

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.