In this notebook, we use the package [`bseries.py`](https://github.com/ketch/bseries) to derive modifying integrators for certain Runge-Kutta methods applied to first-order ODEs, and study how well the solution of the truncated equations approximates the exact solution.


```python
import numpy as np
from BSeries import trees, bs
import matplotlib.pyplot as plt
from nodepy import rk, ivp
from IPython.display import display, Math
import sympy
from sympy import symbols, simplify, lambdify, dsolve, Eq, Function
from sympy import Derivative as D
from sympy.abc import t
cf = trees.canonical_forest
one = sympy.Rational(1)
from sympy import sin
from scipy.integrate import solve_ivp
h = sympy.Symbol('h')
```

# Rigid body Euler equations

Here we reproduce the first example from CHV2007.  First we set up the system of 3 ODEs corresponding to the rigid body problem (CHV2007 eqn. (9)):


```python
from sympy.abc import alpha, beta, gamma
from sympy.abc import x, y, z
y = [symbols('y%d' % i) for i in range(1,4)]

u = [y[0],y[1],y[2]]
u
```


```python
f = np.array([alpha*y[1]*y[2],beta*y[2]*y[0],gamma*y[0]*y[1]])
for rhs in f:
    display(rhs)
```

We will derive a modifying integrator ODE system for the implicit midpoint method.


```python
# Implicit Midpoint Method
A = np.array([[one/2]])
b = np.array([one])
```


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

We can verify that the results match Eqn. (12) of CHV2007:


```python
s3 = sympy.Poly((series[0]-f[0])/(f[0]),h).coeffs()[1]
s3
```


```python
s5 = sympy.Poly((series[0]-f[0])/(f[0]),h).coeffs()[0]
simplify(s5 - 6*one/5*s3**2)
```

# Lotka-Volterra

Next we consider using the explicit Euler method to solve the Lotka-Volterra model:

$$
    p'(t) = (2-q)p \quad \quad q'(t)=(p-1)q.
$$


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

FE1 = rk.loadRKM('FE')

dt = 0.35
T = 15.
IC = [1.5,2.25]

f_ = lambdify([p,q],f)

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

myivp = ivp.IVP(f=f_vec,u0=np.array(IC),T=T)
t0, y0 = FE1(myivp,dt=dt)
y0 = np.array(y0)

f_ex = lambdify([p,q],f)

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

myivp = ivp.IVP(f=f_vec,u0=np.array(IC),T=T)
BS5 = rk.loadRKM('BS5')

t, y = BS5(myivp,errtol=1.e-10,dt=1.e-3)
y_exact = np.array(y)
```

The exact solution of this problem is periodic, but Euler's method produces an unstable trajectory.  Here we use an especially large timestep in order to more clearly illustrate what will follow, but the qualitative behavior is the same for any step size.


```python
plt.figure(figsize=(9,6))
plt.plot(y_exact[:,1],y_exact[:,0],'-k',lw=2)
plt.plot(y0[:,1],y0[:,0],'--b')
plt.xlim(0,4)
plt.ylim(0,2.5)
plt.legend(['Exact solution','Explicit Euler, dt=0.1'],fontsize=15);
```

Now we will derive a "modifying integrator".  What this means is that we will determine a perturbed RHS such that when Euler's method is applied to the perturbed RHS, the result is the exact solution to the original Lotka-Volterra system.  The perturbed system takes the form of a power series in $h$, and in order to compute with it we will truncate it at a certain order.  We can compare the accuracy (and qualitative behavior) obtained by truncating at different orders.


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


```python
ymod = []
for order in range(2,max_order+1):
    fs = simplify(np.array([term.series(h,0,order).removeO() for term in series]))
    f_ = lambdify([p,q,h],fs)

    def f_p_vec(t,u,h=dt):
        return f_(*u,h)

    myivp = ivp.IVP(f=f_p_vec,u0=np.array(IC),T=T)
    _, y1 = FE1(myivp,dt=dt)
    ymod.append(np.array(y1))
```


```python
colors = 'rgc'
plt.figure(figsize=(9,6))
plt.plot(y_exact[:,1],y_exact[:,0],'-k',lw=2)
plt.plot(y0[:,1],y0[:,0],'--b')
for j in range(max_order-1):
    plt.plot(ymod[j][:,1],ymod[j][:,0],'--'+colors[j],lw=3,alpha=0.5)
plt.xlim(0,4)
plt.ylim(0,2.5)
plt.legend(['Exact solution','Explicit Euler, dt='+str(dt),
            'EE with modified flow to O(h)',
           'EE with modified flow to $O(h^2)$',
           'EE with modified flow to $O(h^3)$'],fontsize=15);
```

We see that if we include one additional term, the resulting trajectory still grows, while with two additional terms the solution appears to be dissipative.  With each additional term the solution gets closer to the exact solution of the original problem, and with three added terms it is hard to see the difference between them at this scale.

08c0821dd4b442a297cfed3f1f0e15ed7846b8d6
ipynb
Jupyter Notebook
examples/Modifying integrators.ipynb
ranocha/BSeries
de42bf03280ef7860ae6ade1d20212906b37a3d4
BSD-2-Clause
2021-09-08T08:17:53.000Z
2021-09-18T21:40:28.000Z
examples/Modifying integrators.ipynb
ranocha/BSeries
de42bf03280ef7860ae6ade1d20212906b37a3d4
BSD-2-Clause
2021-09-19T03:05:35.000Z
2021-09-19T04:42:49.000Z
examples/Modifying integrators.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