In this notebook, we use the package [`bseries.py`](https://github.com/ketch/bseries) to derive modified equations for certain Runge-Kutta methods applied to first-order ODEs, and study how well the solution of the modified equations approximates the numerical 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')
```

# Lotka-Volterra

Here we reproduce the example from p. 340 of the book *Geometric Numerical Integration* (Hairer, Lubich, & Wanner), using the explicit Euler method to solve the Lotka-Volterra model:

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

First we define the model:


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

Next, we load the coefficients of the method and generate the modified equations as a B-series:


```python
FE1 = rk.loadRKM('FE')

A = FE1.A
b = FE1.b

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

The numerical solution of the LV model by the explicit Euler method is the exact solution to a system of *modified differential equations*; this system can be expressed as a power series in the step size $h$.  Here we have derived the right had side of that system up to terms of order $h$.  Notice that if we drop the $O(h)$ terms then we have again the original LV system.

We can check that the $O(h)$ terms match what is given in HLW:


```python
-sympy.expand(simplify(series[0]+p*(q-2))*2/(h*p))
```


```python
-simplify(series[1]-q*(p-1))*2/(h*q)
```

Next, we'll solve the modified equations very accurately and compare the result with the numerical solution given by the explicit Euler method with step size $h=0.1$.


```python
dt = 0.1
T = 15.
IC = [1.5,2.25]

fs = simplify(np.array([term.series(h,0,2).removeO() for term in series]))
f_ = lambdify([p,q,h],fs)

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,1000),rtol=1.e-12,atol=1.e-12,method='RK45')

t1, y1 = soln.t, soln.y
```


```python
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)

t, y = FE1(myivp,dt=dt)
y = np.array(y)
```


```python
plt.figure(figsize=(9,6))
plt.plot(y[:,1],y[:,0],'o')
plt.plot(y1[1,:],y1[0,:],'--k')
plt.xlim(0,9)
plt.ylim(0,5.5)
plt.legend(['Explicit Euler, dt=0.1','Modified flow to O(h)'],fontsize=15)
```

The exact solution of the LV model is periodic, but Euler's method generates a solution with growing amplitude.  The modified equations accurately predict this.

Now we go to the next order.


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


```python
dt = 0.12
T = 14.5
IC = [1.,2.75]

fs = simplify(np.array([term.series(h,0,2).removeO() for term in series]))
f_ = lambdify([p,q,h],fs)
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,1000),rtol=1.e-12,atol=1.e-12,method='RK45')
t1, y1 = soln.t, soln.y

fs = simplify(np.array([term.series(h,0,3).removeO() for term in series]))
f_ = lambdify([p,q,h],fs)
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,1000),rtol=1.e-12,atol=1.e-12,method='RK45')
t2, y2 = soln.t, soln.y

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)

t, y = FE1(myivp,dt=dt)
y = np.array(y)
```


```python
plt.figure(figsize=(9,6))
plt.plot(y[:,1],y[:,0],'o')
plt.plot(y1[1,:],y1[0,:],'--')
plt.plot(y2[1,:],y2[0,:],'--k')
plt.xlim(0,9)
plt.ylim(0,5.5)
plt.legend(['Explicit Euler, dt=0.12','Modified flow to $O(h)$','Modified flow to $O(h^2)$'],fontsize=15);
```

Using a larger step size, we see that the 1st-order modified equations are not fully accurate, but by including the $O(h^2)$ terms we get much better accuracy at late times.

Let's keep going.


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


```python
dt = 0.2
T = 10.
IC = [1.,2.75]

fs = simplify(np.array([term.series(h,0,2).removeO() for term in series]))
f_ = lambdify([p,q,h],fs)
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,1000),rtol=1.e-12,atol=1.e-12,method='RK45')
t1, y1 = soln.t, soln.y

fs = simplify(np.array([term.series(h,0,3).removeO() for term in series]))
f_ = lambdify([p,q,h],fs)
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,1000),rtol=1.e-12,atol=1.e-12,method='RK45')
t2, y2 = soln.t, soln.y

fs = simplify(np.array([term.series(h,0,4).removeO() for term in series]))
f_ = lambdify([p,q,h],fs)
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,1000),rtol=1.e-12,atol=1.e-12,method='RK45')
t3, y3 = soln.t, soln.y



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)

t, y = FE1(myivp,dt=dt)
y = np.array(y)
```


```python
plt.figure(figsize=(9,6))
plt.plot(y[:,1],y[:,0],'o')
plt.plot(y1[1,:],y1[0,:],'--')
plt.plot(y2[1,:],y2[0,:],'--')
plt.plot(y3[1,:],y3[0,:],'--k')
plt.xlim(0,15)
plt.ylim(-0.5,6.5)
plt.legend(['Explicit Euler, dt='+str(dt),'Modified flow to $O(h)$','Modified flow to $O(h^2)$','Modified flow to $O(h^3)$'],fontsize=15)
```

Again, with a larger step size we see that additional terms are needed to obtain good accuracy at later times.


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


```python
dt = 0.1
T = 66.4
IC = [1.,2.01]
N = 3000
fs = simplify(np.array([term.series(h,0,2).removeO() for term in series]))
f_ = lambdify([p,q,h],fs)
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')
t1, y1 = soln.t, soln.y

fs = simplify(np.array([term.series(h,0,3).removeO() for term in series]))
f_ = lambdify([p,q,h],fs)
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')
t2, y2 = soln.t, soln.y

fs = simplify(np.array([term.series(h,0,4).removeO() for term in series]))
f_ = lambdify([p,q,h],fs)
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')
t3, y3 = soln.t, soln.y

fs = simplify(np.array([term.series(h,0,7).removeO() for term in series]))
f_ = lambdify([p,q,h],fs)
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')
t5, y5 = soln.t, soln.y


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)

t, y = FE1(myivp,dt=dt)
y = np.array(y)
```


```python
plt.figure(figsize=(9,6))
plt.plot(y[:,1],y[:,0],'o')
plt.plot(y1[1,:],y1[0,:],'--')
plt.plot(y2[1,:],y2[0,:],'--')
plt.plot(y3[1,:],y3[0,:],'--')
plt.plot(y5[1,:],y5[0,:],'--k')
plt.xlim(-0.5,18)
plt.ylim(-0.5,11.5)
plt.legend(['Explicit Euler, dt='+str(dt),'Modified flow to $O(h)$','Modified flow to $O(h^2)$',
            'Modified flow to $O(h^3)$','Modified flow to $O(h^6)$'],fontsize=15);
```

Here we have gone all the way up to the $O(h)^6$ terms and we continue to get improved accuracy for long times.

# Pendulum

Next we consider another simple first-order system of two equations that models a rigid frictionless pendulum (see e.g. p. 4 of HLW).


```python
f = np.array([-sin(u[1]),u[0]])
IC = [1.,0.]
simplify(f)
```

This time we'll consider a more accurate numerical method: a 3-stage, 3rd-order Runge-Kutta method.


```python
rk3 = rk.loadRKM('SSP33')
A = rk3.A
b = rk3.b

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

Notice that the modified equations (which we have derived up to order $h^5$) include no correction terms of order $h$ or $h^2$.  This is true because the method chosen is 3rd-order accurate.

Again, we compare a highly-accurate solution of the modified equations with the approximate solution of the original problem obtained using the Runge-Kutta method.


```python
dt = 1.05
T = 20
N=1000
```


```python
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(7):
    t, y = solve_truncated_modified_equations(order,dt=dt)
    tt.append(t)
    yy.append(y)
```


```python
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_rk3, y = rk3(myivp,dt=dt)
y = np.array(y)
y_rk3 = y[:,0]
```


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

plt.plot(t_rk3,y_rk3,'o')
for i in range(2,6):
    plt.plot(tt[i],yy[i][0,:],'--')

plt.legend(['RK3']+['$O(h^'+str(p)+')$' for p in range(2,6)],fontsize=20)
```

We can see that each successive correction gives a solution that is accurate to later times than the one previous.  Notice that in this case, although the exact solution is periodic, the numerical solution is gradually damped, and this behavior is captured by the more accurate versions of the modified equations.

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