3.11.2. Stationary solutions
By simple time integration, we have seen how the transcritical normal form avoids unstable stationary solutions and tries to approach stable stationary solutions. While time integration can help to approach a stable stationary solution, one is unable to approach unstable stationary solutions that way. Also, it takes a lot of time steps to approach the stable ones. Instead, one can just jump to a stationary solution, both stable and unstable ones, by solving the stationary problem:
Of course, one could implement this equation by hand again, but it is better to reuse the full equation with \(\partial_t x\)-term, i.e. (3.13), and just tell pyoomph to search for a stationary solution instead of integrating the temporal evolution. To that end, we first import the transcritical bifurcation equation and problem from the previous example bifurcation_transient_transcritical.py. Afterwards, we use the solve() method of the problem, which exactly looks for a stationary solution.:
from bifurcation_transient_transcritical import * #Import the equation and problem from the previous script
if __name__=="__main__":
with TranscriticalProblem() as problem:
problem.quiet() # Do not pollute our output with all the messages
problem.r.value=1 # Parameter value
problem.x0=0.001 # Start slightly above the unstable solution x=0
ode = problem.get_ode("transcritical") # Get access the ODE (note: it will initialize the problem!)
xvalue = ode.get_value("x") #Get the current value of x
print(f"We are starting at x={xvalue}")
problem.solve(timestep=None) # Solve without a timestep means stationary solve
xvalue = ode.get_value("x") # Get the current value of x
print(f"Currently, we are at the stationary solution x={xvalue} with r={problem.r.value}")
ode.set_value(x=0.8) # Set the current value of x
problem.solve() # we can omit timestep=None, since it is default
xvalue = ode.get_value("x")
print(F"Currently, we are at the stationary solution x={xvalue} with r={problem.r.value}")
First, we tell the problem class to be quiet(), so that there is no output from the code generation and solution. Thereby, the print statements can be better found in the output. We can get access to the ODE by using the get_ode() method. Note that this will initialize the problem, i.e. generate the code and so on, since before that the problem does not know about the existence of an ODE domain called "transcritical". This domain is added in define_problem(). With the reference to the ODE, we can use set_value() and get_value() to access the current value of \(x\). Initially, we are close to the unstable solution \(x=0\), so the solve() statement will likely run into this solution, which is indeed the case. In the second solve, we are close to the stable solution at \(x=r=1\), so we run into this with the second invocation of solve().
Obviously, one can easily find a stationary solution by the solve() function. However, this only works if one is close to that solution. If one tries to e.g. start at \(x=0.5\), i.e. by ode.set_value(x=0.5) followed by a problem.solve() statement, it will not find any solution. The reason lies in the internal Newton method, which requires the Jacobian, which is singular at \(x=0.5\). In fact, we are exacly between both solutions, so the system cannot decide which way to go. In that case, one can solve a few steps by temporal integration, either via the run() command or via solve(timestep=...) to move a bit away from the singularity. After that, the normal solve() command will converge to a stationary solution again. Note that the numerical stationary solutions are not exactly \(0\) and \(1\). Deviations come from the fact that (i) we have a nonlinear problem, i.e. the solution procedure is stopped after reaching sufficient accuracy and (ii) numerical calculations on computers always have a limited accuracy due to the finite accuracy of floating point numbers. The first issue can be improved by passing newton_solver_tolerance=1e-20 or even smaller numbers to the solve() command. However, in larger systems the latter, i.e. the finite accuracy of float arithmetic, might hamper to reach tiny accuracy thresholds.