lajos/SDIRK p2 s3 s4

Some observations concerning the SDIRK class with \(p=2\), and \(3\le s \le 4\).

Let us fix a positive integer \(s\ge 3\).

Suppose that \(\psi:=\frac{P}{Q^s}\) approximates \(\exp\) at \(0\) to second order \(p=2\) and suppose that \(Q\) has the form \(Q(x)=1-a x\). The Griend-Kraaijevanger article implies that \(a>0\) is necessary to have \(ModRad(\psi)>0\).

The order conditions are equivalent to the following form of \(\psi\):

\[ \psi(x)=\frac{1+\left(1-a \binom{s}{1}\right)x+ \left(\frac{1}{2}-a \binom{s}{1}+a^2 \binom{s}{2}\right)x^2+\sum_{n=3}^{s} a_n x^n}{(1-a x)^{s} } \]

with \(a_n\in\mathbb{R}\) (\(n=3,\ldots,s\)).

Now if \(k\ge 0\) is any integer, we have

\[ \psi^{(k)}=\frac{k!}{Q^{s+k}}\sum_{m=0}^{\min(k,s)} \frac{1}{m!}\binom{s-1+k-m}{s-1} \ P^{(m)}\ Q^m \ \left(-Q\prime\right)^{k-m}\quad\quad\quad (*) \]

Notice that for \(x\le 0\), \(Q(x)>0\), so for the absolute monotonicity of \(\psi\) it is enough to study only the sign of the sum. The sum has a bounded number of terms even if \(k\to \infty\).

Conjecture. For any choice of \(a>0\) and \(a_n\in \mathbb{R}\), we have \(MonRad(\psi)\le 2s\), further, \(MonRad(\psi)=2s\) holds if and only if

\[ P(x)=\left(1+\frac{x}{2s}\right)^s, \quad Q(x)=1-\frac{x}{2s}, \]

that is, if \(a=\frac{1}{2s}\) and \(a_n=\frac{\binom{s}{n}}{(2s)^n}\) for \(n\ge 3\). (Of course, if we denote by \(a_0\), \(a_1\) and \(a_2\) the appropriate coefficients in \(P\), then this last relation holds for \(n=0, 1, 2\) also.)

Remark\(_1\). Notice that in the above we have not used the non-negativity of matrix \(A\). So it seems that under (6.35a\(-\)d.), \(A\ge 0\) is responsible for forcing that the stability function has a pole of maximal possible order \(s\), and although \(A\ge 0\) has some other implications on \(a_n\) (\(n\ge 3\)), the conjecture above seems to be true without any further structural assumptions on the numerator (apart from the ones on \(a_0, a_1\) and \(a_2\) implied by the order conditions).

Remark\(_2\). I have not yet verified whether the given unique form of \(P\) and \(Q\) really implies that \(A\) has the conjectured structure, i.e., \(\frac{1}{2s}\) in the diagonal and \(\frac{1}{s}\) all below, further, \(b=\frac{1}{s}\cdot \mathbf{1}\).

Some notations:

PC(\(x_0, k_0\)) This point condition says that the sum \(\sum_{m=0}^{\min(k,s)}\) in (*) is non-negative if \(x=x_0\) and \(k=k_0\ge 0\).

IC(\([\alpha,\beta], k_0\)) This interval condition says that the sum \(\sum_{m=0}^{\min(k,s)}\) in (*) is non-negative if \(k=k_0\) and for all \(x\in [\alpha,\beta]\).

LC This limit condition says that \[ \lim_{k\to\infty} \frac{1}{k^{s-1}} \sum_{m=0}^{\min(k,s)} \frac{1}{m!}\binom{s-1+k-m}{s-1} \ P^{(m)}\ Q^m \ \left(-Q\prime\right)^{k-m}\ge 0. \] (Since the summand is a polynomial in \(k\) of degree at most \(s-1\), the above limit exists and finite.)

Observations and partial proofs that support the above conjecture.

The \(s=3\) case.

\[ \psi(x)=\frac{1+\left(1-3a\right)x+ \left(\frac{1}{2}-3a+3a^2\right)x^2+c x^3}{(1-a x)^{3} } \]

\(-\) In my earlier Mma calculations (I already reported) I have taken into account the fact that \(A\ge 0\) implies that \(c\) has additional structure.

\(-\) Now I could prove with Mma (no human-readable proof yet) that

PC\((-6,0)\), PC\((-6,1)\), PC\((-6,2)\) and IC\(([-6,0],4)\) imply \(a=\frac{1}{6}\) and \(c=\frac{1}{216}\).

I consider this as a proof of the SDIRK \(s=3\) special case of the main conjecture.

\(-\) Mma also showed that PC\((-6,k)\) for \(k=0,1, 2\) and for all \(k\ge 3\) imply uniqueness. The first \(s\) point conditions have a special role. They imply a gap condition on \(a>0\): PC\((-6,0)\), PC\((-6,1)\) and PC\((-6,2)\) imply that \(a=\frac{1}{6}\) or \(a\ge \frac{1}{2}\) (and also some inequalities on \(c\) in terms of cubic polynomials of \(a\)). Of course, \(a=\frac{1}{6}\) implies \(c=\frac{1}{216}\). Now I’ve used Reduce\([...]\) with \(a\ge \frac{1}{2}\), the cubic inequalities and PC\((-6,k)\) (\(k\ge 3\)) to show that this time the solution set is empty. However, there is a subtle issue due to which I can not consider this as a second proof (yet). In Reduce\([...]\) I’ve used the ForAll\([k,k\ge 3, ...]\) construction, where \(k\ge 3\) is a real number and not only an integer. (Reduce could not cope with the problem if I restricted \(k\) to the Integers.) So in principle it may happen that the solution set is empty for real \(k\ge 3\), but non-empty for integer \(k\ge 3\). (Analogy: the function \(x(x-1)\) is \(\ge 0\) for integers but attains negative values for real \(x\) arguments.) If I use Reduce\([...]\) without separating \(k=0, 1, 2\) and only ForAll\([k,k\ge 0, ...]\), then no such \(a>0\) exists.

\(-\) Mma proved that PC\((-6,0)\), PC\((-6,1)\), PC\((-6,2)\) and IC\(([-\alpha,0],4)\) also imply \(a=\frac{1}{6}\) and \(c=\frac{1}{216}\) even if \(0< \alpha \ne 6\). I determined the set of these \(\alpha\)’s: uniqueness holds if and only if \(\alpha\in [0.103209..., 10.5836...]\), where the upper bound is a simple expression, but the lower bound is the root of a messy (but explicit) polynomial of degree 10.

\(-\) Mma proved (without assuming \(a>0\)) that PC\((-6,0)\), PC\((-6,1)\), PC\((-6,2)\) and PC\((-6,k)\) for all (real) \(k\ge 3\) imply either

\[ a=\frac{1}{6}, \quad\quad c=\frac{1}{216} \]

or

\[ a=-\frac{1}{6}, \quad\quad c=\frac{31}{216}. \]

This latter case corresponds to a stability function with \(x+6\) as a common factor. Apart from the first 2 derivatives (and from the 0th derivative), the sign of \(\psi^{(n)}\) on the whole \([-6,0]\) seems to be alternating (w.r.t. \(n\) even or odd), as in the case of completely/totally monotonic functions. (Notice that now \(\psi\) has a singularity at \(x=-6\), but in the definition of PC we’ve excluded the now singular term \(Q^{-k-s}\) from (*). Now \(Q(x)=1+\frac{x}{6}\).)

The gap conjecture (Mma proved this for \(s=3\), \(s=4\) and \(s=5\)). Suppose PC(\(-2s,k\)) for \(k=0, 1, \ldots, s-1\). Then \(a>0\) implies \(a\ge \frac{1}{2s}\), in other words, the \(x\)-coefficient in \(P\), \(a_1\le \frac{1}{2}\). Of course, in the conjectured optimal case we have equality.

As mentioned before, in the \(s=3\) case, the presence of another gap on the \(a\)-axis is implied by these conditions, namely there is no \(a\) in \((1/6,1/2)\). If \(s=4\), then there is no \(a\) in \((0,1/8)\cup (1/8,7/24)\). If \(s=5\), then there is no \(a\) in \((0,1/10)\cup (1/10,0.1835...)\), where the upper bound is the root of a particular cubic polynomial. Notice that the gap size (i.e. the length of the second interval component) is roughly halved in each \(s\)-step.

The \(s=4\) case. Apart from the “real \(k\) vs. integer \(k\)” subtlety in the ForAll quantifier (which I personally don’t believe to affect the validity of the proof), Mma proved the \(s=4\) case as well. Tests were run on two machines (under different Mma versions and operating systems; different \(a\)-intervals were scheduled on different machines “randomly”, and the proof was repeated one more time with different subdivisions). Some intervals were tested for hours in the beginning before finding a proper approach. Mma could prove there is a unique \(a=1/8\) in \((0,1/2]\) relatively fast (with corresponding \(P\)-coefficients \(c=1/128\) and \(d=1/4096\) as expected). It was also easy to prove no \(a\) exists if \(a\ge 5/2\). But, for example, Mma spent 55 sec in the interval \([204/100, 205/100]\), ten times as much time as in the neighboring interval. The interval \([16037992/10^7,16037993/10^7]\) also proved to be difficult. Of course there was no systematic way to know in advance which interval (and with what length) will require much computing time. As for the intervals \([207/100,208/100]\) and \([209/100, 210/100]\), I was able to prove with Mma that there is no such \(a\) only by using IC\(([-8,0],7)\) and IC\(([-8,0],6)\), respectively, in reasonable time. No point conditions helped here within a few minutes. At some point, the LC worked quite fast, but for some other \(a\) values, it could not prove uniqueness/no such \(a\) exists.

Let us return to the \(s=3\) case. By factoring out the positive \((-Q\prime)^k=a^k\) term, we see that the sign of (*) is determined by the expression

\[ \sum_{m=0}^{\min(k,3)} \frac{1}{m!}\binom{2+k-m}{2} \ P^{(m)}(x)\ \left(\frac{1-ax}{a}\right)^m. \]

A possible goal would be to prove (manually), for example, that if we have \[ \sum_{m=0}^{\min(k,3)} \frac{1}{m!}\binom{2+k-m}{2} \ P^{(m)}(-6)\ \left(\frac{1+6a}{a}\right)^m\ge 0 \] for all \(k\ge 0\) with some \(a>0\) and real \(c\), then \(P^{(m)}(-6)=0\), if \(m=0,1, 2\), and \(P^{(3)}(-6)=\frac{1}{36}\).

There are of course several relations between the numbers \(P^{(m)}(-6)\) (depending on \(a\) and \(c\), in a non-linear and in a linear way, respectively). The crucial thing now is to discover which of these relations would be useful also in the general case. Notice the similarity between the structure of the above expression and the linear combination in the Lemma on this page.

Notice also that the non-linear relation \[\frac{1}{2}\left(a_1^2-\frac{(1-a_1)^2}{s}\right)\ge a_2\] on the same page between \(a_1\) and \(a_2\) is now automatically satisfied by the stucture of \(P\).

A related conjecture, discovered again by pure chance, now due to a sign error is as follows.

It is more general, because it assumes no relations between the numbers \(p_m\), \(q\) and \(\widetilde{q}\). (Of course, I had in mind the connection \(p_m=\frac{P^{(m)}(x)}{m!}\), \(q=Q(x)\) and \(\widetilde{q}=Q\prime=-a\), for some \(x\in [-2s,0]\).) This conjecture seems much easier: Mma proved it quite fast for \(1\le s\le 10\), and I proved it manually for \(s=2\).

Conjecture. Fix any \(q\ne 0\), \(\widetilde{q}>0\) and \(p_m\in\mathbb{R}\) (\(m=0,1,\ldots,s+1\)), and suppose that for each \(k=0,1,\ldots, s+1\) we have \[ \sum_{m=0}^{\min(k,s)} \binom{s-1+k-m}{s-1}\cdot p_m\cdot q^m\cdot (-\widetilde{q})^{k-m}\ge 0. \] Then \(p_0\ge 0\) and for \(m=1,2,\ldots, s+1\), we have \(p_m=p_0\binom{s}{m}\cdot\left(\frac{\widetilde{q}}{q}\right)^m\), moreover, the above sum is equal to \(p_0\), if \(k=0\), and \(0\) for \(k=1, 2, \ldots, s+1\).

Unfortunately, this conjecture can not be used in our case, because, due to the presence of the factor \((-1)^{k-m}\) in (*), we would have \(\widetilde{q}<0\). Nevertheless, the similarity between the above conjecture and, for example, (some parts of the proof of) this Proposition is striking.