# Solving nonlinear PDE in Python with LSODA

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?

• From browsing docs it looks like LSODA will use a finite-difference approximation to the Jacobian of the ODE's right-hand side. You could try supplying an analytical Jacobian to see if that helps. Also, your space discretization doesn't look conservative, but I'm not sure if that's relevant. Jul 15, 2022 at 3:12
• @DanielShapero - Thanks, that's a good thought, but I tried it with the analytic Jacobian and the issue is still there. Jul 15, 2022 at 15:11
• Could the initial condition be a problem here? Maybe some version of the LSODA wrapper sets default parameters that don't handle the jump at the boundary gracefully. Jul 17, 2022 at 14:18
• @Endulum Possibly so; but we'll probably need to look at the wrapper source code to figure it out. Jul 17, 2022 at 14:28

## 1 Answer

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()
$$$$
`
• Wow, that's amazing, thanks for looking into this! I tried it with solve_ivp, using LSODA, and with exactly the same options it is performing much better than calling LSODA through odeint. So, indeed, using different interfaces to the same (?) Fortran ODE package gives rise to quite different results in this case. I am inclined to declare victory on this for now, and pursue it further with the SciPY project developers. Jul 15, 2022 at 15:09
• The analytic solution is unique, at least for the steady state it is easy to see that because it boils down to $\partial_{xx} u^{7/2}$=0. Jul 15, 2022 at 15:12