# How to support or contradict a hypothesis on unconditional stability using numerical optimization

The main motivation behind my next question is that I think I derived a higher order numerical scheme for linear advection equation that is unconditionally stable using Von Neumann stability analysis.

In the stability analysis of the numerical scheme I have to find the maximal value of so called amplification factor $|S|$ where $S$ depends on 5 parameters, it means $S=S(x,y,f,c,d)$ and $$S= \left(2 \left(-c^2+c \cos (x) (c-d (2 f+1) \cos (y)+2 d f+d-1)+d \cos (y) (2 c f+c+d-1)-2 c d f-c d \sin (x) \sin (y)-c d-3 i c \sin (x)+c-d^2-3 i d \sin (y)+d+6\right)\right) \cdot \left(2 i c^2 \sin (x)-i c^2 \sin (2 x)+c^2-2 i c d \sin (x+y)+2 c d \cos (x+y)+2 i c d \sin (x)-2 c (c+d+5) \cos (x)+2 i c d \sin (y)-2 c d \cos (y)+2 c d+10 i c \sin (x)-2 i c \sin (2 x)+(c+2) c \cos (2 x)+8 c+2 i d^2 \sin (y)-i d^2 \sin (2 y)-2 d^2 \cos (y)+d^2 \cos (2 y)+d^2+10 i d \sin (y)-2 i d \sin (2 y)-10 d \cos (y)+2 d \cos (2 y)+8 d+12\right)^{-1}$$ Note that $i$ is imaginary number.

My hypothesis is that $|S| \le 1$ for $-\pi \le x \le \pi$, $-\pi \le y \le \pi$, $-1 \le f \le 0$ and $c \ge 0$ and $d \ge 0$.

Up to now I checked the values of $|S|$ for many tabulated values of input variables e.g. for $c \le 100$ and $d \le 100$. Moreover I used Mathematica to compute numerically the maximal values of $|S|$ with many starting values and several built-in optimization methods, I got always the maximal value of $|S|$ being $1$. Of course I used the numerical scheme to compute several examples of linear advetion equation with very large time steps and observed no instabilities.

I would appreciate any useful thoughts, helps or even your tries ;-) how to find the values of $(x,y,f,c,d)$ bounded by above constraints for which $|S|>1$. Or from the other view - what kind of methods and arguments would you use to support the hypothesis that $|S| \le 1$ ?

P.S. Here is the input for Mathematica

S=(2*(6 + c - c^2 + d - c*d - d^2 - 2*c*d*f +
d*(-1 + c + d + 2*c*f)*Cos[y] +
c*Cos[x]*(-1 + c + d + 2*d*f - d*(1 + 2*f)*Cos[y]) -
3*I*c*Sin[x] - 3*I*d*Sin[y] - c*d*Sin[x]*Sin[y]))/(12 + 8*c + c^2 + 8*d + 2*c*d + d^2 - 2*c*(5 + c + d)*Cos[x] +
c*(2 + c)*Cos[2*x] - 10*d*Cos[y] - 2*c*d*Cos[y] -
2*d^2*Cos[y] + 2*d*Cos[2*y] + d^2*Cos[2*y] +
2*c*d*Cos[x + y] + 10*I*c*Sin[x] + 2*I*c^2*Sin[x] +
2*I*c*d*Sin[x] - 2*I*c*Sin[2*x] - I*c^2*Sin[2*x] +
10*I*d*Sin[y] + 2*I*c*d*Sin[y] + 2*I*d^2*Sin[y] -
2*I*d*Sin[2*y] - I*d^2*Sin[2*y] - 2*I*c*d*Sin[x + y])


P.S.S. If for any reasons you need more insights, see my paper at http://dl.dropboxusercontent.com/u/386482/articles/frolkovic-algoritmy-corr.pdf where the scheme is given in (3.10) for $f=0$.

• I know cylindrical algebraic decomposition sometimes works, but this problem looks like it's too large. – Kirill Oct 16 '16 at 21:35
• I tried it directly, without coming up with anything clever, and QEPCAD failed after an hour with "Prime list exhausted." – Kirill Oct 16 '16 at 23:13
• If you could clean up the equation by clustering function terms, it might be easier to come up with estimates. For instance I can't tell where the inversion begins. Presumably you can also collect the trigonometric terms. They can always be estimated by $\pm 1,0$, whichever gives you the largest or smallest value with the right sign. – Deathbreath Nov 1 '16 at 16:06

In your expression for $S$, is $x = k_x \Delta x$ and $y = k_y \Delta y$ for some "wave number" vector $\mathbf{k} = [k_x \quad k_y]^T$ with $\Delta x$ and $\Delta y$ being grid spacings in $x$ and $y$ directions?

I assumed that $x = \pi$ and $y = \pi$ and obtained the following for $S$: $$-{\frac {2\,cdp+{c}^{2}+{d}^{2}-c-d-3}{{c}^{2}+2\,cd+{d}^{2}+5\,c+5\,d +3}}$$ which satisfies $|S| \le 1$. The quantity $p$ is equal to $2f +1$.

The values $x=\pi$ and $y=\pi$, of course, correspond to the minimum resolvable wavelengths in each direction.

I played with a few other values of $x$ and $y$ as well and in all those cases, $|S| \le 1$.

Here is an example how I have tested the hypothesis that $|S| \le 1$ for given constraints. I have run the following code in Mathematica with many many values of variables nstart and nend and the method $M$:

(* amplification parameter *)
S=
(2*(6+c-c^2+d-c*d-d^2-2*c*d*f+d*(-1+c+d+2*c*f)*Cos[y]+c*Cos[x]*(-1+c+d+2*d*f-d*(1+2*f)*Cos[y])-3*I*c*Sin[x]-3*I*d*Sin[y]-c*d*Sin[x]*Sin[y]))
/
(12+8*c+c^2+8*d+2*c*d+d^2-2*c*(5+c+d)*Cos[x]+c*(2+c)*Cos[2*x]-10*d*Cos[y]-2*c*d*Cos[y]-2*d^2*Cos[y]+2*d*Cos[2*y]+d^2*Cos[2*y]+2*c*d*Cos[x+y]+10*I*c*Sin[x]+2*I*c^2*Sin[x]+2*I*c*d*Sin[x]-2*I*c*Sin[2*x]-I*c^2*Sin[2*x]+10*I*d*Sin[y]+2*I*c*d*Sin[y]+2*I*d^2*Sin[y]-2*I*d*Sin[2*y]-I*d^2*Sin[2*y]-2*I*c*d*Sin[x+y]);

(* numerical global maximum search *)
nstart=1;nend=10;

For[n=nstart,n<=nend,n++,
Print[NMaximize[{Abs[S],-Pi<=x<=Pi, -Pi<= y<= Pi, 0<= c <=n, 0<=d<=n, -1<=f<=0},
{x,y,c,d,f}, Method->{M,"SearchPoints"->20,RandomSeed->n}]]]

(* example of output *)
{1.,{x->3.87721*10^-7,y->-3.14159,c->0.996877,d->0.,f->-0.998816}}
{1.,{x->-3.14159,y->2.19618*10^-8,c->0.,d->1.71028,f->-0.941516}}
{1.,{x->-3.14159,y->-3.14159,c->0.,d->0.,f->-1.}}
{1.,{x->3.14159,y->3.14159,c->0.,d->0.,f->-1.}}
{1.,{x->4.71742*10^-7,y->3.14159,c->3.3804,d->0.,f->-0.997559}}
{1.,{x->-3.14159,y->-3.14159,c->0.,d->0.,f->-1.}}
{1.,{x->-3.14159,y->-2.4762*10^-7,c->0.,d->1.34533,f->-0.922121}}
{1.,{x->-0.000101745,y->-0.0000189795,c->3.48753,d->4.20828,f->-0.84976}
{1.,{x->3.14159,y->1.0698*10^-6,c->0.,d->8.98382,f->-0.63621}}
{1.,{x->0.000268716,y->-0.0000310003,c->0.593874,d->3.95438,f->-0.832499}}