Hostname: page-component-788cddb947-jbkpb Total loading time: 0 Render date: 2024-10-11T01:46:32.675Z Has data issue: false hasContentIssue false

Linear PDE with constant coefficients

Published online by Cambridge University Press:  10 November 2021

Rida Ait El Manssour
Affiliation:
MPI-MiS, Leipzig, Germany
Marc Härkönen
Affiliation:
Georgia Institute of Technology, Atlanta, GA, USA
Bernd Sturmfels
Affiliation:
MPI-MiS, Leipzig, Germany, and UC Berkeley, Berkeley, CA, USA E-mail: bernd@mis.mpg.de
Rights & Permissions [Opens in a new window]

Abstract

We discuss practical methods for computing the space of solutions to an arbitrary homogeneous linear system of partial differential equations with constant coefficients. These rest on the Fundamental Principle of Ehrenpreis–Palamodov from the 1960s. We develop this further using recent advances in computational commutative algebra.

Type
Research Article
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited.
Copyright
© The Author(s), 2021. Published by Cambridge University Press on behalf of Glasgow Mathematical Journal Trust

1. Introduction

Our calculus class taught us how to solve ordinary differential equations (ODE) of the form

(1.1) \begin{equation} c_0 \phi + c_1 \phi ' + c_2 \phi '' + \cdots + c_m \phi^{(m)} = 0 .\end{equation}

Here we seek functions $\phi = \phi(z)$ in one unknown z. The ODE is linear of order m, it has constant coefficients $c_i \in \mathbb{C}$ , and it is homogeneous, meaning that the right-hand side is zero. The set of all solutions is a vector space of dimension m. A basis consists of m functions

(1.2) \begin{equation} \qquad \qquad \phi(z) = z^a \cdot {\textrm{exp}}(u_i z) .\end{equation}

Here $u_i$ is a complex zero with multiplicity larger than $a \in \mathbb{N}$ of the characteristic polynomial

(1.3) \begin{equation} p(x) = c_0 + c_1 x + c_2 x^2 + \cdots + c_m x^m .\end{equation}

Thus solving the ODE (1.1) means finding all the zeros of (1.3) and their multiplicities.

We next turn to a partial differential equation (PDE) for functions $\phi\,:\, \mathbb{R}^2 \rightarrow \mathbb{R}$ that is familiar from the undergraduate curriculum, namely the one-dimensional wave equation

(1.4) \begin{equation} \qquad \phi_{tt}(z,t) = c^2 \phi_{zz}(z,t), \qquad {\textrm{where}}\ c \in \mathbb{R} \backslash \{0\} .\end{equation}

D’Alembert found in 1747 that the general solution is the superposition of traveling waves,

(1.5) \begin{equation}\phi(z,t) = f(z+ct) + g(z-ct) ,\end{equation}

where f and g are twice differentiable functions in one variable. For the special parameter value $c=0$ , the PDE (1.4) becomes $\phi_{tt} = 0$ , and the general solution has still two summands

(1.6) \begin{equation} \phi(z,t) = f(z) + t \cdot h'(z). \end{equation}

We get this from (1.5) by replacing $ g(z-ct)$ with $ \frac{1}{2c}(h(z{+}ct) - h(z{-}ct)) $ and taking the limit $c \rightarrow 0$ . Here, the role of the characteristic polynomial (1.3) is played by the quadratic form

(1.7) \begin{equation} q_c(u,v) \quad = \quad v^2 - c^2 u^2 = (v-cu)(v+cu).\end{equation}

The solutions (1.5) and (1.6) mirror the algebraic geometry of the conic $\{q_c=0\}$ for any $c \in \mathbb{R}$ .

Our third example is a system of three PDE for unknown functions

\begin{equation*} \psi \,:\, \mathbb{R}^4 \rightarrow \mathbb{C}^2 , (x,y,z,w) \mapsto \bigl(\alpha(x,y,z,w),\beta(x,y,z,w) \bigr) . \end{equation*}

Namely, we consider the following linear PDE with constant coefficients:

(1.8) \begin{equation} \alpha_{xx} + \beta_{xy} = \alpha_{yz} + \beta_{zz} = \alpha_{xxz} + \beta_{xyw}= 0.\end{equation}

The general solution to this system has nine summands, labeled $a,b, \ldots,h$ and $(\tilde \alpha,\tilde \beta)$ :

(1.9) \begin{align} \alpha & = a_z(y,z,w)- b_y(x,y)+c(y,w)+xd(y,w)+xg(z,w)-xyh_z(z,w)+ \tilde \alpha(x,y,z,w), \nonumber\\[4pt] \beta & = -a_y(y,z,w)+b_x(x,y)+e(x,w)+zf(x,w) +xh(z,w)+ \tilde \beta(x,y,z,w) .\end{align}

Here, a is any function in three variables, b, c, d, e, f, g, h are functions in two variables, and $\tilde \psi = (\tilde \alpha, \tilde \beta)$ is any function $\mathbb{R}^4 \rightarrow \mathbb{C}^2$ that satisfies the following four linear PDE of first order:

(1.10) \begin{equation} \tilde \alpha_{x} + \tilde \beta_{y} = \tilde \alpha_{y} + \tilde \beta_{z} = \tilde \alpha_{z} - \tilde \alpha_{w}= \tilde \beta_{z} - \tilde \beta_{w} = 0 .\end{equation}

We note that all solutions to (1.10) also satisfy (1.8), and they admit the integral representation

(1.11) \begin{equation}\tilde \alpha = \int\!\! t({\textrm{exp}}(s^2x+sty+t^2(z{+}w))) d\mu(s,t) ,\ \ \tilde \beta =-\!\! \int\!\! s({\textrm{exp}}(s^2x+sty+t^2(z{+}w))) d\mu(s,t),\end{equation}

where $\mu$ is a measure on $\mathbb{C}^2$ . All functions in (1.9) are assumed to be suitably differentiable.

Our aim is to present methods for solving arbitrary systems of homogeneous linear PDE with constant coefficients. The input is a system like (1.1), (1.4), (1.8), or (1.10). We seek to compute the corresponding output (1.2), (1.5), (1.9), or (1.11), respectively. We present techniques that are based on the Fundamental Principle of Ehrenpreis and Palamodov, as discussed in the classical books [Reference Björk7, Reference Ehrenpreis17, Reference Hörmander23, Reference Palamodov32]. We utilize the theory of differential primary decomposition [Reference Cid-Ruiz and Sturmfels12]. While deriving (1.5) from (1.4) is easy by hand, getting from (1.8) to (1.9) requires a computer.

This article is primarily expository. One goal is to explain the findings in [Reference Chen, Cid-Ruiz, Härkönen, Krone and Leykin8Reference Cid-Ruiz and Sturmfels12], such as the differential primary decompositions of minimal size, from the viewpoint of analysis and PDE. In addition to these recent advances, our development rests on a considerable body of earlier work. The articles [Reference Damiano, Sabadini and Struppa15, Reference Oberst29, Reference Oberst31] are especially important. However, there are also some new contributions in the present article, mostly in Sections 4, 5, and 6. We describe the first universally applicable algorithm for computing Noetherian operators.

This presentation is organized as follows. Section 2 explains how linear PDE are represented by polynomial modules. The Fundamental Principle (Theorem 2.2) is illustrated with concrete examples. In Section 3, we examine the support of a module and how it governs exponential solutions (Proposition 3.7) and polynomial solutions (Proposition 3.9). Theorem 3.8 characterizes PDE whose solution space is finite dimensional. Section 4 features the theory of differential primary decomposition [Reference Chen and Cid-Ruiz9, Reference Cid-Ruiz and Sturmfels12]. Theorem 4.4 shows how this theory yields the integral representations promised by Ehrenpreis–Palamodov. This result appeared implicitly in the analysis literature, but the present algebraic form is new. It is the foundation of our algorithm for computing a minimal set of Noetherian multipliers. This is presented in Section 5, along with its implementation in the command solvePDE in Macaulay2 [Reference Grayson and Stillman21].

The concepts of schemes and coherent sheaves are central to modern algebraic geometry. In Section 6, we argue that linear PDE are an excellent tool for understanding these concepts and for computing their behaviors in families. Hilbert schemes and Quot schemes make an appearance along the lines of [Reference Chen and Cid-Ruiz9, Reference Cid-Ruiz, Homs and Sturmfels11]. Section 7 is devoted to directions for further study and research in the subject area of this paper. It also features more examples and applications.

2. PDE and polynomials

Our point of departure is the observation that homogeneous linear partial differential equations with constant coefficients are the same as vectors of polynomials. The entries of the vectors are elements in the polynomial ring $R = K[\partial_{1}, \partial_{2}, \ldots,\partial_{n}]$ , where K is a subfield of the complex numbers $\mathbb{C}$ . In all our examples, we use the field $K = \mathbb{Q}$ of rational numbers. This has the virtue of being amenable to exact symbolic computation, e.g. in Macaulay2 [Reference Grayson and Stillman21].

For instance, in (1.1), we have $n=1$ . Writing $\partial = \dfrac{\partial}{\partial z}$ for the generator of R, our ODE is given by one polynomial $p(\partial) = c_0 + c_1 \partial + \cdots + c_m \partial^m$ , where $c_0,c_1,\ldots,c_m \in K$ . For $n \geq 2$ , we write ${\bf z} = (z_1,\ldots,z_n)$ for the unknowns in the functions we seek, and the partial derivatives that act on these functions are $\partial_i = \partial_{z_i} = \dfrac{\partial}{\partial z_i}$ . With this notation, the wave equation in (1.4) corresponds to the polynomial $q_c(\partial) = \partial_2^2 - c^2 \partial_1^2 = (\partial_2 - c \partial_1)(\partial_2 + c \partial_1)$ with $n=2$ . Finally, the PDE in (1.8) has $n=4$ and is encoded in three polynomial vectors

(2.1) \begin{equation}\begin{pmatrix} \partial_1^2 \\[4pt] \partial_1 \partial_2 \end{pmatrix} , \quad \begin{pmatrix} \partial_2 \partial_3 \\[4pt] \partial_3^2 \end{pmatrix} \quad {\textrm{and}} \quad \begin{pmatrix} \partial_1^2 \partial_3 \\[4pt] \partial_1 \partial_2 \partial_4 \end{pmatrix} .\end{equation}

The system (1.8) corresponds to the submodule of $R^2$ that is generated by these three vectors.

We shall study PDE that describe vector-valued functions from n-space to k-space. To this end, we need to specify a space $\mathcal{F}$ of sufficiently differentiable functions such that $\mathcal{F}^k$ contains our solutions. The scalar-valued functions in $\mathcal{F}$ are either real-valued functions $\psi\,:\,\Omega \rightarrow \mathbb{R}$ or complex-valued functions $\psi\,:\,\Omega \rightarrow \mathbb{C}$ , where $\Omega$ is a suitable subset of $\mathbb{R}^n$ or $\mathbb{C}^n$ . Later we will be more specific about the choice of $\mathcal{F}$ . One requirement is that the space $\mathcal{F}^k$ should contain the exponential functions

(2.2) \begin{equation} q({\bf z}) \cdot {\textrm{exp}}( {\bf u}^t {\bf z}) = q(z_1,\ldots,z_n) \cdot {\textrm{exp}}( u_1 z_1 + \cdots + u_n z_n).\end{equation}

Here ${\bf u} \in \mathbb{C}^n$ and q is any vector of length k whose entries are polynomials in n unknowns.

Remark 2.1 ( $k=1$ ) A differential operator $p(\partial) $ in R annihilates the function ${\textrm{exp}}({\bf u}^t {\bf z})$ if and only if $p({\bf u}) = 0$ . This is the content of [Reference Michałek and Sturmfels27, Lemma 3.25]. See also Lemma 3.26. If $p(\partial)$ annihilates a function $ q({\bf z}) \cdot {\textrm{exp}}( {\bf u}^t {\bf z})$ , where q is a polynomial of positive degree, then u is a point of higher multiplicity on the hypersurface $\{ p = 0 \}$ . In the case $n=1$ , when p is the characteristic polynomial (1.3), we have a solution basis of exponential functions (1.2).

Another requirement for the space $\mathcal{F}$ is that it is closed under differentiation. In other words, if $\phi = \phi(z_1,\ldots,z_n)$ lies in $\mathcal{F}$ , then so does $\partial_i \bullet \phi = \dfrac{\partial \phi}{\partial z_i}$ for $i=1,2,\ldots,n$ . The elements of $\mathcal{F}^k$ are vector-valued functions $\psi = \psi({\bf z})$ . Their coordinates $\psi_i$ are scalar-valued functions in $\mathcal{F}$ . All in all, $\mathcal{F}$ should be large, in the sense that it furnishes enough solutions. Formulated algebraically, we want $\mathcal{F}$ to be an injective R-module [Reference Lomadze25]. A more precise desideratum, formulated by Oberst [Reference Oberst28Reference Oberst30], is that $\mathcal{F}$ should be an injective cogenerator.

Examples of injective cogenerators include the ring $\mathbb{C}[[z_1,\ldots,z_n]]$ of formal power series, the space $C^\infty(\mathbb{R}^n)$ of smooth complex-valued functions over $\mathbb{R}^n$ , or more generally, the space $\mathcal{D}'(\mathbb{R}^n)$ of complex-valued distributions on $\mathbb{R}^n$ . If $\Omega $ is any open convex domain in $\mathbb{R}^n$ , then we can also take $\mathcal{F}$ to be $C^\infty(\Omega)$ or $\mathcal{D}'(\Omega)$ . In this paper, we focus on algebraic methods. Analytic difficulties are mostly swept under the rug.

Our PDE are elements in the free R-module $R^k$ , that is, they are column vectors of length k whose entries are polynomials in $\partial = (\partial_1,\ldots,\partial_n)$ . Such a vector acts on $\mathcal{F}^k$ by coordinate-wise application of the differential operator and then adding up the results in $\mathcal{F}$ . In this manner, each element in $R^k$ defines an R-linear map $\mathcal{F}^k \rightarrow \mathcal{F}$ . For instance, the third vector in (2.1) is an element in $R^2$ that acts on functions $\psi \,:\, \mathbb{R}^4 \rightarrow \mathbb{C}^2$ in $\mathcal{F}^2$ as follows:

(2.3) \begin{equation}\begin{pmatrix} \partial_1^2 \partial_3 \\[4pt] \partial_1 \partial_2 \partial_4 \end{pmatrix} \bullet(\psi_1 ({\bf z} ) ,\psi_2({\bf z})) \quad = \quad \frac{\partial^3 \psi_1}{ \partial z_1^2 \partial z_3} + \frac{\partial^3 \psi_2}{ \partial z_1 \partial z_2 \partial z_4}.\end{equation}

The right-hand side is a scalar-valued function $\mathbb{R}^4 \rightarrow \mathbb{C}$ , that is, it is an element of $\mathcal{F}$ .

Our systems of PDE are submodules M of the free module $R^k$ . By Hilbert’s Basis Theorem, every module M is finitely generated, so we can write $M = {\textrm{image}}_R(A)$ , where A is a $k \times l$ matrix with entries in R. Each column of A is a generator of M, and it defines a differential operator that maps $\mathcal{F}^k$ to $\mathcal{F}$ . The solution space to the PDE given by M equals

(2.4) \begin{equation}{\textrm{Sol}}(M) \,:\!=\, \bigl\{ \psi \in \mathcal{F}^k \,:\, m \bullet \psi = 0\ \hbox{for all} \ m \in M \bigr\}.\end{equation}

It suffices to take the operators m from a generating set of M, such as the l columns of A. The case $k=1$ is of special interest, since we often consider PDE for scalar-valued functions. In that case, the submodule is an ideal in the polynomial ring R, and we denote this by I. The solution space ${\textrm{Sol}}(I)$ of the ideal $I \subseteq R$ is the set of functions $\phi$ in $\mathcal{F}$ such that $p(\partial) \bullet \phi = 0 $ for all $p \in I$ . Thus, ideals are instances of modules, with their own notation.

The solution spaces ${\textrm{Sol}}(M)$ and ${\textrm{Sol}}(I)$ are $\mathbb{C}$ -vector spaces and R-modules. Indeed, any $\mathbb{C}$ -linear combination of solutions is again a solution. The R-module action means applying the same differential operator $p(\partial)$ to each coordinate, which leads to another vector in $\mathcal{F}^k$ . This action takes solutions to solutions because the ring of differential operators with constant coefficients $R = \mathbb{C}[\partial_1,\ldots,\partial_n]$ is commutative.

The purpose of this paper is to present practical methods for the following task:

(2.5) \begin{equation} \begin{matrix}{ Given\, a\, k \times l\, matrix\, A \,with\, entries\, in\, R = K[\partial_1,\ldots,\partial_n],\, compute\, a\, good} \\[4pt]{ representation\, for\, the\, solution\, space\, {\textrm{Sol}}(M) \,of\, the\, module\, M = {\textrm{image}}_R(A).}\end{matrix}\end{equation}

If $k=1$ then we consider the ideal I generated by the entries of A and we compute ${\textrm{Sol}}(I)$ .

This raises the question of what a “good representation” means. The formulas in (1.2), (1.5), (1.9) and (1.11) are definitely good. They guide us to what is desirable. Our general answer stems from the following important result at the crossroads of analysis and algebra. It involves two sets of unknowns ${\bf z} = (z_1,\ldots,z_n)$ and ${\bf x}= (x_1,\ldots,x_n)$ . Here x gives coordinates on certain irreducible varieties $V_i$ in $\mathbb{C}^n$ that are parameter spaces for solutions. Our solutions $\psi$ are functions in z. We take $\mathcal{F} = C^\infty(\Omega)$ where $\Omega \subset \mathbb{R}^n$ is open, convex, and bounded.

Theorem 2.2 (Ehrenpreis--Palamodov Fundamental Principle). Consider a module $M \subseteq R^k$ , representing linear PDE for a function $\psi\,:\, \Omega \rightarrow \mathbb{C}^k$ . There exist irreducible varieties $V_1,\ldots,V_s$ in $\mathbb{C}^n$ and finitely many vectors $B_{ij}$ of polynomials in 2n unknowns $({\bf x},{\bf z})$ , all independent of the set $\Omega$ , such that any solution $\psi \in \mathcal{F}$ admits an integral representation

(2.6) \begin{equation} \psi(\mathbf{z}) = \sum_{i=1}^s \sum_{j=1}^{m_i} \int_{V_i} B_{ij} \left(\mathbf{x},\mathbf{z}\right) \exp\left( \mathbf{x}^t \mathbf{z} \right) d\mu_{ij} (\mathbf{x}).\end{equation}

Here $m_i$ is a certain invariant of $(M,V_i)$ and each $\mu_{ij}$ is a bounded measure supported on the variety $V_i $ .

Theorem 2.2 appears in different forms in the books by Björk [Reference Björk7, Theorem 8.1.3], Ehrenpreis [Reference Ehrenpreis17], Hörmander [Reference Hörmander23, Section 7.7], and Palamodov [Reference Palamodov32]. Other references with different emphases include [Reference Berndtsson and Passare5, Reference Lomadze25, Reference Oberst29]. For a perspective from commutative algebra see [Reference Cid-Ruiz, Homs and Sturmfels11, Reference Cid-Ruiz and Sturmfels12].

In the next sections, we will study the ingredients in Theorem 2.2. Given the module M, we compute each associated variety $V_i$ , the arithmetic length $m_i$ of M along $V_i$ , and the Noetherian multipliers $B_{i,1},B_{i,2},\ldots,B_{i,m_i}$ . We shall see that not all n of the unknowns $z_1,\ldots,z_n$ appear in the polynomials $B_{i,j}$ but only a subset of ${\textrm{codim}}(V_i)$ of them.

The most basic example is the ODE in (1.1), with $l=n=k=1$ . Here $V_i = \{u_i\}$ is the ith root of the polynomial (1.3), which has multiplicity $m_i$ , and $B_{i,j} = z^{j-1}$ . The measure $\mu_{ij}$ is a scaled Dirac measure on $u_i$ , so the integrals in (2.6) are multiples of the basis functions (1.2).

In light of Theorem 2.2, we now refine our computational task in (2.5) as follows:

(2.7) \begin{equation} \begin{matrix}{ Given \,a \,k \times l\, matrix \,A\, with\, entries\, in\, R = K[\partial_1,\ldots,\partial_n],\ \ compute\, the\, varieties\, V_i} \\[4pt] { and\, the\, Noetherian\, multipliers\, B_{ij}({\bf x},{\bf z}).\ \ This\, encodes\, {\textrm{Sol}}(M)\, for\, M = {\textrm{image}}_R(A).}\end{matrix}\end{equation}

In our introductory examples, we gave formulas for the general solution, namely (1.5) and (1.9). We claim that such formulas can be read off from the integrals in (2.6). For instance, for the wave equation (1.4), we have $s=2$ , $B_{1,1} = B_{1,2} = 1$ , and (1.5) is obtained by integrating ${\textrm{exp}}( {\bf x}^t {\bf z})$ against measures $ d\mu_{i1}({\bf x})$ on two lines $V_1$ and $V_2$ in $ \mathbb{C}^2$ . For the system (1.8), we find $s=6$ , with $m_1=m_2=m_3=1$ and $m_4=m_5=m_6=2$ , and the nine integrals in (2.6) translate into (1.9). We shall explain such a translation in full detail for two other examples.

Example 2.3 ( $n=3,k=1,l=2$ ) The ideal $I = \langle \partial_1^2 - \partial_2 \partial_3, \partial_3^2 \rangle $ represents the PDE

(2.8) \begin{equation} \frac{\partial^2 \phi}{\partial z_1^2} = \frac{\partial^2 \phi} {\partial z_2 \partial z_3}\qquad {\textrm{and}} \qquad \frac{\partial^2 \phi}{\partial z_3^2} = 0\end{equation}

for a scalar-valued function $\phi = \phi(z_1,z_2,z_3)$ . This is [Reference Chen, Härkönen, Krone and Leykin10, Example 4.2]. A Macaulay2 computation as in Section 5 shows that $s=1, m_1 =4$ . It reveals the Noetherian multipliers

\begin{equation*} B_1 = 1, B_2 = z_1, B_3 = z_1^2 x_2 + 2 z_3, B_4 = z_1^3 x_2 + 6 z_1 z_3 . \end{equation*}

Arbitrary functions $f(z_2) = \int {\textrm{exp}}( t z_2 ) dt$ are obtained by integrating against suitable measures on the line $V_1= \{ (0,t,0) \,:\,t \in \mathbb{C} \} \subset \mathbb{C}^3$ . Their derivatives are found by differentiating under the integral sign, namely $f'(z_2) = \int t \cdot {\textrm{exp}}( t z_2 )dt$ . Consider four functions a,b,c,d, each specified by a different measure. Thus, the sum of the four integrals in (2.6) evaluates to

(2.9) \begin{equation} \phi({\bf z}) = a(z_2) + z_1 b(z_2) + ( z_1^2 c'(z_2) + 2 z_3 c(z_2) ) + ( z_1^3 d'(z_2) + 6 z_1z_3 d(z_2)).\end{equation}

According to Ehrenpreis–Palamodov, this sum is the general solution of the PDE (2.8).

Our final example uses concepts from primary decomposition, to be reviewed in Section 3.

Example 2.4 ( $n=4,k=2,l=3$ ). Let $M \subset R^4$ be the module generated by the columns of

(2.10) \begin{equation} A = \begin{bmatrix} \partial_{1} \partial_{3} & \quad\partial_{1} \partial_{2} & \quad\partial_{1}^2 \partial_{2} \\[4pt] \partial_{1}^2 & \quad\partial_{2}^2 & \quad\partial_{1}^2 \partial_{4} \end{bmatrix}. \end{equation}

Computing ${\textrm{Sol}}(M)$ means solving $ \dfrac{\partial^2 \psi_1}{\partial z_1 \partial z_3} + \dfrac{\partial^2 \psi_2}{\partial z_1^2} = \dfrac{\partial^2 \psi_1}{\partial z_1 \partial z_2} + \dfrac{\partial^2 \psi_2}{\partial z_2^2} = \dfrac{\partial^3 \psi_1}{\partial z_1^2 \partial z_2} + \dfrac{\partial^3 \psi_2}{\partial z_1^2 \partial z_4} =0$ . Two solutions are $\psi({\bf z}) = \bigl(\phi(z_2,z_3,z_4) , 0\bigr)$ and $\psi({\bf z}) = {\textrm{exp}}(s^2 t z_1 + s t^2 z_2 + s^3 z_3 + t^3 z_4) \cdot \bigl( t , -s \bigr)$ .

We apply Theorem 2.2 to derive the general solution to (2.10). The module M has six associated primes, namely $P_1 = \langle \partial_{1} \rangle$ , $P_2 = \langle \partial_{2}, \partial_{4} \rangle $ , $P_3 = \langle \partial_{2}, \partial_{3} \rangle $ , $P_4 = \langle \partial_{1}, \partial_{3} \rangle $ , $P_5 = \langle \partial_{1}, \partial_{2} \rangle $ , and $P_6 = \langle \partial_{1}^2 - \partial_2 \partial_3, \partial_1 \partial_2 - \partial_3 \partial_4,\partial_2^2 - \partial_1 \partial_4 \rangle $ . Four of them are minimal and two are embedded. We find that $m_1 = m_2 = m_3 = m_4 = m_6 = 1$ and $m_5 = 4$ . A minimal primary decomposition

(2.11) \begin{equation} M = M_1 \cap M_2 \cap M_3 \cap M_4 \cap M_5 \cap M_6 \end{equation}

is given by the following primary submodules of $R^4$ , each of which contains M:

\begin{align*} M_1 & = {\textrm{im}}_R \begin{bmatrix} \partial_1 & \quad 0 \\[4pt] 0 & \quad 1 \end{bmatrix},\quad M_2 = {\textrm{im}}_R \begin{bmatrix} \partial_2 & \quad \partial_4 & \quad 0 & \quad 0 & \quad \partial_3 \\[4pt] 0 & \quad 0 & \quad \partial_2 & \quad \partial_4 & \quad \partial_1\end{bmatrix},\qquad M_3 = {\textrm{im}}_R \begin{bmatrix} \partial_2 & \quad \partial_3 & \quad 0 \\[4pt] 0 & \quad 0 & \quad 1 \end{bmatrix},\\[5pt]M_4 & = {\textrm{im}}_R \begin{bmatrix} \partial_3^5 & \quad \partial_1 & \quad 0 & \quad 0 \\[4pt] 0 & \quad \partial_2 & \quad \partial_1 & \quad \partial_3 \end{bmatrix},\ \ M_5 = {\textrm{im}}_R \begin{bmatrix} \partial_1 & \quad \partial_2^5 & \quad 0 & \quad 0 \\[4pt] 0 & \quad 0 & \quad \partial_1^2 & \quad \partial_2^2 \end{bmatrix},\ \ M_6 = {\textrm{im}}_R \begin{bmatrix}\partial_1 & \quad \partial_2 & \quad \partial_3 \\[4pt]\partial_2 & \quad \partial_4 & \quad \partial_1\end{bmatrix}. \end{align*}

The number of Noetherian multipliers $B_{ij}$ is $\sum_{i=1}^6 m_i = 9$ . We choose them to be

\begin{equation*} B_{1,1} {=} \begin{bmatrix} 1 \\[4pt] 0 \end{bmatrix} , B_{2,1} {=} \begin{bmatrix} \phantom{-}x_1 \\[4pt] -x_3 \end{bmatrix} , B_{3,1} {=} \begin{bmatrix} 1 \\[4pt] 0 \end{bmatrix} , B_{4,1} = \begin{bmatrix} x_2 z_1 \\[4pt] -1 \end{bmatrix} , B_{5,i} = \begin{bmatrix} 0 \\[4pt] z_1 z_2 \end{bmatrix} , \begin{bmatrix} 0 \\[4pt] z_1 \end{bmatrix} , \begin{bmatrix} 0 \\[4pt] z_2 \end{bmatrix} , \begin{bmatrix} 0 \\[4pt] 1 \end{bmatrix} , B_{6,1} = \begin{bmatrix} \phantom{-}x_4 \\[4pt] -x_2 \end{bmatrix}.\end{equation*}

These nine vectors describe all solutions to our PDE. For instance, $B_{3,1}$ gives the solutions $ \Big[\begin{array}{c} \alpha(z_1,z_4) \\[4pt] 0 \end{array}\Big]$ , and $B_{5,1}$ gives the solutions $ \Big[\begin{array}{c} 0 \\[4pt] z_1 z_2 \beta(z_3,z_4) \end{array}\Big]$ , where $\alpha,\beta$ are bivariate functions. Furthermore, $B_{1,1}$ and $B_{6,1}$ encode the two families of solutions mentioned after (2.10).

For the latter, we note that $V_6 = V(P_6) $ is the surface in $\mathbb{C}^4$ with parametric representation $(x_1,x_2,x_3,x_4) = (s^2 t, st^2, s^3,t^3)$ for $s,t \in \mathbb{C}$ . This surface is the cone over the twisted cubic curve, in the same notation as in [Reference Cid-Ruiz, Homs and Sturmfels11, Section 1]. The kernel under the integral in (2.6) equals

\begin{equation*} \begin{bmatrix} \phantom{-}x_4 \\[4pt] -x_2 \end{bmatrix}{\textrm{exp}}\bigl(x_1 z_1 + x_2 z_2 + x_3 z_3 + x_4 z_4\bigr) \quad = \quad t^2 \begin{bmatrix} \phantom{-} t \\[4pt] - s \end{bmatrix} {\textrm{exp}}\bigl( s^2t z_1 + st^2 z_2 + s^3 z_3 + t^3 z_4 \bigr).\end{equation*}

This is a solution to $M_6$ , and hence to M, for any values of s and t. Integrating the left- hand side over ${\bf x} \in V_6$ amounts to integrating the right-hand side over $(s,t) \in \mathbb{C}^2$ . Any such integral is also a solution. Ehrenpreis–Palamodov tells us that these are all the solutions.

3. Modules and varieties

Our aim is to offer practical tools for solving PDE. The input is a $k \times l$ matrix A with entries in $R = K[\partial_1,\ldots,\partial_n]$ , and $M = {\textrm{image}}_R(A)$ is the corresponding submodule of $R^k = \bigoplus_{j=1}^k Re_j$ . The output is the description of ${\textrm{Sol}}(M)$ sought in (2.7). That description is unique up to basis change, in the sense of [Reference Cid-Ruiz and Sturmfels12, Remark 3.8], by the discussion in Section 4. Our method is implemented in a Macaulay2 command, called solvePDE and to be described in Section 5.

We now explain the ingredients of Theorem 2.2 coming from commutative algebra (cf. [Reference Eisenbud18]). For a vector $m \in R^k$ , the quotient $(M\,:\,m)$ is the ideal $\{f \in R \,:\, fm \in M\}$ . A prime ideal $P_i \subseteq R$ is associated to M if there exists $m \in R^k$ such that $(M\,:\,m) = P_i$ . Since R is Noetherian, the list of associated primes of M is finite, say $P_1,\ldots,P_s$ . If $s=1$ then the module M is called primary or $P_1$ -primary. A primary decomposition of M is a list of primary submodules $M_1,\ldots,M_s \subseteq R^k$ where $M_i$ is $P_i$ -primary and $ M = M_1 \cap M_2 \cap \cdots \cap M_s $ .

Primary decomposition is a standard topic in commutative algebra. It is usually presented for ideals $(k=1)$ , as in [Reference Michałek and Sturmfels27, Chapter 3]. The case of modules is analogous. The latest version of Macaulay2 has an implementation of primary decomposition for modules, as described in [Reference Chen and Cid-Ruiz9, Section 2]. Given M, the primary module $M_i$ is not unique if $P_i$ is an embedded prime.

The contribution of the primary module $M_i$ to M is quantified by a positive integer $m_i$ , called the arithmetic length of M along $P_i$ . To define this, we consider the localization $(R_{P_i})^k/M_{P_i}$ . This is a module over the local ring $R_{P_i}$ . The arithmetic length is the length of the largest submodule of finite length in $(R_{P_i})^k/M_{P_i}$ ; in symbols, $m_i = {\textrm{length}} \bigl( H^0_{P_i} ((R_{P_i})^k/ M_{P_i})\bigr)$ . The sum $m_1 + \cdots + m_s$ is an invariant of the module M, denoted ${\textrm{amult}}(M)$ , and known as the arithmetic multiplicity of M. These numbers can be computed in Macaulay2 as in [Reference Cid-Ruiz and Sturmfels12, Remark 5.1]. We return to these invariants in Theorem 4.3.

To make the connection to Theorem 2.2, we set $V_i = V(P_i)$ for $i=1,2,\ldots,s$ . Thus, $V_i$ is the irreducible variety in $\mathbb{C}^n$ defined by the prime ideal $P_i$ in $R = K[\partial_1,\ldots,\partial_n]$ . The integer $m_i$ is an invariant of the pair $(M,V_i)$ : it measures the thickness of the module M along $V_i$ .

By taking the union of the irreducible varieties $V_i$ we obtain the variety

\begin{equation*} V(M) \quad \,:\!=\, \quad V_1 \cup V_2 \cup \cdots \cup V_s\quad \subset \mathbb{C}^n.\end{equation*}

Algebraists refer to V(M) as the support of M, while analysts call it the characteristic variety of M. The support is generally reducible, with $\leq s$ irreducible components. For instance, the module M in Example 2.4 has six associated primes, and an explicit primary decomposition was given in (2.11). However, the support V(M) has only four irreducible components in $\mathbb{C}^4$ , namely one hyperplane, two-dimensional planes, and one nonlinear surface (twisted cubic).

Remark 3.1. If $k=1$ and $M=I$ , then the support V(M) coincides with the variety V(I) attached as usual to an ideal I, namely the common zero set in $\mathbb{C}^n$ of all polynomials in I.

The relationship between modules and ideals mirrors the relationship between PDE for vector-valued functions and related PDE for scalar-valued functions. To pursue this a bit further, we now define two ideals that are naturally associated with a given module $M\subseteq R^k$ .

The first ideal is the annihilator of the quotient module $R^k/M = {\textrm{coker}}_R(A)$ , which is

\begin{equation*} I \,:\!=\, {\textrm{Ann}}_R(R^k/M) = \big\{ f \in R \,:\, f m \in M\ \hbox{for all}\ m \in R^k \bigr\} . \end{equation*}

The second is the zeroth Fitting ideal of $R^k/M$ , which is the ideal in R generated by the $k \times k$ minors of the presentation matrix A. It is independent of the choice of A, and we write

\begin{equation*} J \,:\!=\, {\textrm{Fitt}}_0(R^k/M) = \bigl\langle\hbox{$k \times k$ subdeterminants of A}\bigr\rangle. \end{equation*}

We are interested in the affine varieties in $\mathbb{C}^n$ defined by these ideals. They are denoted by V(I) and V(J), respectively. The following is a standard result in commutative algebra.

Proposition 3.2. The three varieties above are equal for every submodule M of $R^k$ , that is,

(3.1) \begin{equation} V(M) = V(I) = V(J) \subseteq \mathbb{C}^n. \end{equation}

Proof. This follows from [Reference Eisenbud18, Proposition 20.7].

Remark 3.3. It can happen that ${\textrm{rank}}(A) < k$ , for instance when $k > l$ . In that case, $I = J = \{ 0 \}$ and $V(M) = \mathbb{C}^n$ . Geometrically, the module M furnishes a coherent sheaf that is supported on the entire space $\mathbb{C}^n$ . For instance, let $k=n=2,l=1$ , and $A = \displaystyle\binom{\phantom{-}\partial_1}{-\partial_2}$ . The PDE asks for pairs $(\psi_1,\psi_2)$ such that $\partial \psi_1 /\partial z_1 = \partial \psi_2 /\partial z_2 $ . We see that $\textrm{Sol}(M)$ consists of all pairs $\bigl( \partial \alpha/ \partial z_2,\partial \alpha/ \partial z_1 \big)$ , where $\alpha= \alpha(z_1,z_2)$ runs over functions in two variables. In general, the left kernel of A furnishes differential operators for creating solutions to M.

The following example shows that (3.1) is not true at the level of schemes (cf. Section 6).

Example 3.4. ( $n=k=3,l=5$ ) Let $R = \mathbb{C}[\partial_1,\partial_2,\partial_3]$ and M the submodule of $R^3$ given by

\begin{equation*} A = \begin{pmatrix}\partial_1 & \quad 0 & \quad 0 & \quad 0 & \quad 0 \\[4pt] 0 & \quad \partial_1^2 & \quad \partial_2 & \quad 0 & \quad 0 \\[4pt] 0 & \quad 0 & \quad 0 & \quad \partial_1 & \quad \partial_3 \end{pmatrix} . \end{equation*}

We find $I = \langle \partial_1^2, \partial_1 \partial_2 \rangle \supset J = \langle \partial_1^4, \partial_1^3 \partial_3,\partial_1^2 \partial_2 , \partial_1 \partial_2 \partial_3 \rangle$ . The sets of associated primes are

\begin{equation*} \begin{matrix}& \textrm{Ass}(I) & = & \bigl\{ \langle \partial_1 \rangle, \langle \partial_1, \partial_2 \rangle \bigr\}&\qquad & {\textrm{with}}\ {\textrm{amult}}(I) =2 \\[4pt] \subset & \textrm{Ass}(M) & = & \bigl\{ \langle \partial_1 \rangle, \langle \partial_1, \partial_2 \rangle , \langle \partial_1, \partial_3 \rangle \bigr\} & \qquad & {\textrm{with}}\ {\textrm{amult}}(M) = 4 \\[4pt] \subset & \textrm{Ass}(J) & = & \bigl\{ \langle \partial_1 \rangle, \langle \partial_1, \partial_2 \rangle , \langle \partial_1, \partial_3 \rangle, \langle \partial_1, \partial_2, \partial_3 \rangle \bigr\} & \qquad & {\textrm{with}}\ {\textrm{amult}}(J) = 5 \end{matrix} \end{equation*}

The support V(M) is a plane in 3-space, on which I and J define different scheme structures. Our module M defines a coherent sheaf on that plane that lives between these two schemes. We consider the PDE in each of the three cases, we compute the Noetherian multipliers, and from this we derive the general solution. To begin with, functions in $ {\textrm{Sol}}(J)$ have the form

\begin{equation*} \alpha(z_2,z_3) + z_1 \beta(z_3) + z_1^2 \gamma(z_3) +z_1 \delta(z_2) + c \cdot z_1^3 . \end{equation*}

The first two terms give functions in the subspace ${\textrm{Sol}}(I)$ . Elements in ${\textrm{Sol}}(M)$ are vectors

\begin{equation*} \bigl( \rho(z_2,z_3) , \sigma(z_3) + z_1 \tau(z_3), \omega(z_2) \bigr) . \end{equation*}

These represent all functions $\mathbb{C}^3 \rightarrow \mathbb{C}^3$ that satisfy the five PDE given by the matrix A.

Remark 3.5. The quotient $R/I$ embeds naturally into the direct sum of k copies of $R^k/M$ , via $1 \mapsto e_j$ . This implies ${\textrm{Ass}}(I) \subseteq {\textrm{Ass}}(M)$ . It would be worthwhile to understand how the differential primary decompositions of I,J and M are related, and to study implications for the solution spaces ${\textrm{Sol}}(I)$ , ${\textrm{Sol}}(J)$ , and ${\textrm{Sol}}(M)$ . What relationships hold between these?

Lemma 3.6. Fix a $k \times l$ matrix $A(\partial)$ and its module $M \subseteq R^k$ as above. A point ${\bf u}\in \mathbb{C}^n$ lies in V(M) if and only if there exist constants $c_1,\ldots,c_k \in \mathbb{C}$ , not all zero, such that

(3.2) \begin{align} \begin{pmatrix} c_1 \\[4pt] \vdots \\[4pt] c_k \end{pmatrix} \exp(u_1z_1 + \dotsb + u_nz_n) \in \textrm{Sol}(M). \end{align}

More precisely, (3.2) holds if and only if $ (c_1,\ldots,c_k) \cdot A({\bf u}) = 0 $ .

Proof. Let $a_{ij}(\partial)$ denote the entries of the matrix $A(\partial)$ . Then (3.2) holds if and only if

\begin{align*} \sum_{i = 1}^k a_{ij}(\partial) \bullet (c_i \exp(u_1z_1+\dotsb + u_nz_n)) = 0 \quad \text{ for all } j = 1,\ldots, l. \end{align*}

This is equivalent to

\begin{align*} \sum_{i=1}^k c_i a_{ij}({\bf u}) \exp(u_1z_1+\dotsb + u_nz_n) = 0 \quad \text{ for all } j = 1,\ldots,l. \end{align*}

This condition holds if and only if $ (c_1,\ldots,c_k) \cdot A({\bf u}) $ is the zero vector in $\mathbb{C}^l$ . We conclude that, for any given ${\bf u} \in \mathbb{C}^n$ , the previous condition is satisfied for some $c \in \mathbb{C}^k \backslash \{0\}$ if and only if ${\textrm{rank}}(A({\bf u})) < k$ if and only if ${\bf u} \in V(M) = V(I)$ . Here we use Proposition 3.2.

Here is an alternative way to interpret the characteristic variety of a system of PDE:

Proposition 3.7. The solution space $\textrm{Sol}(M)$ contains an exponential solution $ q({\bf z}) \cdot {\textrm{exp}}( {\bf u}^t {\bf z}) $ if and only if ${\bf u} \in V(M)$ . Here q is some vector of k polynomials in n unknowns, as in (2.2).

Proof. One direction is clear from Lemma 3.6. Next, suppose $ q(\mathbf{z}) \exp({\bf u}^t {\bf z}) \in \textrm{Sol}(M)$ . The partial derivative of this function with respect to any unknown $z_i$ is also in ${\textrm{Sol}}(M)$ . Hence,

\begin{align*} \partial_i \bullet (q({\bf z}) \exp({\bf u}^t {\bf z})) = (\partial_i \bullet q({\bf z})) \exp({\bf u}^t {\bf z})+ u_i q({\bf z}) \exp({\bf u}^t {\bf z}) \in \textrm{Sol}(M) \quad \hbox{for $i=1,\ldots,n$.} \end{align*}

Hence, the exponential function $ (\partial_i \bullet q({\bf z})) \exp({\bf u}^t {\bf z})$ is in ${\textrm{Sol}}(M)$ . Since the degree of $\partial_i \bullet q({\bf z})$ is less than that of $q({\bf z})$ , we can find a sequence $D = \partial_{i_1} \partial_{i_2} \dotsb \partial_{i_s}$ such that $D \bullet q$ is a nonzero constant vector and $(D \bullet q) \exp({\bf u}^t {\bf z}) \in \textrm{Sol}(M)$ . Lemma 3.6 now implies that ${\bf u} \in V(M)$ .

The solution space ${\textrm{Sol}}(M)$ to a submodule $M \subseteq R^k$ is a vector space over $\mathbb{C}$ . It is infinite-dimensional whenever V(M) is a variety of positive dimension. This follows from Lemma 3.6 because there are infinitely many points u in V(M). However, if V(M) is a finite subset of $\mathbb{C}^n$ , then ${\textrm{Sol}}(M)$ is finite-dimensional. This is the content of the next theorem.

Theorem 3.8. Consider a module $M \subseteq R^k$ , viewed as a system of linear PDE. Its solution space $\textrm{Sol}(M)$ is finite-dimensional over $\mathbb{C}$ if and only if V(M) has dimension 0. In this case, ${\textrm{dim}}_\mathbb{C} \textrm{Sol}(M) = {\textrm{dim}}_K(R^k/M) = {\textrm{amult}}(M)$ . There is a basis of ${\textrm{Sol}}(M)$ given by vectors $ q({\bf z}) {\textrm{exp}}({\bf u}^t {\bf z})$ , where ${\bf u} \in V(M)$ and $q({\bf z})$ runs over a finite set of polynomial vectors, whose cardinality is the length of M along the maximal ideal $\langle x_1 - u_1,\ldots,x_n-u_n \rangle$ . There exist polynomial solutions if and only if $\mathfrak{m} = \langle x_1,\ldots,x_n \rangle$ is an associated prime of M. The polynomial solutions are found by solving the PDE given by the $\mathfrak{m}$ -primary component of M.

Proof. This is the main result in Oberst’s article [Reference Oberst30], proved in the setting of injective cogenerators $\mathcal{F}$ . The same statement for $\mathcal{F} = C^\infty(\Omega)$ appears in [Reference Björk7, Ch. 8, Theorem 7.1]. The scalar case $(k=1)$ is found in [Reference Michałek and Sturmfels27, Theorem 3.27]. The proof given there uses solutions in the power series ring, which is an injective cogenerator, and it generalizes to modules.

By a polynomial solution we mean a vector $q({\bf z})$ whose coordinates are polynomials. The $\mathfrak{m}$ -primary component in Theorem 3.8 is computed by a double saturation step. When $M=I$ is an ideal, then this double saturation is $I\,:\,(I\,:\,\mathfrak{m}^\infty)$ , as seen in [Reference Michałek and Sturmfels27, Theorem 3.27]. For submodules M of $R^k$ with $k \geq 2$ , we would compute $ M \,:\, {\textrm{Ann}}(R^k / (M \,:\, \mathfrak{m}^\infty) ) $ . The inner colon $(M\,:\,\mathfrak{m}^\infty)$ is the intersection of all primary components of M whose variety $V_i$ does not contain the origin 0. It is computed as $(M\,:\,f) = \{m \in R^k\,:\, fm \in M \}$ , where f is a random homogeneous polynomial of large degree. The outer colon is the module $(M\,:\,g)$ , where g is a general polynomial in the ideal $\textrm{Ann}(R^k/(M\,:\,f))$ . See also [Reference Chen and Cid-Ruiz9, Proposition 2.2].

It is an interesting problem to identify polynomial solutions when V(M) is no longer finite and to decide whether these are dense in the infinite-dimensional space of all solutions. Here “dense” refers to the topology on $\mathcal{F}$ used by Lomadze in [Reference Lomadze26]. The following result gives an algebraic characterization of the closure in ${\textrm{Sol}}(M)$ of the subspace of polynomial solutions.

Proposition 3.9. The polynomial solutions are dense in ${\textrm{Sol}}(M)$ if and only if the origin 0 lies in every associated variety $V_i$ of the module M. If this fails, then the topological closure of the space of polynomial solutions $q({\bf z})$ to M is the solution space of $M \,:\, \textrm{Ann}(R^k/(M \,:\, \mathfrak{m}^\infty))$ .

Proof. This proposition is our reinterpretation of Lomadze’s result in [Reference Lomadze26, Theorem 3.1].

The result gives rise to algebraic algorithms for answering analytic questions about a system of PDE. The property in the first sentence can be decided by running the primary decomposition algorithm in [Reference Chen and Cid-Ruiz9]. For the second sentence, we need to compute a double saturation as above. This can be carried out in Macaulay2 as well.

4. Differential primary decomposition

We now shift gears and pass to a setting that is dual to the one we have seen so far. Namely, we discuss differential primary decompositions [Reference Chen and Cid-Ruiz9, Reference Cid-Ruiz and Sturmfels12]. That duality is subtle and can be confusing at first sight. To mitigate this, we introduce new notation. We set $x_i = \partial_i = \partial_{z_i}$ for $i=1,\ldots,n$ . Thus, R is now the polynomial ring $K[x_1,\ldots,x_n]$ . This is the notation we are used to from algebra courses (such as [Reference Michałek and Sturmfels27]). We write $\partial_{x_1},\ldots,\partial_{x_n}$ for the differential operators corresponding to $x_1,\ldots,x_n$ . Later on, we also identify $z_i = \partial_{x_i}$ , and we think of the unknowns x and z in the multipliers $B_i({\bf x},{\bf z})$ as dual in the sense of the Fourier transform.

The ring of differential operators on the polynomial ring R is the Weyl algebra

\begin{equation*} D_n = R \langle \partial_{x_1},\ldots,\partial_{x_n} \rangle =K \langle x_1,\ldots,x_n, \partial_{x_1},\ldots,\partial_{x_n} \rangle . \end{equation*}

The 2n generators commute, except for the n relations $\partial_{x_i} x_i - x_i \partial_{x_i} = 1$ , which expresses the Product Rule from Calculus. Elements in the Weyl algebra $D_n$ are linear differential operators with polynomial coefficients. We write $\delta \bullet p$ for the result of applying $\delta \in D_n$ to a polynomial $p = p({\bf x})$ in R. For instance, $x_i \bullet p = x_i p$ and $\partial_{x_i} \bullet p = \partial p/\partial x_i$ . Let $D_n^k$ denote the k-tuples of differential operators in $D_n$ . These operate on the free module $R^k$ as follows:

\begin{equation*} D_n^k \times R^k \rightarrow R \,:\, (\delta_1,\ldots,\delta_k) \bullet (p_1,\ldots,p_k)= \sum_{j=1}^k \delta_j \bullet p_j. \end{equation*}

Fix a submodule M of $R^k$ and let $P_1,\ldots,P_s$ be its associated primes, as in Section 3. A differential primary decomposition of M is a list $\mathcal{A}_1,\ldots,\mathcal{A}_s$ of finite subsets of $D_n^k$ such that

(4.1) \begin{equation}M = \bigl\{ m \in R^k \,:\, \delta \bullet m \in P_i\ \ \hbox{for all}\ \delta \in \mathcal{A}_i\ \hbox{and}\ i=1,2, \ldots,s \bigr\}.\end{equation}

This is a membership test for the module M using differential operators. This test is geometric since the polynomial $\delta \bullet m $ lies in $P_i $ if and only if it vanishes on the variety $V_i=V(P_i)$ .

Theorem 4.1 Every submodule M of $R^k$ has a differential primary decomposition. We can choose the sets $\mathcal{A}_1,\ldots,\mathcal{A}_s$ such that $|\mathcal{A}_i|$ is the arithmetic length of M along the prime $P_i$ .

Proof and discussion. The result is proved in [Reference Cid-Ruiz and Sturmfels12] and further refined in [Reference Chen and Cid-Ruiz9]. These sources also develop an algorithm. We shall explain this in Section 5, along with a discussion of the Macaulay2 command solvePDE, which computes differential primary decompositions.

The differential operators in $\mathcal{A}_1,\ldots,\mathcal{A}_s$ are known as Noetherian operators in the literature; see [Reference Chen, Härkönen, Krone and Leykin10, Reference Cid-Ruiz, Homs and Sturmfels11, Reference Damiano, Sabadini and Struppa15, Reference Oberst31]. Theorem 4.1 says that we can find a collection of $\textrm{amult}(M) = m_1 + \cdots + m_s$ Noetherian operators in $D_n^k$ to characterize membership in the module M.

Remark 4.2 The construction of Noetherian operators is studied in [Reference Björk7, Reference Chen, Cid-Ruiz, Härkönen, Krone and Leykin8, Reference Chen, Härkönen, Krone and Leykin10, Reference Cid-Ruiz, Homs and Sturmfels11, Reference Hörmander23, Reference Oberst31]. Some of these sources offer explicit methods, while others remain at an abstract level. All previous methods share one serious shortcoming, namely they yield operators separately for each primary component $M_i$ of M. They do not take into account how one primary component is embedded into another. This leads to a number of operators that can be much larger than amult (M). We refer to [Reference Cid-Ruiz and Sturmfels12, Example 5.6] for an instance from algebraic statistics where the previous methods require 1044 Noetherian operators, while ${\textrm{amult}}(M) = 207$ suffice.

While Theorem 4.1 makes no claim of minimality, it is known that $\textrm{amult}(M)$ is the minimal number of Noetherian operators required for a differential primary decomposition of a certain desirable form. To make this precise, we begin with a few necessary definitions. For any given subset $\mathcal{S}$ of $\{x_1,\ldots,x_n\}$ , the relative Weyl algebra is defined as the subring of the Weyl algebra $D_n$ using only differential operators corresponding to variables not in $\mathcal{S}$ :

(4.2) \begin{equation} D_n(\mathcal{S}) \,:\!=\, R \langle \partial_{x_i} \,:\, x_i \not \in \mathcal{S} \rangle \subseteq D_n.\end{equation}

Thus, if $\mathcal{S} = \emptyset$ , then $D_n(\mathcal{S}) = D_n$ , and if $\mathcal{S} = \{x_1,\ldots,x_n\}$ , then $D_n(\mathcal{S}) = R= K[x_1,\ldots,x_n]$ .

For any prime ideal $P_i$ in R we fix a set $\mathcal{S}_i \subseteq \{x_1,\ldots,x_n\}$ that satisfies $K[\mathcal{S}_i] \cap P_i = \{0\}$ and is maximal with this property. Thus, $\mathcal{S}_i$ is a maximal independent set of coordinates on the irreducible variety $V(P_i)$ . Equivalently, $\mathcal{S}_i$ is a basis of the algebraic matroid defined by the prime $P_i$ ; cf. [Reference Michałek and Sturmfels27, Example 13.2]. The cardinality of $\mathcal{S}_i$ equals the dimension of $V(P_i)$ .

Theorem 4.3 The differential primary decomposition in Theorem 4.1 can be chosen so that $\mathcal{A}_i \subset D_n(\mathcal{S}_i)^k$ . The arithmetic length of M along $P_i$ is a lower bound for the cardinality of $\mathcal{A}_i$ in any differential primary decomposition of M such that $\mathcal{A}_i \subset D_n(\mathcal{S}_i)^k$ for $i = 1,\ldots,s$ .

Proof and discussion. This was shown in [Reference Cid-Ruiz and Sturmfels12, Theorem 4.6]. The case of ideals $(k=1)$ appears in [Reference Cid-Ruiz and Sturmfels12, Theorem 3.6]. See also [Reference Chen and Cid-Ruiz9]. The theory developed in [Reference Cid-Ruiz and Sturmfels12] is more general in that R can be any Noetherian K-algebra. In this paper, we restrict to polynomial rings $R = K[x_1,\ldots,x_n]$ where K is a subfield of $\mathbb{C}$ . That case is treated in detail in [Reference Chen and Cid-Ruiz9].

We next argue that Theorems 2.2 and 4.1 are really two sides of the same coin. Every element A in the Weyl algebra $D_n$ acts as a differential operator with polynomial coefficients on functions in the unknowns ${\bf x} = (x_1,\ldots,x_n)$ . Such a differential operator has a unique representation where all derivatives are moved to the right of the polynomial coefficients:

(4.3) \begin{equation}\qquad A({\bf x},\partial_{\bf x}) = \sum_{{\bf r},{\bf s} \in \mathbb{N}^n} c_{{\bf r},{\bf s}} x_1^{r_1} \cdots x_n^{r_n} \partial_{x_1}^{s_1} \cdots \partial_{x_n}^{s_n} , \qquad{\textrm{where}} c_{{\bf r},{\bf s}} \in K .\end{equation}

There is a natural K-linear isomorphism between the Weyl algebra $D_n$ and the polynomial ring $K[{\bf x},{\bf z}]$ which takes the operator A in (4.3) to the following polynomial B in 2n variables:

(4.4) \begin{equation} B({\bf x},{\bf z}) = \sum_{{\bf r},{\bf s} \in \mathbb{N}^n} c_{{\bf r},{\bf s}} x_1^{r_1} \cdots x_n^{r_n} z_1^{s_1} \cdots z_n^{s_n} . \qquad \qquad \quad \end{equation}

In Sections 1, 2, and 3, polynomials in $R = K[x_1,\ldots,x_n] = K[\partial_1,\ldots,\partial_n] $ act as differential operators on functions in the unknowns ${\bf z} = (z_1,\ldots,z_n)$ . For such operators, polynomials in x are constants. By contrast, in the current section, we introduced the Weyl algebra $D_n$ . Its elements act on functions in ${\bf x} = (x_1,\ldots,x_n)$ , with polynomials in z being constants. These two different actions of differential operators, by $D_n$ and R on scalar-valued functions, extend to actions by $D_n^k$ and $R^k$ on vector-valued functions. We highlight the following key point:

(4.5) \begin{equation} {Our\, distinction\, between\, the\, {z} -variables\, and\, {x}-variables\, is \, absolutely\, essential.}\end{equation}

The following theorem is the punchline of this section. It allows us to identify Noetherian operators (4.3) with Noetherian multipliers (4.4). This was assumed tacitly in [Reference Cid-Ruiz, Homs and Sturmfels11, Section 3].

Theorem 4.4 Consider any differential primary decomposition of the module M as in Theorem 4.3. Then this translates into an Ehrenpreis–Palamodov representation of the solution space ${\textrm{Sol}}(M)$ . Namely, if we replace each operator $A({\bf x},\partial_{\bf x})$ in $\mathcal{A}_i$ by the corresponding polynomial $B({\bf x},{\bf z})$ , then these ${\textrm{amult}}(M)$ polynomials satisfy the conclusion of Theorem 2.2.

Example 4.5 ( $k=l=n=1$ ) We illustrate Theorem 4.4 and the warning (4.5) for an ODE (1.1) with $m=3$ . Set $ p(x) = x^3 + 3 x^2 - 9x + 5 = (x-1)^2 (x+5)$ in (1.3). The ideal $I = \langle p \rangle$ has $s=2$ associated primes in $R = \mathbb{Q}[x] $ , namely $P_1 = \langle x-1 \rangle $ and $P_2 = \langle x+5 \rangle$ , with $m_1 = 2$ and $m_2=1$ , so ${\textrm{amult}}(I) = 3$ . A differential primary decomposition of I is given by $\mathcal{A}_1 = \{1,\partial_x \}$ and $\mathcal{A}_2 = \{1\}$ . The three Noetherian operators translate into the Noetherian multipliers $B_{11} = 1, B_{12} = z, B_{21} = 1$ . The integrals in (2.6) now furnish the general solution $ \phi(z) = \alpha {\textrm{exp}}(z) + \beta z {\textrm{exp}}(z) + \gamma {\textrm{exp}}({-}5z) $ to the differential equation $\phi''' + 3 \phi'' - 9 \phi' + 5 \phi = 0$ .

The derivation of Theorem 4.4 rests on the following lemma on duality between x and z.

Lemma 4.6 Let p and q be polynomials in n unknowns with coefficients in K. We have

(4.6) \begin{equation} q(\partial_{\bf z}) \bullet \bigl(p({\bf z}) \exp({\bf x}^t {\bf z}) \bigr) = p(\partial_{\bf x}) \bullet \bigl( q({\bf x}) \exp({\bf x}^t {\bf z}) \bigr).\end{equation}

Proof. The parenthesized expression on the left equals $p(\partial_{\bf x}) \bullet {\textrm{exp}}({\bf x}^t {\bf z})$ , while that on the right equals $q(\partial_{\bf z}) \bullet {\textrm{exp}}({\bf x}^t {\bf z})$ . Therefore, the expression in (4.6) is the result of applying the operator $p(\partial_{\bf x}) q(\partial_{\bf z}) = q(\partial_{\bf z}) p(\partial_{\bf x}) $ to $ \exp({\bf x}^t {\bf z}) $ , when viewed as a function in 2n unknowns.

We now generalize this lemma to $k \geq 2$ , we replace p by a polynomial vector that depends on both x and z, and we rename that vector using the identification between (4.3) and (4.4).

Proposition 4.7 Let $B({\bf x},{\bf z})$ be a k-tuple of polynomials in 2n variables and $A({\bf x},\partial_{\bf x}) \in D_n^k$ the corresponding k-tuple of differential operators in the Weyl algebra. Then we have

(4.7) \begin{align} q(\partial_\mathbf{z}) \bullet (B(\mathbf{x}, \mathbf{z}) \exp(\mathbf{x}^t \mathbf{z})) = A(\mathbf{x}, \partial_\mathbf{x}) \bullet (q(\mathbf{x}) \exp(\mathbf{x}^t \mathbf{z})) . \end{align}

Proof. If $k=1$ , we write $A(\mathbf{x}, \partial_\mathbf{x}) = \sum_{\alpha} c_\alpha(\mathbf{x}) \partial_\mathbf{x}^\alpha$ as in (4.3) and $B({\bf x},{\bf z}) = \sum_{\alpha} c_\alpha(\mathbf{x}) \mathbf{z}^\alpha$ as in (4.4). Only finitely many of the polynomials $c_\alpha(\mathbf{x})$ are nonzero. Applying Lemma 4.6 gives

\begin{equation*} \begin{matrix} A(\mathbf{x}, \partial_\mathbf{x}) \bullet (q(\mathbf{x}) \exp(\mathbf{x}^t \mathbf{z})) = \sum_\alpha c_\alpha(\mathbf{x}) q(\partial_\mathbf{z}) \bullet ({\bf z}^\alpha \exp(\mathbf{x}^t \mathbf{z})) = q(\partial_\mathbf{z}) \bullet (B(\mathbf{x}, \mathbf{z}) \exp(\mathbf{x}^t \mathbf{z})). \end{matrix}\end{equation*}

The extension from $k=1$ to $k \geq 2$ follows because the differential operation $\bullet$ is K-linear.

We now take a step toward proving Theorem 4.4 in the case $s=1$ . Let M be a primary submodule of $R^k$ with $\textrm{Ass}(M) = \{P\}$ . Its support $V(M) = V(P)$ is an irreducible affine variety in $\mathbb{C}^n$ . Consider the sets of all Noetherian operators and all Noetherian multipliers:

(4.8) \begin{align} \mathfrak{A} & \,:\!= \bigl\{ A \in D_n^k\,:\, A \bullet m \in P\ \hbox{for all}\ m \in M \bigr\} \qquad \qquad \qquad \qquad \qquad {\textrm{and}} \nonumber\\[4pt] \mathfrak{B} & \,:\!= \bigl\{B \in K[{\bf x},{\bf z}]\,:\, B({\bf x},{\bf z}) {\textrm{exp}}({\bf x}^t {\bf z}) \in {\textrm{Sol}}(M)\hbox{for all ${\bf x} \in V(P) $} \bigr\}.\end{align}

Proposition 4.8 The bijection between $D_n^k$ and $K[{\bf x},{\bf z}]^k$ , given by identifying the operator A in (4.3) with the polynomial B in (4.4), restricts to a bijection between the sets $\mathfrak{A}$ and $\mathfrak{B}$ .

Proof. Let $m_1,\ldots,m_l \in K[{\bf x}]^k$ be generators of M. Suppose $A \in \mathfrak{A}$ . Then

\begin{align*} \sum_{i = 1}^k A_i(\mathbf{x}, \partial_\mathbf{x}) \bullet \sum_{j = 1}^l m_{ij}(\mathbf{x}) f_j(\mathbf{x}) \end{align*}

vanishes for all $\mathbf{x} \in V(P)$ and all polynomials $f_1,\ldots,f_l \in \mathbb{C}[\mathbf{x}]$ . Since the space of complex-valued polynomials is dense in the space of all entire functions on $\mathbb{C}^n$ , the preceding implies

\begin{equation*} \sum_{i = 1}^k A_i(\mathbf{x}, \partial_\mathbf{x}) \bullet m_{ij}(\mathbf{x}) \exp(\mathbf{x}^t\mathbf{z}) = 0 \quad \hbox{for all $\mathbf{z} \in \mathbb{C}^n$, $\mathbf{x} \in V(P)$ and $j = 1,\ldots,l$.}\end{equation*}

Using Proposition 4.7, this yields

\begin{equation*} \sum_{i = 1}^k m_{ij}(\partial_\mathbf{z}) \bullet B_i(\mathbf{x}, \mathbf{z}) \exp(\mathbf{x}^t\mathbf{z})= 0 \quad \hbox{for all $\mathbf{z} \in \mathbb{C}^n$, $\mathbf{x} \in V(P)$ and $j = 1,\ldots,l$.} \end{equation*}

We conclude that the polynomial vector $B(\mathbf{x},\mathbf{z}) $ corresponding to $A({\bf x},\partial_{\bf x})$ lies in $\mathfrak{B}$ .

To prove the converse, we note that the implications above are reversible. Thus, if $B({\bf x},{\bf z}) $ is in $\mathfrak{B}$ , then $A({\bf x},\partial_{\bf x})$ is in $\mathfrak{A}$ . This uses the fact that linear combinations of the exponential functions $\mathbf{x} \to \exp(\mathbf{x}^t \mathbf{z})$ , for ${\bf z} \in \mathbb{C}^n$ , are also dense in the space of entire functions.

Proof of Theorem 4.4. Let $\mathcal{A}$ be any finite subset of $\mathfrak{A}$ which gives a differential primary decomposition of the P-primary module M. This exists and can be chosen to have cardinality equal to the length of M along P. Let $\mathcal{B}$ be the set of Noetherian multipliers (4.4) corresponding to the set $\mathcal{A}$ of Noetherian operators (4.3). Proposition 4.8 shows that the exponential function ${\bf z} \to B({\bf x},{\bf z}) {\textrm{exp}}({\bf x}^t {\bf z})$ is in ${\textrm{Sol}}(M)$ whenever ${\bf x} \in V(P)$ and $B \in \mathcal{B}$ . Hence all $\mathbb{C}$ -linear combinations of such functions are in ${\textrm{Sol}}(M)$ . More generally, by differentiating under the integral sign, we find that all functions of the following form are solutions of M:

\begin{equation*} \psi({\bf z}) = \sum_{B \in \mathcal{B}} \int_{V(P)} B({\bf x},{\bf z}) {\textrm{exp}}({\bf x}^t {\bf z}) d\mu_B({\bf x}).\end{equation*}

We need to argue that all solutions in $\mathcal{F} = C^\infty(\Omega)$ admit such an integral representation. Suppose first that all associated primes of M are minimal. Then each $\mathcal{A}_i$ spans a bimodule in the sense of [Reference Chen and Cid-Ruiz9, Theorem 3.2 (d)]. Hence, for each associated prime $P_i$ , the module

\begin{align*} M_i = \{ m \in R^k \,:\, \delta \bullet m \in P_i \text{ for all } \delta \in \mathcal{A}_i \}\end{align*}

is $P_i$ -primary, and $M = M_1 \cap \dotsb \cap M_s$ is a minimal primary decomposition. The operators in $\mathcal{A}_i$ are in the relative Weyl algebra $D_n(\mathcal{S}_i)$ and fully characterize the $P_i$ -primary component of M. We may thus follow the classical analytical constructions in the books [Reference Björk7, Reference Hörmander23, Reference Palamodov32] to patch together the integral representation of ${\textrm{Sol}}(M_i)$ for $i=1,\ldots,s$ , under the correspondence of Noetherian operators and Noetherian multipliers. Therefore, all solutions have the form (2.6).

Things are more delicate when M has embedded primes. Namely, if $P_i$ is embedded, then the operators in $\mathcal{A}_i$ only characterize the contribution of the $P_i$ -primary component relative to all other components contained in $P_i$ . We see this in Section 5. One argues by enlarging $\mathcal{A}_i$ to vector space generators of the relevant bimodule. Then the previous patching argument applies. And, afterward one shows that the added summand in the integral representation are redundant because they are covered by associated varieties $V(P_j)$ containing $V(P_i)$ .

5. Software and algorithm

In this section, we present an algorithm for solving linear PDE with constant coefficients. It is based on the methods for ideals given in [Reference Chen, Cid-Ruiz, Härkönen, Krone and Leykin8, Reference Chen, Härkönen, Krone and Leykin10, Reference Cid-Ruiz, Homs and Sturmfels11]. The case of modules appears in [Reference Chen and Cid-Ruiz9]. We note that the computation of Noetherian operators has a long history, going back to work in the 1990’s by Ulrich Oberst [Reference Oberst28, Reference Oberst29, Reference Oberst30, Reference Oberst31], who developed a construction of Noetherian operators for primary modules. This was further developed by Damiano, Sabadini and Struppa [Reference Damiano, Sabadini and Struppa15] who presented the first Gröbner-based algorithm. It works for primary ideals under the restrictive assumption that the characteristic variety has a rational point after passing to a (algebraically nonclosed) field of fractions. Their article also points to an implementation in CoCoA, but we were unable to access that code. Since these early approaches rely on the ideals or modules being primary, using them in practice requires first computing a primary decomposition. If there are embedded primes, the number of Noetherian operators output by these methods will not be minimal either.

We here present a new algorithm that is universally applicable, to all ideals and modules over a polynomial ring. There are no restrictions on the input and the output is minimal. The input is a submodule M of $R^k$ , where $R = K[x_1,\ldots,x_n]$ . The output is a differential primary decomposition of size ${\textrm{amult}}(M)$ as in Theorem 4.3. A first step is to find ${\textrm{Ass}}(M) = \{P_1,\ldots,P_s\}$ . For each associated prime $P_i$ , the elements $A({\bf x},\partial_{\bf x})$ in the finite set $\mathcal{A}_i \subset D_n(\mathcal{S}_i)$ are rewritten as polynomials $B({\bf x},{\bf z})$ , using the identification of (4.3) with (4.4). Only the ${\textrm{codim}}(P_i)$ many variables $z_i $ with $x_i \not\in \mathcal{S}_i$ appear in these Noetherian multipliers B.

We now describe our implementation for (2.7) in Macaulay2 [Reference Grayson and Stillman21]. The command is called solvePDE, as in [Reference Cid-Ruiz and Sturmfels12, Section 5]. It is distributed with Macaulay2 starting from version 1.18 in the package NoetherianOperators [Reference Chen, Cid-Ruiz, Härkönen, Krone and Leykin8]. The user begins by fixing a polynomial ring $R = K[x_1,\ldots,x_n]$ . Here K is usually the rational numbers $ \texttt{QQ}$ . Fairly arbitrary variable names $x_i$ are allowed. The argument of solvePDE is an ideal in R or a submodule of $R^k$ . The output is a list of pairs $\bigl\{P_i,\{B_{i1},\ldots,B_{i,m_i}\} \bigr\}$ for $i=1,\ldots,s$ , where $P_i$ is a prime ideal given by generators in R, and each $B_{ij}$ is a vector over a newly created polynomial ring $K[x_1,\ldots,x_n,z_1,\ldots,z_n]$ . The new variables $z_i$ are named internally by Macaulay2. The system writes $ \texttt{d} x_i$ for $z_i$ . To be precise, each new variable is created from an old variable by prepending the character $ \texttt{d}$ . This notation can be confusing at first, but one gets used to it. The logic comes from the differential primary decompositions described in [Reference Cid-Ruiz and Sturmfels12, Section 5].

Each $B_{ij}$ in the output of solvePDE encodes an exponential solution $B_{ij}({\bf x},{\bf z}) {\textrm{exp}}({\bf x}^t {\bf z})$ to M. Here x are the old variables chosen by the user, and x denotes points in the irreducible variety $V(P_i) \subseteq \mathbb{C}^n$ . The solution is a function in the new unknowns ${\bf z} = (\texttt{d}x_1,\ldots, \texttt{d}x_n)$ . For instance, if $n=3$ and the input is in the ring $ \texttt{QQ[u,v,w]}$ , then the output lives in the ring $ \texttt{QQ[u,v,w,du,dv,dw]}$ . Each solution to the PDE is a function $\psi(\texttt{du}, \texttt{dv}, \texttt{dw})$ , and these functions are parametrized by a variety $V(P_i)$ in a 3-space whose coordinates are $(\texttt{u}, \texttt{v}, \texttt{w})$ .

We now demonstrate how this works for two examples featured in the introduction.

Example 5.1. Consider the third order ODE (1.1) in Example 4.5. We solve this as follows:

R = QQ[x]; I = ideal(x^3 + 3*x^2 - 9*x + 5); solvePDE(I)

{{ideal(x - 1), {| 1 |, | dx |}}, {ideal(x + 5), {| 1 |}}}

The first line is the input. The second line is the output created by solvePDE. This list of $s=2$ pairs encodes the general solution $\phi(z)$ . Remember: z is the newly created symbol dx.

Example 5.2. We solve the PDE (1.8) by typing the $ 2 \times 3$ matrix whose columns are (2.1):

R = QQ[x1,x2,x3,x4];

M = image matrix {{x1^2,x2*x3,x1^2*x3},{x1*x2,x3^2,x1*x2*x4}}; solvePDE(M)

The reader is encouraged to run this code and to check that the output is the solution (1.9).

The method in solvePDE is described in Algorithm 1 below. A key ingredient is a translation map. We now explain this in the simplest case, when the module is supported in one point. Suppose $V(M) = \{\mathbf{u}\}$ for some $\mathbf{u} \in K^n$ . We set $ \mathfrak{m}_\mathbf{u} = \langle x_1 - u_1, \ldots, x_n - u_n \rangle$ and

(5.1) \begin{equation} \gamma_\mathbf{u} \,:\, R \to R , x_i \mapsto x_i + u_i \quad \text{ for } i=1,\ldots,n. \end{equation}

The following two results are straightforward. We will later use them when M is any primary module, u is the generic point of V(M), and $\mathbb{K} = K({\bf u})$ is the associated field extension of K.

Proposition 5.3 A constant coefficient operator $A(\partial_\mathbf{x})$ is a Noetherian operator for the $\mathfrak{m}_\mathbf{u}$ -primary module M if and only if $A(\partial_\mathbf{x})$ is a Noetherian operator for the $\mathfrak{m}_0$ -primary module $\hat M \,:\!=\, \gamma_\mathbf{u}(M)$ . Dually, $B(\mathbf{z}) \exp(\mathbf{u}^t\mathbf{z}) $ is in $ \textrm{Sol}(M)$ if and only if $B(\mathbf{z})$ is in $ \textrm{Sol}(\hat M)$ .

We note that all Noetherian operators over a K-rational point can be taken to have constant coefficients. This follows from Theorem 3.8. This observation reduces the computation of solutions for a primary module to finding the polynomial solutions of the translated module. Next, we bound the degrees of these polynomials.

Proposition 5.4 Let $\hat M \subseteq R^k$ be an $\mathfrak{m}_0$ -primary module. There exists an integer r such that $\mathfrak{m}_0^{r+1} R^k \subseteq \hat M$ . The space $\textrm{Sol}(\hat M)$ consists of k-tuples of polynomials of degree $\leq r$ .

Propositions 5.3 and 5.4 furnish a method for computing solutions of an $\mathfrak{m}_\mathbf{u}$ -primary module M. We start by translating M so that it becomes the $\mathfrak{m}_0$ -primary module $\hat M$ . The integer r provides an ansatz $ \sum_{j=1}^k\sum_{|\alpha| \leq r} v_{\alpha,j} \mathbf{z}^\alpha e_j $ for the polynomial solutions. The coefficients $v_{\alpha,j}$ are computed by linear algebra over the ground field K. Here are the steps:

  1. 1. Let r be the smallest integer such that $\mathfrak{m}_0^{r+1} R^k \subseteq \hat M$ .

  2. 2. Let $\textrm{Diff}( \hat M)$ be the matrix whose entries are the polynomials $\hat m_i \bullet ({\bf z}^\alpha e_j) \in R$ . The row labels are the generators $\hat m_1, \ldots, \hat m_l $ of $\hat M$ , and the column labels are the ${\bf z}^\alpha e_j$ .

  3. 3. Let $\ker_K(\textrm{Diff}(\hat M))$ denote the K-linear subspace of the R-module $\ker_R(\textrm{Diff}(\hat M))$ consisting of vectors $(v_{\alpha,j}) $ with all entries in K. Every such vector gives a solution

    (5.2) \begin{align} \sum_{j=1}^k \sum_{|\alpha| \leq r} v_{\alpha, j} {\bf z}^\alpha \exp(\mathbf{u}^t \mathbf{z}) e_j \in \textrm{Sol}(M).\end{align}

Example 5.5 $[n=k=r=2]$ The following module is $\mathfrak{m}_0$ -primary of multiplicity three:

(5.3) \begin{equation}M = {{image}}_R \begin{bmatrix}\partial_1^3 & \quad \partial_2 - c_1 \partial_1^2 - c_2 \partial_1 & \quad c_3 \partial_1^2 + c_4 \partial_1 + c_5 \\[4pt] 0 & \quad 0 & \quad 1 \\[4pt]\end{bmatrix}.\end{equation}

Here $c_1,c_2,c_3,c_4,c_5$ are arbitrary constants in K. The matrix ${\textrm{Diff}}(M)$ has three rows, one for each generator of M, and it has 12 columns, indexed by $e_1,z_1e_1,\ldots , z_2^2 e_1,e_2,z_1e_2,\ldots , z_2^2 e_2$ . The space ${\textrm{ker}}_K({\textrm{Diff}}(M))$ is 3-dimensional. A basis furnishes the three polynomial solutions

(5.4) \begin{equation}\begin{bmatrix}-1 \\[4pt] c_5,\end{bmatrix} ,\begin{bmatrix} -(z_1+c_2 z_2) \\[4pt] c_5 z_1+c_2 c_5 z_2+c_4 \end{bmatrix} , \begin{bmatrix} -((z_1+c_2 z_2)^2 + 2 c_1 z_2) \\[4pt] c_5 (z_1{+}c_2 z_2)^2 + 2 c_4z_1 + 2 (c_1 c_5{+}c_2 c_4)z_2 + 2c_3 \end{bmatrix}.\end{equation}

We now turn to Algorithm 1. The input and output are as described in (2.7). The method was introduced in [Reference Chen and Cid-Ruiz9, Algorithm 4.6] for computing differential primary decompositions. We use it for solving PDE. It is implemented in Macaulay2 under the command solvePDE. In our discussion, the line numbers refer to the corresponding lines of pseudocode in Algorithm 1.

  • Line 1 We begin by finding all associated primes of M. These define the irreducible varieties $V_i$ in (2.7). By [Reference Eisenbud, Huneke and Vasconcelos19, Theorem 1.1], the associated primes of codimension i coincide with the minimal primes of $\textrm{Ann}\ \textrm{Ext}^i_R(M,R)$ . This reduces the problem of finding associated primes of a module to the more familiar problem of finding minimal primes of a polynomial ideal. This method is implemented and distributed with Macaulay2 starting from version 1.17 via the command associatedPrimes $\texttt{R^k/M}$ . See [Reference Chen and Cid-Ruiz9, Section 2].

Algorithm 1 SolvePDE

The remaining steps are repeated for each $P \in \textrm{Ass}(M)$ . For a fixed associated prime P, our goal is to identify the contribution to $\textrm{Sol}(M)$ of the P-primary component of M.

  • Lines 2–3 To achieve this goal, we study solutions for two different R-submodules of $R^k$ . The first one, denoted U, is the intersection of all $P_i$ -primary components of M, where $P_i$ are the associated primes contained in P. Thus $U = MR^k_P \cap R^k$ , which is the extension-contraction module of M under localization at P. It is computed as $U = (M \,:\, f^\infty)$ , where $f \in R$ is contained in every associated prime $P_j$ not contained in P.

    The second module, denoted V, is the intersection of all $P_i$ -primary components of M, where $P_i$ is strictly contained in P. Hence, $V = (U \,:\, P^\infty)$ is the saturation of U at P. We have $U = V \cap Q$ , where Q is a P-primary component of M. Thus, the difference between the solution spaces $\textrm{Sol}(U)$ and $\textrm{Sol}(V)$ is caused by the primary module Q.

    When P is a minimal prime, U is the unique P-primary component of M, and $V = R^k$ .

  • Line 4 The integer r bounds the degree of Noetherian multipliers associated to U but not V. Namely, if the function $\phi(\mathbf{z}) = B(\mathbf{x}, \mathbf{z}) \exp(\mathbf{x}^t \mathbf{z})$ lies in $\textrm{Sol}(U) \backslash \textrm{Sol}(V)$ for all $\mathbf{x} \in V(P)$ , then the $\mathbf{z}$ -degree of the polynomial $B(\mathbf{x}, \mathbf{z})$ is at most r. This will lead to an ansatz for the Noetherian multipliers responsible for the difference between $\textrm{Sol}(U)$ and $\textrm{Sol}(V)$ .

  • Lines 5–8 The modules U and V are reduced to simpler modules $\hat U$ and $\hat V$ with similar properties. Namely, $\hat U$ and $\hat V$ are primary and their characteristic varieties are the origin. This reduction involves two new ingredients: a new polynomial ring T in fewer variables over a field $\mathbb{K}$ that is a finite extension of K, and a ring map $\gamma \,:\, R \to T$ .

    Fix a maximal set $\mathcal{S} = \{x_{i_1},\ldots,x_{i_{n-c}}\}$ with $P\cap K[x_{i_1},\ldots,x_{i_{n-c}}]=\{0\}$ . We define $T\,:\!=\, \mathbb{K}[y_i\,:\, x_i\notin \mathcal{S}]$ , where $\mathbb{K}=\textrm{Frac}(R/P)$ . This is a polynomial ring in $n - |\mathcal{S}| = c$ new variables $y_i$ , corresponding to the $x_i$ not in the set $\mathcal{S}$ of independent variables. Writing $u_i$ for the image of $x_i$ in $\mathbb{K}=\textrm{Frac}(R/P)$ , the ring map $\gamma$ is defined as follows:

    (5.5) \begin{equation}\gamma \,:\, R \to T, \quad x_i\mapsto \begin{cases} y_i+u_i, & \text{ if } x_i\notin S,\\[4pt] \quad u_i, & \text{ if } x_i\in S. \end{cases} \end{equation}
    By abuse of notation, we denote by $\gamma$ the extension of (5.5) to a map $R^k \to T^k$ .
  • Lines 9–11 Let $\mathfrak{m} \,:\!=\, \langle y_i \,:\, x_i \not \in \mathcal{S} \rangle $ be the irrelevant ideal of T. We define the T-submodules

    \begin{equation*} \hat U \,:\!=\, \gamma(U) + \mathfrak{m}^{r+1} T^k \quad {\textrm{and}} \quad \hat V \,:\!=\, \gamma(U) + \mathfrak{m}^{r+1} T^k \quad {\textrm{of}}\ T^k. \end{equation*}
    These modules are $\mathfrak{m}$ -primary: their solutions are finite-dimensional $\mathbb{K}$ -vector spaces consisting of polynomials of degree $\leq r$ . The polynomials in $\textrm{Sol}(\hat U) \backslash \textrm{Sol}(\hat V)$ capture the difference between $\hat U$ and $\hat V$ , and also the difference between U and V after lifting.
  • Lines 12–14 We construct matrices $\textrm{Diff}(\hat U)$ and $\textrm{Diff}(\hat V)$ with entries in $\mathbb{K} [z_i \,:\, x_i \not\in \mathcal{S}] $ . As in (5.2), their kernels over $\mathbb{K}$ correspond to polynomial solutions of $\hat U$ and $\hat V$ . The set $N = \{{\bf z}^\alpha e_j \,:\, |\alpha| \leq r, j=1,\ldots, k\}$ is a $\mathbb{K}$ -basis for elements of degree $\leq r$ in $\mathbb{K}[z_i \,:\, x_i \not\in \mathcal{S}]^k$ . The $y_i$ -variables act on the $z_i$ variables as partial derivatives, i.e. $y_i = \dfrac{\partial}{\partial z_i}$ . We define the matrix $\textrm{Diff}(\hat U)$ as follows. Let $\hat{U}_1,\ldots, \hat{U}_\ell$ be generators of $\hat{U}$ . The rows of $\textrm{Diff}(\hat{U})$ are indexed by these generators, the columns are indexed by N, and the entries are the polynomials $\hat{U}_i\bullet {\bf z}^\alpha e_j$ . In the same way, we construct $\textrm{Diff}(\hat{V})$ .

  • Lines 15–16 Let $\ker_\mathbb{K}(\textrm{Diff}(\hat U))$ be the space of vectors in the kernel of $\textrm{Diff}(\hat U)$ whose entries are in $\mathbb{K}$ . The $\mathbb{K}$ -vector space $\ker_\mathbb{K}(\textrm{Diff}(\hat U))$ parametrizes the polynomial solutions

    \begin{align*} \sum_{j=1}^k \sum_{|\alpha| \leq r} v_{\alpha,j} {\bf z}^\alpha e_j \in \textrm{Sol}(\hat U). \end{align*}
    The same holds for $\hat V$ . The quotient space $\mathcal{K}\,:\!=\, \ker_\mathbb{K}(\textrm{Diff}(\hat U)) / \ker_\mathbb{K}(\textrm{Diff}(\hat V))$ characterizes excess solutions in $\textrm{Sol}(\hat U)$ relative to $\textrm{Sol}(\hat V)$ . Write $\mathcal{A}$ for a $\mathbb{K}$ -basis of $\mathcal{K}$ .
  • Lines 17–18 We interpret $\mathcal{A}$ as a set of Noetherian multipliers for M by performing a series of lifts and transformations. For each element $\bar{\mathbf{v}} \in \mathcal{A}$ , we choose a representative $\mathbf{v} \in \ker_\mathbb{K}(\textrm{Diff}(\hat U))$ . The entries of $\mathbf{v}$ are in $\mathbb{K} = \textrm{Frac}(R/P)$ and may contain denominators. Multiplying $\mathbf{v}$ by a common multiple of the denominators yields a vector with entries in $R/P$ , indexed by N. We lift this to a vector $\mathbf{u} = (u_{\alpha,j})$ with entries in R. The Noetherian multiplier corresponding to u is the following vector in $R[\texttt{d}x_i \,:\, x_i \not \in \mathcal{S}]^k$ :

    \begin{align*} B(\mathbf{x}, \mathbf{\texttt{d}x}) = \sum_{j=1}^k \sum_{|\alpha| \leq r} u_{\alpha,j}(\mathbf{x}) \mathbf{\texttt{d}x}^\alpha e_j . \end{align*}
    Applying the map $\bar{\mathbf{v}} \mapsto \mathbf{u}$ to each $\bar{\mathbf{v}} \in \mathcal{A}$ yields a set $\mathcal{B}$ of Noetherian multipliers. These multipliers describe the contribution of the P-primary component of M to ${\textrm{Sol}}(M)$ .

The output of Algorithm 1 is a list of pairs $(P, \mathcal{B})$ , where P ranges over ${\textrm{Ass}}(M)$ and $\mathcal{B} = \{B_1,\ldots,B_m\}$ is a subset of $R[\texttt{d}x_1, \ldots, \texttt{d}x_n]^k$ . The cardinality m is the multiplicity of M along P. The output describes the solutions to the PDE given by M. Consider the functions

\begin{equation*}\phi_P(\texttt{d}x_1,\ldots,\texttt{d}x_n)= \sum_{i=1}^m \int_{V(P)}B_i(\mathbf{x}, \mathbf{\texttt{d}x})\exp(x_1\texttt{d}x_1+\cdots+x_n\texttt{d}x_n){d\mu_i}(\mathbf{x}).\end{equation*}

Then the space of solutions to M consists of all functions

\begin{equation*}\sum_{P \in {\textrm{Ass}}(M)} \phi_P(\texttt{d}x_1,\ldots,\texttt{d}x_n).\end{equation*}

A differential primary decomposition of M is obtained from this by reading $\texttt{d}x_i$ as $\partial_{x_i}$ . Indeed, the command differentialPrimaryDecomposition described in [Reference Chen and Cid-Ruiz9] is identical to our command solvePDE. All examples in [Reference Chen and Cid-Ruiz9, Section 6] can be interpreted as solving PDE.

6. Schemes and coherent sheaves

The concepts of schemes and coherent sheaves are central to modern algebraic geometry. These generalize varieties and vector bundles, and they encode geometric structures with multiplicities. The point is that the supports of coherent sheaves and other schemes are generally nonreduced. We here argue that our linear PDE offer a useful way to think about the geometry of these objects. That perspective motivated the writing of [Reference Michałek and Sturmfels27, Section 3.3].

The affine schemes we consider are defined by ideals I in a polynomial ring R. Likewise, submodules M of $R^k$ represent coherent sheaves on $\mathbb{C}^n$ . We study the affine scheme ${\textrm{Spec}}(R/I)$ and the coherent sheaf given by the module $R^k/M$ . The underlying geometric objects are the affine varieties V(I) and V(M) in $\mathbb{C}^n$ . The latter was discussed in Section 3. The solution spaces ${\textrm{Sol}}(I)$ or ${\textrm{Sol}}(M)$ furnish nonreduced structures on these varieties, encoded in the integral representations due to Ehrenpreis–Palamodov. According to Section 4, these are dual to differential primary decompositions. Coherent sheaves were a classical tool in the analysis of linear PDE, but in the analytic category, where their role was largely theoretical. The Ehrenpreis–Palamodov Fundamental Principle appears in Hörmander’s book under the header Coherent analytic sheaves on Stein manifolds [Reference Hörmander23, Chapter VII]. Likewise, Treves’ exposition, in the title of [Reference Treves35, Section 3.2], calls for Analytic sheaves to the rescue. By contrast, sheaves in this paper are concrete and algebraic: they are modules in Macaulay2.

One purpose of this section is to explore how PDE and their solutions behave under degenerations. We consider ideals and modules whose generators depend on a parameter $\epsilon$ . This is modelled algebraically by working over the field $K = \mathbb{C}(\epsilon)$ of rational functions in the variable $\epsilon$ . Algorithm 1 can be applied to the polynomial ring $R = K[x_1,\ldots,x_n]$ over that field. We think of $\epsilon$ as a small quantity and we are interested in what happens when $\epsilon \rightarrow 0$ .

Our discussion in this section is very informal. This is by design. We present a sequence of examples that illustrates the geometric ideas. The only formal result is Theorem 6.6, which concerns the role of the Quot scheme in parametrizing systems of linear PDE.

Example 6.1 ( $n=2$ ) Consider the prime ideal $I_\epsilon = \langle \partial_1^2 - \epsilon^2 \partial_2 \rangle$ . For nonzero parameters $\epsilon $ , by Theorem 2.2, the solutions to this PDE are represented as one-dimensional integrals

\begin{equation*} \alpha_\epsilon(z_1,z_2) = \int {\textrm{exp}} (\epsilon\ t\ z_1 + t^2\ z_2) dt \quad \in {\textrm{Sol}}(I). \end{equation*}

By taking the limit for $\epsilon \rightarrow 0$ , this yields arbitrary functions $a(z_2)$ . These are among the solutions to $I_0 = \langle \partial_1^2 \rangle $ . Other limit solutions are obtained via the reflection $t \mapsto -t$ . Set

\begin{equation*} \beta_\epsilon(z_1,z_2) = \int {\textrm{exp}} (-\epsilon\ t\ z_1 + t^2\ z_2) dt \quad \in {\textrm{Sol}}(I). \end{equation*}

Note the similarity to the one-dimensional wave equation (1.5) with $c = \epsilon$ . The solution for $\epsilon = 0$ is given in (1.6). This is found from the integrals above by taking the following limit:

(6.1) \begin{equation}\begin{aligned} {\textrm{lim}}_{\epsilon \rightarrow 0} \frac{1}{2 \epsilon} \bigl( \alpha_\epsilon(z_1,z_2) - \beta_\epsilon(z_1,z_2) \bigr) & = \int {\textrm{lim}}_{\epsilon \rightarrow 0}\frac{{\textrm{exp}}(\epsilon tz_1+t^2z_2)-{\textrm{exp}}(-\epsilon tz_1+t^2z_2)}{2\epsilon}dt\\[4pt] &= \int tz_1{\textrm{exp}}(t^2z_2) dt\\[4pt]& = z_1 b(z_2). \end{aligned}\end{equation}

We conclude that the general solution to $I_0$ equals $ \phi(z_1,z_2) = a(z_2) + z_1 b(z_2)$ , where b is any function in one variable. The calculus limit in (6.1) realizes a scheme-theoretic limit in the sense of algebraic geometry. Namely, two lines in (1.7) converge to a double line in $\mathbb{C}^2$ .

Example 6.2 ( $n=3$ ). For $\epsilon \not= 0$ consider the curve $t \mapsto (\epsilon t^3, t^4, \epsilon^2 t^2)$ in $\mathbb{C}^3$ . Its prime ideal equals $I_\epsilon = \langle \partial_1^2 - \partial_2 \partial_3, \partial_3^2 - \epsilon^4 \partial_2 \rangle$ . The solution space ${\textrm{Sol}}(I_\epsilon)$ consists of the functions

(6.2) \begin{equation} \phi(z_1,z_2,z_3) = \int {\textrm{exp}}( \epsilon t^3 z_1 + t^4 z_2 + \epsilon^2 t^2 z_3) dt. \end{equation}

What happens to these functions when $\epsilon $ tends to zero? We address this question algebraically. The scheme-theoretic limit of the given ideal $I_\epsilon$ is the ideal in Example 2.3. This is verified by a Gröbner basis computation (cf. [Reference Eisenbud18, Section 15.8]). Passing from ideals to their varieties, we see a toric curve in $\mathbb{C}^3$ that degenerates to a line with multiplicity four.

We claim that the formula in (2.9) arises from (6.2), just as in Example 6.1. Namely, set $i = \sqrt{-1}$ and let $\phi_s \in {\textrm{Sol}}(I_\epsilon)$ be the function that is obtained from $\phi$ in (6.2) by replacing the parameter t with $ i^s t$ . Then the following four functions on the left are solutions to $I_\epsilon$ :

\begin{equation*}\begin{matrix} \phi_0+\phi_1+\phi_2+\phi_3 & \longrightarrow & a(z_2), \\[4pt] \epsilon^{-1} (\phi_0+i \phi_1+i^2 \phi_2+i^3 \phi_3) & \longrightarrow & z_1 b(z_2), \\[4pt] \epsilon^{-2} (\phi_0+i^2 \phi_1+i^4 \phi_2+i^6 \phi_3) & \longrightarrow & z_1^2 c'(z_2) + 2 z_3 c(z_2), \\[4pt] \epsilon^{-3} (\phi_0+i^3 \phi_1+i^6 \phi_2+i^9 \phi_3) & \longrightarrow & z_1^3 d'(z_2) + 6 z_1z_3 d(z_2).\end{matrix}\end{equation*}

The functions obtained as limits on the right are precisely the four summands seen in (2.9). Thus, the solution spaces to this family of PDE reflect the degeneration of the toric curve.

Such limits make sense also for modules. If a module $M_\epsilon \subseteq R^k$ depends on a parameter $\epsilon$ then we study its solution space ${\textrm{Sol}}(M_\epsilon)$ as $\epsilon $ tends to zero. Geometrically, we examine flat families of coherent sheaves on $\mathbb{C}^n$ or on $\mathbb{P}^{n-1}$ . A typical scenario comes from the action of the torus $(\mathbb{C}^*)^n$ , where Gröbner degenerations arise as limits under one-parameter subgroups. The limit objects are monomial ideals (for $k=1$ ) or torus-fixed submodules (for $k \geq 2$ ). The next example illustrates their rich structure with an explicit family of torus-fixed submodules.

Example 6.3 ( $n=2, k=3, l=6$ ) Given a $3 \times 6$ matrix A with random real entries, we set

\begin{equation*} M = {\textrm{image}}_R \bigl( A \cdot {\textrm{diag}}(\partial_1,\partial_1^2,\partial_1^3, \partial_2,\partial_2^2,\partial_2^3)\bigr) \quad \subset \quad R^3. \end{equation*}

Then M is torus-fixed and $\mathfrak{m}$ -primary, where $\mathfrak{m} = \langle \partial_1,\partial_2 \rangle$ , and ${\textrm{amult}}(M) = 10$ . A basis of ${\textrm{Sol}}(M)$ is given by ten polynomial solutions, namely the standard basis vectors $e_1,e_2,e_3$ , four vectors that are multiples of $z_1,z_1,z_2, z_2$ , and three vectors that are multiples $z_1^2, z_1 z_2, z_2^2$ . The reader is invited to verify this with Macaulay2. Here is the input for one concrete instance:

R = QQ[x1,x2]

M = image matrix {{7*x1,5*x1^2,8*x1^3, 5*x2,9*x2^2,5*x2^3},

{8*x1,9*x1^2,8*x1^3, 4*x2,2*x2^2,4*x2^3},

{3*x1,2*x1^2,6*x1^3, 4*x2,4*x2^2,7*x2^3}}

solvePDE(M)

By varying the matrix A, and by extracting the vector multipliers of $1,z_1$ and