I have trouble finding guidance on choosing the timestep for solving non-linear ODEs with an external time-varying input.
I have a non-linear ODE which I need to solve: $$ \dot y = f( y(t),u(t) ) $$ $u(t)$ is a known external forcing function. (FYI my application is inelastic deformation and $u$ is an external applied displacement)
If I use a numerical method such as $$ y(t+dt) = y(t)+dt\cdot f(y(t+\theta), u(t+\theta)) $$ with $\theta$ some value in the interval $0<\theta<dt$, then it’s intuitively clear to me (and my numerical experiments agree) that if $u(t)$ changes “a lot” on the interval $0<\theta<dt$ then the numerical solution diverges. I need some guidance on choosing the time step $dt$ based on $u(t)$. I tried some things like "$u(t)$ changes less than p%" but it doesn't work consistently, so I need some theoretical guidance on how to choose dt.
Does someone have a reference to a book chapter on that? I tried checking the library on numerical methods for ODEs or state space methods but I couldn’t find anything, I could only find information on ODEs without the forcing function $u(t)$ but I’m sure I’m not the first one trying to solve numerically a non-linear state space system. Does someone have any reference for me?
Every non-autonomous ODE can be transformed into an autonomous system adding a new unknown. In your case, $y'(t)=f(y(t),u(t))$, can be rewritten as $y'(t)=f(y(t),z(t))$, and $z'(t)=u'(t)$. Every textbook on numerical analysis discusses the stability of your numerical method in this case. Obviously, first, you have to specify the meaning of $y(t+\theta)$ in terms of $y(t+\Delta t)$ and $y(t)$; for example, $y(t+\theta \Delta t)\approx \theta y(t+\Delta t)+(1-\theta) y(t)$.
The stability condition for the numerical method of an autonomous ODE system is given in terms of $\lambda \Delta t$, where $\lambda$ is a complex number representing the eigenvalues of the Jacobian of the nonlinearity (in your case, $\lambda$ is the largest of $\partial f/\partial y$, and $d u/d t$; note that $\partial f/\partial u$ does not appear in the eigenvalues). The stability condition requires that $|1+\theta \lambda \Delta t|<1$, hence $\Delta t$ must be small enough to cope with both $\partial f/\partial y$ and $d u/d t$.
See any textbook in numerical analysis for further details in the stability analysis of your specific numerical method.