```python
import numpy as np
import matplotlib.pyplot as plt
from sympy import symbols, simplify, lambdify, dsolve, Eq, Function, expand
import sympy
one = sympy.Rational(1)
from BSeries import bs
from scipy.integrate import solve_ivp
from nodepy import rk, ivp

h = sympy.Symbol('h')
```

As described in [this notebook](https://nbviewer.jupyter.org/gist/ketch/de14e33acabce66fd0bda2b979b5d16f) and [this paper](https://epubs.siam.org/doi/abs/10.1137/19M1290346), when the explicit 2-stage midpoint Runge-Kutta method is applied to the nonlinear oscillator problem

$$
\begin{bmatrix} p \\ q \end{bmatrix} = \frac{1}{p^2 + q^2}\begin{bmatrix} -q \\ p \end{bmatrix}
$$

the numerical solution energy $E=u^2+v^2$ is constant regardless of the step size.  Here we use the method of modified equations to get some insight into this behavior.

First, we set up the right hand side of the ODE system:


```python
p, q = symbols('p,q')
u = [p,q]
f = np.array([-u[1]/(u[0]**2+u[1]**2), u[0]/(u[0]**2+u[1]**2)])
IC = [1.,0.]
simplify(f)
```

and the midpoint Runge-Kutta method coefficients:


```python
# Runge's Method
A = np.array([[0,0],[one/2,0]])
b = np.array([0,one])
```

# Modified equation

Now we generate the modified equation.  This is a differential equation that is satisfied exactly by the numerical solution.  In principle it is in an infinite series (in the step size $h$), so we must truncate it at some order.


```python
series = bs.modified_equation(u, f, A, b, order=5)
series = simplify(series)
```


```python
series
```

As expected for a 2nd-order method, the modified equations contain no terms of order $h$.  Notice that in this case, there are also no terms of order $h^3$.  In fact, because of the symmetry of this method, only even powers of $h$ will appear in the modified equations.

## Solutions of the modified equation

Next we will solve the modified equations directly and compare with the exact solution.  We consider the modified equations truncated to different orders in $h$.


```python
dt = 1.
T = 20
N=1000
f = simplify(np.array([term.series(h,0,1).removeO() for term in series]))

def solve_truncated_modified_equations(order,dt):
    f = simplify(np.array([term.series(h,0,order+1).removeO() for term in series]))
    f_ = lambdify([p,q,h],f)
    
    def f_p_vec(t,u,h=dt):
        return f_(*u,h)

    soln = solve_ivp(f_p_vec,[0,T],IC,t_eval=np.linspace(0,T,N),rtol=1.e-12,atol=1.e-12,method='RK45')

    return soln.t, soln.y

tt = []
yy = []
for order in range(5):
    t, y = solve_truncated_modified_equations(order,dt=dt)
    tt.append(t)
    yy.append(y)
    
rk2 = rk.ExplicitRungeKuttaMethod(A,b)

f_ex = lambdify([p,q],f)
f_ex(0.,1.)

def f_vec(t,u):
    return f_ex(*u)


myivp = ivp.IVP(f=f_vec,u0=np.array(IC),T=T)

t_rk2, y = rk2(myivp,dt=dt)
y = np.array(y)
y_rk2 = y[:,0]

plt.figure(figsize=(16,12))

plt.plot(t_rk2,y_rk2,'o')
for i in [0,2,4]:
    plt.plot(tt[i],yy[i][0,:],'--')

plt.legend(['RK2']+['$O(h^'+str(p)+')$' for p in [0,2,4]],fontsize=20)
```

As expected, we see that including more terms in the modified equations yields a solution that is accurate to longer times.  But what is remarkable is that all solutions of the modified equations (like the numerical solution from the RK method itself) seem to be indeed periodic.  This suggests that the truncated modified equations are energy-conserving at every order.  Let's check this by looking more closely at the modified equations.

## Structure of the modified equation

First, we generate the modified equations to a higher order than before; this may take a few minutes if you are running the notebook yourself.  Then we extract the numerator and denominator of each of the series.


```python
series = bs.modified_equation(u, f, A, b, order=7)
series = simplify(series)
series
```


```python
rhs_p=simplify(series[0])
numerator_p = rhs_p.as_numer_denom()[0]
denominator_p = sympy.factor(rhs_p.as_numer_denom()[1])
rhs_q=simplify(series[1])
numerator_q = rhs_q.as_numer_denom()[0]
denominator_q = sympy.factor(rhs_q.as_numer_denom()[1])
```

Next, we loop over the terms (by powers of $h$) and simplify (symbolically) the ratios:


```python
Numer_p = sympy.Poly(numerator_p,h)
coeffs = Numer_p.all_coeffs()[::-1]
series_p = 0
for j, coeff in enumerate(coeffs):
    series_p += h**j*sympy.factor(coeff)/denominator_p
series_p
```

These are the terms in the right-hand side for $p'(t)$; we see that all the odd orders of $h$ vanish identically.  Meanwhile, the even order terms have a simple structure that is obvious except for the values of the coefficients appearing in each denominator.

This is just to double-check that the manipulations above actually gave us back the correct right-hand side:


```python
simplify(series_p - rhs_p)
```

We can simplify the right-hand side series for $q$ in the same way, and find a complementary structure:


```python
Numer_q = sympy.Poly(numerator_q,h)
coeffs = Numer_q.all_coeffs()[::-1]
series_q = 0
for j, coeff in enumerate(coeffs):
    series_q += h**j*sympy.factor(coeff)/denominator_q
series_q
```

We conjecture that the modified equation for this numerical solution has the form

$$
\begin{bmatrix} p'(t) \\ q'(t) \end{bmatrix} = \sum_{j=0}^\infty \alpha_j \frac{(ih)^j}{(p^2+q^2)^{j+1}} \begin{bmatrix} -q \\ p \end{bmatrix},
$$

where $\alpha_j=0$ for every odd value of $j$.

Since each term of the series on the RHS is energy-conservative, the full modified equation is energy-conservative.

Note also that this means the modified equation has the form
$$
    \tilde{u}'(t) = f(\tilde{u}) P(h/\|u\|^2_2)
$$
where $P$ is a polynomial and $u'(t)=f(u)$ is the original ODE system.

**Open problem 1**: Use the theory of B-series to prove the above conjecture, possibly also giving a formula for the coefficients $\alpha_j$.

**Open problem 2**: Find another combination of ODE system and explicit Runge-Kutta method (or other B-series integrator) that yields a modified equation with this kind of structure?  I.e., find another example of unconditionally stable explicit integration?

Note that for unconditional stability, is not necessary that each term in the modified equation be orthogonal to the vector $[p,q]$; it only needs to have non-positive inner product with that vector.

# Modifying integrator

Complementary to the above approach, we can instead derive a modified system of ODEs such that when the explicit midpoint method is applied to the modified system it gives the exact solution of the original oscillator equation.  It may be revealing to examine the structure of this *modifying integrator*.

For some other examples of modifying integrators, see [this notebook](cac0b451425f6e76afb533d61adbca8c).


```python
series = bs.modifying_integrator(u, f, A, b, order=5)
simplify(series)
```


```python
rhs_p=simplify(series[0])
numerator_p = rhs_p.as_numer_denom()[0]
denominator_p = sympy.factor(rhs_p.as_numer_denom()[1])
rhs_q=simplify(series[1])
numerator_q = rhs_q.as_numer_denom()[0]
denominator_q = sympy.factor(rhs_q.as_numer_denom()[1])
```


```python
Numer_p = sympy.Poly(numerator_p,h)
coeffs = Numer_p.all_coeffs()[::-1]
series_p = 0
for j, coeff in enumerate(coeffs):
    series_p += h**j*sympy.factor(coeff)/denominator_p
series_p
```


```python
Numer_q = sympy.Poly(numerator_q,h)
coeffs = Numer_q.all_coeffs()[::-1]
series_q = 0
for j, coeff in enumerate(coeffs):
    series_q += h**j*sympy.factor(coeff)/denominator_q
series_q
```

We see a similar structure again.  Namely, the modifying integrator system seems to take the form

$$
    u'(t) = f(u) \hat{P}(h/\|u\|^2_2)
$$


though the polynomial $\hat{P}$ has different coefficients compared to $P$ that appeared above in the modified equation.

We might easily conjecture that all terms take the form given above.  However, if we go ahead and compute more terms, we find a surprise:


```python
series = bs.modifying_integrator(u, f, A, b, order=6)
simplify(series)
```


```python
rhs_p=simplify(series[0])
numerator_p = rhs_p.as_numer_denom()[0]
denominator_p = sympy.factor(rhs_p.as_numer_denom()[1])
rhs_q=simplify(series[1])
numerator_q = rhs_q.as_numer_denom()[0]
denominator_q = sympy.factor(rhs_q.as_numer_denom()[1])
```


```python
Numer_p = sympy.Poly(numerator_p,h)
coeffs = Numer_p.all_coeffs()[::-1]
series_p = 0
for j, coeff in enumerate(coeffs):
    series_p += h**j*sympy.factor(coeff)/denominator_p
series_p
```


```python
Numer_q = sympy.Poly(numerator_q,h)
coeffs = Numer_q.all_coeffs()[::-1]
series_q = 0
for j, coeff in enumerate(coeffs):
    series_q += h**j*sympy.factor(coeff)/denominator_q
series_q
```

Surprisingly, we get a term of order $h^5$, and it is proportional to $u$!

If we go further, we have a similar $h^7$ term:


```python
series = bs.modifying_integrator(u, f, A, b, order=8)
simplify(series)
```


```python
rhs_p=simplify(series[0])
numerator_p = rhs_p.as_numer_denom()[0]
denominator_p = sympy.factor(rhs_p.as_numer_denom()[1])
rhs_q=simplify(series[1])
numerator_q = rhs_q.as_numer_denom()[0]
denominator_q = sympy.factor(rhs_q.as_numer_denom()[1])
```


```python
Numer_p = sympy.Poly(numerator_p,h)
coeffs = Numer_p.all_coeffs()[::-1]
series_p = 0
for j, coeff in enumerate(coeffs):
    series_p += h**j*sympy.factor(coeff)/denominator_p
series_p
```


```python
Numer_q = sympy.Poly(numerator_q,h)
coeffs = Numer_q.all_coeffs()[::-1]
series_q = 0
for j, coeff in enumerate(coeffs):
    series_q += h**j*sympy.factor(coeff)/denominator_q
series_q
```

e81dee00febc9a44e580bf417732d056edfaaf97
ipynb
Jupyter Notebook
examples/Analysis of an explicit energy-conserving discretization.ipynb
ranocha/BSeries
de42bf03280ef7860ae6ade1d20212906b37a3d4
BSD-2-Clause
2021-09-08T08:17:53.000Z
2021-09-18T21:40:28.000Z
examples/Analysis of an explicit energy-conserving discretization.ipynb
ranocha/BSeries
de42bf03280ef7860ae6ade1d20212906b37a3d4
BSD-2-Clause
2021-09-19T03:05:35.000Z
2021-09-19T04:42:49.000Z
examples/Analysis of an explicit energy-conserving discretization.ipynb
ranocha/BSeries
de42bf03280ef7860ae6ade1d20212906b37a3d4
BSD-2-Clause
2021-09-18T16:06:24.000Z
2021-09-18T16:06:24.000Z
Qwen/Qwen-72B
1. YES
2. YES