Operator exponential

[mjyshin]

First of all, thank you for developing this package, it's a very instructive way to learn about function decomposition methods!

And I do like the way time dependent PDEs can be solved by extending the state-space to include the time domain (see this example) and defining the IC as a Dirichlet BC.

1) But I was wondering if there was any discussion/plan of solving linear PDEs (u' = A*u) with the exponential of operators (u = exp(A*t)*u0) by overloading exp? Matlab's chebfun has the expm function (see examples here and here), but it seems like the operators there have BCs baked in. There is an ApproxFun example too, but it's not very intuitive as one is not applying the propagator to a Fun type directly, and you need to explicitly specify how many coefficients to use.

2) Also for more general PDEs, what is the status on using Fun in ODEProblem? I noticed that naively defining something like

function model!(du,u,p,t)
    du = F(u,p)
end

and plugging it into ODEProblem(model!,u0,tspan,p) won't work. I noticed this question has been asked before, but since it's been a few years, I was wondering if there was any progress on this front?

Thank you again!

[dlfivefifty]

Operator exponentials are hard, particularly with boundary conditions. You can define them with Cauchy integral formula or with rational approximations, but to my knowledge this requires knowledge of the spectrum.

[mjyshin]

I tried using your method in your heat equation example with both BCs Neumann, but I get a LAPACKException(1) error: See the last cell of this. Am I doing something wrong?

Also what exactly is going on in these lines?

B = [ldirichlet(S); rneumann(S)]
    B = B[1:2,1:n]
    B = B[1:2,1:2]\B

M = C[1:n-2,1:n]
    M = M - M[:,1:2]*B # remove degrees of freedom
    A = M[:,3:end]\L[1:n-2,3:n]

For nonlinear problems and interfacing with DifferentialEquations, would @ChrisRackauckas be more familiar with the progress on that front?

[dlfivefifty]

It’s because Neumann conditions have a zero first column. To make it work you’d have to permute the first 3 columns first

[mjyshin]

I am having some trouble getting that workaround to work. I've done

perm = [2;3;1;4:n]
B = B[1:2,1:n][:,perm]

but what other columns/vectors should I permute and permute back? And is B = B[1:2,1:2]\B (making the first sublock identity) necessary?

[dlfivefifty]

The method works as follows. First we make sure B[1:2,1:2] == I. Then we want to solve a time-evolution equation with boundary conditions:

B*u = 0
M*u_t = A*u

Differentiating we get B*u_t = 0 which implies M[:,1:2]*B*u_t = 0 and A[:,1:2]*B*u = 0. Thus we also have:

(M - M[:,1:2]*B)*u_t = (A - A[:,1:2]*B)*u

The new mass matrix and RHS have zeros in the first two columns so we can drop those rows.

But if B has a zero first column we need to incorporate a permutation. That is, we let v = P'*u so that we have

B[1:2,1:2] \ B*P*v = 0
M*P*v_t = A*P*v

Relabelling, this is the same as before hence one can proceed.

I'll try to get an example working in ContinuumArrays.jl which will be cleaner

[mjyshin]

Thank you, this is fantastic! I updated my little heat equation tutorial with this implementation, along with some basic checks. I'm new to this having only recently seen your JuliaCon 2018 talk, so not everything is clear, but I'm slowly working through your 2013 SIAM paper and learning a lot from it.

I tried the same PDE with the extended state space (X x T) formulation here, but I noticed that on my (fairly old) laptop, the backslash solve gets stuck even at tolerance=1e-2, so this method seems to be preferable in this case.

[dlfivefifty]

I think in this setting time stepping is preferable: solutions to PDEs in rectangles usually have corner singularities which a tensor based method will try to resolve. time-stepping avoids this because these singularities aren’t present on slices.

Note the reason your simple example does not hit this issue is not because of the number of coefficients, but because (I believe) every even derivative vanishes, which is precisely the condition needed to avoid corner singularities.

[MikaelSlevinsky]

Hey @mjyshin, you can also do this with a symmetric-definite and banded eigendecomposition when propagating with self-adjoint linear differential operators. A reference on the code below is here, and this was followed up with more detail by my former M Sc student in his thesis. Sample code for your problem:

using ApproxFun, LinearAlgebra
import ApproxFun: bandwidth

d = Segment(0..1)
S = Ultraspherical(0.5, d)
NS = NormalizedPolynomialSpace(S)
D = Derivative(S)
θ = (D=0.05, )
Δ = D^2
L = -θ.D*Δ # We'll negate this later
ϵ = 10^-1.5
u0 = Fun(x->(tanh((x-0.3)/ϵ)-tanh((x-0.7)/ϵ))/2, S)

C = Conversion(domainspace(L), rangespace(L))
B = [ldirichlet(S); rneumann(S)]

QS = QuotientSpace(B)
Q = Conversion(QS, S)
D1 = Conversion(S, NS)
D2 = Conversion(NS, S)
R = D1*Q
P = cache(PartialInverseOperator(C, (0, bandwidth(L, 1) + bandwidth(R, 1) + bandwidth(C, 2))))
AA = R'D1*P*L*D2*R
BB = R'R

n = ncoefficients(u0)
SB = Symmetric(BB[1:n,1:n], :L)
SA = Symmetric(AA[1:n,1:n], :L) + 0*SB # + 0*SB is a bug workaround for bandwidths in A smaller than bandwidths in B.

F = eigen(SA, SB);
F.values .= - F.values # negation complete!

w = Fun(x-> Number((B*u0)[1]) + Number((B*u0)[2])*x, S, 2) # A degree-1 polynomial that exactly captures nonzero boundary data

v0 = Q\(u0-w) # We'll represent u(x,t) = v(x,t) + w(x), where v satisfies 0 boundary data for all time
pad!(v0, n)

# u' = L u with B u = c
# Let u = v + w, where v evoles and w contains the nonzero boundary data.
# Then v' = L (v + w) = Lv since Lw = 0.
# L = B^{-1} A = VΛV^{-1},  VᵀBV = I => V^{-1} = VᵀB
# L = VΛVᵀB
# exp(tL) = Vexp(tΛ)VᵀB

vt = t -> Fun(QS, F.vectors*(Diagonal(exp.(t.*F.values))*(F.vectors'*(SB*v0.coefficients))))

t = 0:0.5:5
v = vt.(t)

u0(0) - v[end](0) - w(0) # left boundary condition
u0'(1) - (Q*v[end])'(1) - w'(1) # right boundary condition

[mjyshin]

Wow, thank you both for your inputs. @MikaelSlevinsky I appreciate all the resources. I certainly have a long reading list for the next several weekends to parse through your code. A few thoughts: 1) S = Ultraspherical(0.5, d) seems to be required (e.g., using Chebychev(d), Ultraspherical(1, d), etc. would not conserve the "heat"). 2) For some reason B = [ldirichlet(S); rneumann(S)] runs faster than [lneumann(S); rneumann(S)] and much faster than [ldirichlet(S); rdirichlet(S)]. Maybe it's my machine... 3) Question: Do you have a few cases in mind where this method would be preferred over the one @dlfivefifty outlined above (or expm in Chebfun, if you're aware of how that's implemented)?

Thanks again! I think it'll be nice to someday have this type of propagator operation baked into ApproxFun for everyone.

[MikaelSlevinsky]

1. Yes, Legendre series are required for this method. Not really a big deal since Chebyshev and Legendre series can be converted to and from quite rapidly.
2. I had hardcoded the left-dirichlet right-neumann boundary condition into the w(x). So the generic version would be

   w = Fun(S, B[1:2,1:2]\[Number((B*u0)[1]); Number((B*u0)[2])]) # A degree-1 polynomial that exactly captures nonzero boundary data

   The general boundary condition check is B*u0 - B*(Q*v[end] + w) # check boundary conditions. Probably the code would otherwise hang (or work slowly) on v0 = Q\(u0-w). That line works unless, of course, you hit Neumann conditions. Then, as with @dlfivefifty's method, a degree-1 polynomial can't have two different slopes, so you need to modify two lines now: firstly, the w(x), which uses a Moore-Penrose pseudo inverse with three columns of the boundary conditions to get minimum 2-norm Legendre coefficients, is now

   w = Fun(S, B[1:2,1:3]\[Number((B*u0)[1]); Number((B*u0)[2])]) # A degree-2 polynomial that exactly captures nonzero boundary data

   Secondly, since Lw ≠ 0, the general solution is now

   u = exp(tL) v0 + exp(tL) w = exp(tL) v0 + w + tLw, because (tL)^k w = 0 for k ≥ 2, since L is diffusion and w is at most quadratic.

   Now to check the boundary conditions in the Neumann case B*u0 - B*(Q*v[end] + w + t[end]*L*w) # check Neumann boundary conditions.

3. As for which method is best, I don't know. Both methods have subtleties which have probably prevented a generic exp implementation. This method gives you a structured discretization of a self-adjoint linear differential operator which allows finite sections of exp(tL) to be computed in O(n^2) flops. It looks like your implementation of @dlfivefifty's method gives either a dense regular eigenvalue problem or an almost-banded generalized eigenvalue problem for A = M[:,3:n]\L[1:n-2,3:n] # BC-merged generator. I think the complexity of exp(A*t) is O(n^3) flops in either case. Does n get really large in the evolution of the heat equation? Normally, no. So maybe it's an engineering problem which method is better. What if L is not self-adjoint? If it's skew-adjoint then you can recover that structure (with purely imaginary finite-section spectrum). If it's neither nor, then there's no point trying to get more structure out of your solver.

[mjyshin]

@MikaelSlevinsky, your method works as advertised for Dirichlet and Neumann (and mixed) BCs. If we define L using differential operators on, say Chebfun(...), is there a way to convert it to a mapping from Ultraspherical(0.5,...) to the appropriate range?

Also, do you know of a way to propagate u0 forward to u given a nonlinear operator (u' = N(u))?

[MikaelSlevinsky]

If we define L using differential operators on, say Chebfun(...), is there a way to convert it to a mapping from Ultraspherical(0.5,...) to the appropriate range?

This can be done with connection coefficients, but since these are generally upper-triangular and dense, the bandwidths of banded operators are hard to preserve (they would require some sophisticated numerical test for small coefficients worth truncating). In that sense, it's best to work with the original differential operator. The operator's polynomial variable coefficients don't need to be expanded in Legendre series; they can be in Chebyshev expansions in case that helps.

Also, do you know of a way to propagate u0 forward to u given a nonlinear operator (u' = N(u))?

I think this is just general autonomous time-stepping, so I don't know anything specific to orthogonal polynomials in this regard.

[mjyshin]

This is a somewhat unrelated question, but is there any meaning to the adjoint of an ApproxFun operator in Julia L' in relation to the formal adjoint L* of an operator? I'm guessing not, since L would be an operator from, say, Ultraspherical(0.5,d) to Ultraspherical(2.5,d) for L = Δ on Ultraspherical(0.5,d), while L' would be from Ultraspherical(2.5,d) to Ultraspherical(0.5,d) and functions like Conversion(domainspace(L'),rangespace(L')) in either implementations of exp above would not work. Just wanted to know your thoughts on this.

[dlfivefifty]

This is a motivation for using ContinuumArrays.jl as it makes adjoints well-defined

[mjyshin]

Wow, very interesting (ContinuumArrays I mean). Looks kind of Chebfun-y, given that a function f has size(f) that isn't () like in ApproxFun and is accessed with an index. Are you planning for it to be integrated into ApproxFun or for it to replace it?

I've been using ApproxFun on the Kolmogorov equations (here), and it would be cool to be able to just define the generator L for the backward equation and get the adjoint generator for the forward equation with L'.