Access to (Nonlinear) Differential Operator

[gnorman7]

Is there a reasonably straight-forward way to access the differential operator? My end goal is to look at the eigenvalues. To clarify further, for the semi-discrete system, with operator $A$:

$$ u_t = A(u) u $$

which in an explicit finite difference time step scheme gives, for time step $i$:

$$ u^{i+1} = L^i u^i $$

I am hoping to find a way to access $L^{i}$, ideally for all $i$, but this may be a little too in-the-weeds. As an example, this approach for 1D advection would give $L^{i}$ (constant for all $i$ as the operator is linear) such that all its eigenvalues are strictly imaginary. Similarly, 1D diffusion would give negative real eigenvalues. For more explanation, see this set of notes, although disregarding the spectra for the time steppers.

For reference, here is the code setup that I hope to modify. I also want to do this for another application (see #355 ) with a similar structure (just changing the StringPDE class). As long as this structure stays similar, performance is not a concern.

class StringPDE(PDEBase):
    def __init__(self, N_string):
        self.N_string = N_string

    def evolution_rate(self, state, t=0):
        du_dx = state.gradient(bc='periodic').to_scalar()
        d2u_dx2 = state.laplace(bc='periodic')
        du_dt = eval(self.N_string)

        return du_dt

# Viscous Burgers Equation: u_t + u*ux = nu*u_xx, nu=0.1
N = "-state*du_dx+0.1*d2u_dx2"

dt = 0.001
dx = dt/0.03

n_x = int(16/dx)

# Boundary and Initial Conditions
grid = CartesianGrid([[-8.0, 8.0]], [n_x], periodic=True)
state = ScalarField.from_expression(grid, "-1.0*sin(pi*x/8)")
eq = StringPDE(N)

# Need a tracker, as we want the time history, not just the steady state.
storage = MemoryStorage()

# Default is explicit time stepper, which we probably need.
eq.solve(state, t_range=10.0, tracker=storage.tracker(0.1), dt=dt)

[david-zwicker]

I'm afraid this is not really possible with py-pde since this is not how we represent PDEs internally. Instead, we directly evaluate the operator O(u) describing the right-hand side (which you write as A(u)u). For general non-linear PDEs, the separation is not too meaningful (unless you perform a linear stability analysis, where A(u) is then the Jacobian).

[gnorman7]

Okay, that makes sense. To clarify, I ultimately want to perform a linear stability analysis for my O(u). It sounds like this analysis is really just dependent on the operator O, and the linearization state. Therefore this analysis should not be dependent on the solution scheme.

[david-zwicker]

True, but I guess it requires determining the Jacobian, which we currently don't do. I don't think py-pde provides efficient mechanisms to calculate what you require.