Solving nonlinear PDE in Python with LSODA

Question

I am attempting to solve a nonlinear diffusion equation of the form $\partial_t u = \partial_x (\kappa(u) \partial_x u)$, where the conductivity function $\kappa(u)$ is a power law $\kappa = u^{5/2}$, using the LSODA time integrator in Python interfaced through SciPy. The equation is discretized by finite difference on a spatial grid, and the resulting system of ODEs is time-integrated by LSODA. The spatial domain is $x \in [0,1]$, the boundary values are fixed as $u(0,t)=u_l$ and $u(1,t)=u_{r}$. Here is a little Python code implementing all this.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint


def kappa(z):
    #-heat conductivity
    return z**2.5
    

def rhs(u, t, dx, param):
    N = len(u) - 1
    rhs = np.zeros(N+1)

    #-fixed value on the left end
    rhs[0] = 0.0
    
    for i in range(1, N):
        rhs[i] = (1/dx**2)*(
            (kappa(u[i])*(u[i+1] - 2*u[i] + u[i-1]) +
             (kappa(u[i+1])-kappa(u[i-1]))*(u[i+1]-u[i-1])/4.0)
            )

    #-fixed value on the right end
    rhs[N] = 0.0    
        
    return rhs


def run_hc1d(L = 1.0, nx = 40, lval=1.0, rval=3.0, nt=10, tmax=1e0):

    #-spatial grid
    x = np.linspace(0, L, nx+1)
    dx = x[1] - x[0]

    
    U_0 = np.zeros(nx+1)

    #-boundary values
    U_0[0] = lval
    U_0[nx] = rval

    #-initial values
    U_0[1:nx] = np.min([lval,rval])
    
    #-time grid
    t = np.linspace(0,tmax,nt)

    #-solving ODEs
    u = odeint(rhs, U_0, t, args=(dx,0.))

    plt.plot(x,u[0,:])
    plt.title("Nonlinear diffusion equation w/ LSODE")
    plt.xlabel("x")
    for it in range(0,nt):
        plt.plot(x,u[it,:])
    plt.show()


run_hc1d(tmax=0.1)

Using $u_l$=1 and $u_r$=3 I can obtain the solution with my Python code, calling it as

>>>run_hc1d(tmax=0.1, rval=3.0)

And it produces a reasonably looking solution, here are several snapshots:

However, using a slightly larger value for $u_r$ results in failure of the time integrator, e.g., it fails if the code is called as

 >>> run_hc1d(tmax=0.1, rval=4.0)

This behavior is puzzling since this equation does not look hard to solve, after all. Why is the time integrator failing? What can be done to make it better behaving, in terms of the choices for the solver parameters, finite difference etc?

Answer 1

I've run your code, and indeed there seems to be a strange issue whit for instace rval=4.. A very suspicious observation is that if your urn the code multiple times, the solution profile differ, i.e. the divergence does not always occur at the same time (or did so every 3-4 runs)...

The odeint package is old and has been replaced since a few years by the solve_ivp one, which is more user-friendly as well. It also has an interface to Fortran LSODA ODE solver.

I've changed your code to use this library instead (and also vectorised the rhsfunction). It now works flawlessly, with all the solvers available in the solve_ivppackage. Not sure what was the origin of your issue with LSODA, but maybe something is not handled well in the odeint library.

EDIT: I have added the comparison with the odepackage as well, which is an old cousin of odeint. It works slightly better, however it does not yield the correct solution it seems, even with tight error tolerances. I would stick with solve_ivp for safety. But that behaviour is really astonishing !

EDIT2: just thinking, are you sure your problem is well-posed ? If there are multiple solutions, this could explain the different solutions obtained with ode and solve_ivp.

Here is my modification of the code:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint, solve_ivp, ode

def kappa(z):
    #-heat conductivity
    return z**2.5

def rhs(u, t, dx, param):
    rhs = np.zeros_like(u)
    #-fixed value on both ends
    rhs[0] = 0.0
    rhs[-1] = 0.0    
    
    rhs[1:-1] = (1/dx**2)*(
        (kappa(u[1:-1])*(u[2:] - 2*u[1:-1] + u[:-2]) +
         (kappa(u[2:])-kappa(u[:-2]))*(u[2:]-u[:-2])/4.0)
        )
        
    return rhs

# default values
L = 1.0; nx = 40; lval=1.0; rval=3.0; nt=10; tmax=1e0

# test case
tmax=0.1; rval=4.

if __name__=='__main__':

    #-spatial grid
    x = np.linspace(0, L, nx+1)
    dx = x[1] - x[0]

    #-initial solution
    U_0 = np.zeros(nx+1)
    U_0[0] = lval
    U_0[nx] = rval
    U_0[1:nx] = np.min([lval,rval])
    
    #-time grid
    t = np.linspace(0,tmax,nt)

    #-solving ODEs
    nChoice=3
    if nChoice==1:
      u = odeint(rhs, U_0, t, args=(dx,0.), atol=1e-9, rtol=1e-9)
    elif nChoice==2:
      sol = solve_ivp(fun=lambda t,x,*args: rhs(x,t,dx,0.), t_span=[t[0],t[-1]], y0=U_0,
                      t_eval=t, method='DOP853', first_step=1e-10, # to avoid warnings with explicit methods
                      atol=1e-9, rtol=1e-9)
      u = sol.y.T
    else:
      r = ode(lambda t,x,*args: rhs(x,t,dx,0.)).set_integrator('vode', method='adams',
                                                             atol=1e-9, rtol=1e-9)
      r.set_initial_value(U_0, t[0])
      u=[U_0]
      for i in range(len(t)-1):
        u.append( r.integrate(t[i+1]) )
      u=np.array(u)
    
    plt.plot(x,u[0,:])
    plt.title("Nonlinear diffusion equation")
    plt.xlabel("x")
    for it in range(0,nt):
        plt.plot(x,u[it,:], label='t={:.2e}'.format(t[it]))
    plt.legend()
```