How can I solve steady-state problems?

[Rufnacous]

Hi, I am looking at moving an existing computational model to a PyElastica implementation.
An important consideration in my problem is finding steady states. That is, x''=x'=0.

Imagining a Forward-Euler integration, one might generate a normal solution like so:
x' <- dt * f(x,x',t)
x <- dt * x'
Along a similar vein, I have previously used this method to converge on steady state solutions:
x' <- 0
x <- dt * f(x,x',t)
This might fly in the face of traditional iterative methods (I tried several "regular" ones before landing on this) but it works well for my problem.

Of course, it requires quite unorthodox changes to the integrator code. I've just tried changing PyElastica to achieve this, but have gotten lost in the class definitions.

Can anyone point me to where in PyElastica's integrators would be a good place for trying these changes?
If not, does anyone have alternate suggestions for steady state solution with PyElastica?

Thanks :)
Rufus

[armantekinalp]

Hi @Rufnacous

In PyElastica we are using Position Verlet time integrator and if you want to see/hack in here is the time-stepper. do step method does one Position verlet time step. You can hack in and change the kinematic and dynamic stepper parts.

But I think there is a fundamental problem. I don't think you can compute steady-state by just setting up velocities to zero, because then positions are not going to be updated at anytime. In order to update positions without updating velocities you will need an iterative solver.

Let me know if this helps.

[Rufnacous]

Thanks @armantekinalp I'll look at the do step method today and will try to remember to report back if I have success.

RE your second point, I can see why you would say this - physically it doesn't make any sense to expect position updates with velocity at zero. This misunderstanding is just from my explanation being a little lacking. In the iterative method I am describing, the position vector is updated not by dt * velocity, but rather dt * acceleration. Of course, this is not correct physically or even dimensionally.

There are probably lots of reasons why this method shouldn't be used, but my working theory is that for continuous and conservative force-vector fields, the position iteration will always move the system towards static equilibrium. Potentially this method is problematic if dx''/dx =/= 0 at equilibrium, but for my application this is not the case.

[Rufnacous]

So I've written my integrator:

`
class IterativeEquilibirum:
Tag = SymplecticStepperTag()
def _first_prefactor(self, dt):
return 0.5*dt
def _first_kinematic_step(self, System, time: np.float64, dt: np.float64):
pass

def _first_dynamic_step(self, System, time: np.float64, dt: np.float64):
    rates = 0 * System.dynamic_states.rate_collection

    overload_operator_dynamic_numba( 
        rates,
        System.dynamic_rates(time, dt),
    )

    prefac = self._first_prefactor(dt)
    
    overload_operator_kinematic_numba(
        System.n_nodes,
        prefac,
        System.kinematic_states.position_collection,
        System.kinematic_states.director_collection,
        np.reshape(rates[0,:],(3,-1)),
        np.reshape(rates[1,:-3],(3,-1)),
    )

`

Overall it's very hacky, and is limited incredibly by the stiffness of PyElastica's equations compared to the stiffness of my previous implementation (which did not include shear deformations) so ultimately I am not sure I can use PyElastica for my particular problem.

Thanks for the pointer @armantekinalp - I will mark this as closed if I can work out how.