title: Regularizing viscosity with scipy's ode solvers
body: <p><strong>Consider</strong> for the sake of simplicity <strong>the following equation</strong> (Burgers equation):</p>
  
  <p><a href="https://i.stack.imgur.com/5TyM4.gif" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/5TyM4.gif" alt="Burgers equation"></a></p>
  
  <p>Let's solve it using scipy (in my case <code>scipy.integrate.ode.set_integrator("zvode", ..).integrate(T)</code>) with a variable time-step solver.</p>
  
  <p>The issue is the following: <strong>if we use the naïve implementation</strong> in Fourier space</p>
  
  <p><a href="https://i.stack.imgur.com/SjoSm.gif" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/SjoSm.gif" alt="enter image description here"></a></p>
  
  <p>then the viscosity term <code>nu * d2x(u[t])</code> can cause an overshoot if the time step is too big. This can lead to a fair amount of noise in the solutions, or even to (fake) <strong>diverging solutions</strong> (even with stiff solvers, on slightly more complex version of this equation).</p>
  
  <p>One way <strong>to regularize this</strong> is to evaluate the viscosity term at step <code>t+dt</code>, and the update step becomes</p>
  
  <p><a href="https://i.stack.imgur.com/EA34t.gif" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/EA34t.gif" alt="enter image description here"></a></p>
  
  <p><strong>This solution works well when programmed explicitly. How can I use scipy's variable-step ode solver to implement it ?</strong> To my surprise I haven't found any documentation on this fairly elementary thorny issue...</p>