Regularizing viscosity with scipy's ode solvers

Question

Consider for the sake of simplicity the following equation (Burgers equation):

Let's solve it using scipy (in my case scipy.integrate.ode.set_integrator("zvode", ..).integrate(T)) with a variable time-step solver.

The issue is the following: if we use the naïve implementation in Fourier space

then the viscosity term nu * d2x(u[t]) 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) diverging solutions (even with stiff solvers, on slightly more complex version of this equation).

One way to regularize this is to evaluate the viscosity term at step t+dt, and the update step becomes

This solution works well when programmed explicitly. How can I use scipy's variable-step ode solver to implement it ? To my surprise I haven't found any documentation on this fairly elementary thorny issue...

Answer 1

You actually can't, or on the other extreme, odeint or ode->zvode already does that to any given problem.

To the first, you would need to give the two parts of the equation separately. Obviously, that is not part of the solver interface. Look at DDE and SDE solvers where such a partition of the equation is actually required.

To the second, odeint and ode->zvode use implicit multi-step methods, which means that the values of u(t+dt) and the right side there enter the computation and the underlying local approximation.

You could still try to hack your original approach into the solver by providing a Jacobian function that only contains the second derivative term, but quite probably you will not achieve an improvement.

You could operator-partition the ODE and solve the linear part separately introducing

vhat(k,t) = exp(nu*k^2*t)*uhat(k,t)

so that

d/dt vhat(k,t) = -i*k*exp(nu*k^2*t)*conv(uhat(.,t),uhat(.,t))(k)