Update notes

Author: cpraveen

myst.yml

@@ -65,6 +65,7 @@ project:
       - file: ode/multistep.md
       - file: ode/fe_be_stability.md
       - file: ode/euler2.md
+      - file: ode/rk.md
       - file: ode/ode_stability.md
       - file: ode/stiff1.md
       - file: ode/stiff2.md

ode/euler.md

@@ -18,8 +18,7 @@ numbering:
 ```
 
 ```{code-cell}
-import numpy as np
-from matplotlib import pyplot as plt
+from pylab import *
 ```
 
 The following function implements Euler method.
@@ -31,8 +30,8 @@ $$
 ```{code-cell}
 def euler(f,t0,T,y0,h):
     N = int((T-t0)/h)
-    y = np.zeros(N)
-    t = np.zeros(N)
+    y = zeros(N)
+    t = zeros(N)
     y[0] = y0
     t[0] = t0
     for n in range(1,N):
@@ -64,29 +63,28 @@ def f(t,y):
     return y
 
 def yexact(t):
-    return np.exp(t)
+    return exp(t)
 ```
 
 We run the Euler method for decreasing values of step size.
 
 ```{code-cell}
 t0,y0,T = 0.0,1.0,5.0
 
-te = np.linspace(t0,T,100)
+te = linspace(t0,T,100)
 ye = yexact(te)
-plt.plot(te,ye,'--',label='Exact')
+plot(te,ye,'--',label='Exact')
 
 H = [0.2,0.1,0.05]
 
 for h in H:
     t,y = euler(f,t0,T,y0,h)
-    plt.plot(t,y,label="h="+str(h))
+    plot(t,y,label="h="+str(h))
 
-plt.legend(loc=2)
-plt.xlabel('t')
-plt.ylabel('y')
-plt.grid(True);
+legend(loc=2), xlabel('t'), ylabel('y'), grid(True);
 ```
+
+Even visually we observe that the numerical solution is converging towards the exact solution as $h \to 0$.
 :::
 
 +++
@@ -109,10 +107,10 @@ Right hand side function and exact solution
 
 ```{code-cell} 
 def f(t,y):
-    return -y + 2.0*np.exp(-t)*np.cos(2.0*t)
+    return -y + 2.0*exp(-t)*cos(2.0*t)
 
 def yexact(t):
-    return np.exp(-t)*np.sin(2.0*t)
+    return exp(-t)*sin(2.0*t)
 ```
 
 **Solve for a given h**
@@ -122,40 +120,36 @@ t0,y0 = 0.0,0.0
 T  = 10.0
 h  = 1.0/20.0
 t,y = euler(f,t0,T,y0,h)
-te = np.linspace(t0,T,100)
+te = linspace(t0,T,100)
 ye = yexact(te)
-plt.plot(t,y,te,ye,'--')
-plt.legend(('Numerical','Exact'))
-plt.xlabel('t')
-plt.ylabel('y')
-plt.grid(True)
-plt.title('Step size = ' + str(h));
+plot(t,y,te,ye,'--')
+legend(('Numerical','Exact'))
+xlabel('t'), ylabel('y'), grid(True)
+title('Step size = ' + str(h));
 ```
 
 **Solve for several h**
 
 Study the effect of decreasing step size. The error is plotted in log scale.
 
 ```{code-cell}
-hh = 0.1/2.0**np.arange(5)
-err= np.zeros(len(hh))
+hh = 0.1/2.0**arange(5)
+err= zeros(len(hh))
 for (i,h) in enumerate(hh):
     t,y = euler(f,t0,T,y0,h)
     ye = yexact(t)
-    err[i] = np.abs(y-ye).max()
-    plt.semilogy(t,np.abs(y-ye))
+    err[i] = abs(y-ye).max()
+    semilogy(t,abs(y-ye))
     
-plt.legend(hh)
-plt.xlabel('t')
-plt.ylabel('log(error)')
+legend(hh), xlabel('t'), ylabel('log(error)')
 
 for i in range(1,len(hh)):
-    print('%e %e %e'%(hh[i],err[i],np.log(err[i-1]/err[i])/np.log(2.0)))
+    print('%e %e %e'%(hh[i],err[i],log(err[i-1]/err[i])/log(2.0)))
 
 # plot error vs h
-plt.figure()
-plt.loglog(hh,err,'o-',label="Error")
-plt.loglog(hh,hh,'--',label="y=h")
-plt.legend(), plt.xlabel('h'), plt.ylabel('Error norm');
+figure()
+loglog(hh,err,'o-',label="Error")
+loglog(hh,hh,'--',label="y=h")
+legend(), xlabel('h'), ylabel('Error norm');
 ```
 :::

ode/rk.md

@@ -0,0 +1,95 @@
+---
+exports:
+  - format: pdf
+kernelspec:
+  display_name: Python 3 (ipykernel)
+  language: python
+  name: python3
+numbering:
+  code: false
+  equation: true
+  title: true
+  headings: true
+---
+
+# Runge-Kutta method
+
+```{include} ../math.md
+```
+
+```{code-cell}
+from pylab import *
+```
+
++++
+
+:::{prf:example} Harmonic oscillator
+We repeat the previous example with RK method.
+
+```{code-cell}
+y0    = 1.0
+omega = 4.0
+h     = 0.1
+
+A = matrix([[0.0,       1.0], 
+            [-omega**2, 0.0]])
+hA = h * A
+A2 = eye(2) + hA + (1.0/2.0)*hA**2
+A4 = eye(2) + hA + (1.0/2.0)*hA**2 + (1.0/6.0)*hA**3 + (1.0/24.0)*hA**4
+
+n = int(2*pi/h)
+y2, y4 = zeros((n,2)), zeros((n,2))
+y2[0,0] = y0
+y4[0,0] = y0 
+t = zeros(n)
+for i in range(1,n):
+    y2[i,:] = A2 @ y2[i-1,:]
+    y4[i,:] = A4 @ y4[i-1,:]
+    t[i]     = t[i-1] + h
+
+plot(t, cos(omega*t),'--',label="Exact")
+plot(t,y2[:,0],label="RK2")
+plot(t,y4[:,0],label="RK4")
+axis([0, 2*pi,-2,2])
+legend(), xlabel("t"), ylabel("$\\theta$");
+```
+
+The solution from RK2 is growing with time, no matter how small a time step we take. RK4 will be conditionally stable, if $h < 2\sqrt{2}/\omega \approx 0.707$.
+:::
+
+:::{prf:example} Stability regions of RK
+
+```{code-cell}
+r2 = lambda z: 1 + z + z**2/2
+r3 = lambda z: 1 + z + z**2/2 + z**3/6
+r4 = lambda z: 1 + z + z**2/2 + z**3/6 + z**4/24
+
+n = 100
+x = linspace(-3,0.5,n)
+y = linspace(-3,3,n)
+X,Y = meshgrid(x,y)
+Z2 = r2(X + 1j*Y)
+Z3 = r3(X + 1j*Y)
+Z4 = r4(X + 1j*Y)
+contour(X,Y,abs(Z2),levels=[1],colors="black")
+contour(X,Y,abs(Z3),levels=[1],colors="blue")
+contour(X,Y,abs(Z4),levels=[1],colors="red")
+xlabel('real(z)'), ylabel('imag(z)')
+title('Stability domain of RK2, RK3, RK4')
+grid(True), axis('equal');
+```
+
+Zoomed view around the imaginary axis
+
+```{code-cell}
+contour(X,Y,abs(Z2),levels=[1],colors="black")
+contour(X,Y,abs(Z3),levels=[1],colors="blue")
+contour(X,Y,abs(Z4),levels=[1],colors="red")
+axis([-0.3, 0.3, -3, 3])
+xlabel('real(z)'), ylabel('imag(z)')
+title('Stability domain of RK2, RK3, RK4')
+grid(True);
+```
+
+RK2 does not enclose any portion of the imaginary axis, it is tangential to it, while RK3 and RK4 enclose a part of the imaginary axis. In case of RK4, the  interval $[ -\ii 2\sqrt{2}, \ii 2\sqrt{2}]$ is inside the stability domain.
+:::

ode/shooting.md

@@ -15,6 +15,7 @@ kernelspec:
 
 +++
 
+:::{prf:example}
 Consider the BVP
 
 $$
@@ -33,11 +34,7 @@ $$
 y(x) = (e^x + e^{-x})^{-1}
 $$
 
-+++
-
-## Formulate as an initial value problem
-
-+++
+**Formulate as an initial value problem**
 
 $$
 y'' = -y + \frac{2 (y')^2}{y}, \qquad -1 < x < +1
@@ -55,11 +52,7 @@ $$
 \phi(s) = y(1;s) - (e + e^{-1})^{-1} = 0
 $$
 
-+++
-
-## Newton method
-
-+++
+**Newton method**
 
 To find the root of $\phi(s)$, we use Newton method with an initial guess $s_0$. Then the Newton method updates the guess by
 
@@ -84,32 +77,37 @@ The derivative of $\phi(s)$ is given by
 $$
 \frac{d}{ds}\phi(s) = \frac{\partial}{\partial s} y(1;s) = z_s(1)
 $$
+
 We can write an equation for $z_s(x)$
+
 $$
 z_s'' = \left[ -1 - 2 \left( \frac{y'}{y} \right)^2 \right] z_s + 4 \frac{y'}{y} z_s'
 $$
+
 with initial conditions
+
 $$
 z_s(-1) = 0, \qquad z_s'(-1) = 1
 $$
 
-+++
-
-## First order ODE system
-
-+++
+**First order ODE system**
 
 Define the vector
+
 $$
 u = \begin{bmatrix}
 u_1 \\ u_2 \\ u_3 \\ u_4 \end{bmatrix} = \begin{bmatrix}
 y \\ y' \\ z_s \\ z_s' \end{bmatrix}
 $$
+
 Then
+
 $$
 u_1' = u_2, \qquad u_2' = -u_1 + \frac{2 u_2^2}{u_1}, \qquad u_3' = u_4, \qquad u_4' = \left[ -1 - 2 \left( \frac{u_2}{u_1} \right)^2 \right] u_3 + 4 \frac{u_2}{u_1} u_4
 $$
+
 Hence we get the first order ODE system
+
 $$
 u' = f(u) = \begin{bmatrix}
 u_2 \\
@@ -118,7 +116,9 @@ u_4 \\
 \left[ -1 - 2 \left( \frac{u_2}{u_1} \right)^2 \right] u_3 + 4 \frac{u_2}{u_1} u_4
 \end{bmatrix}
 $$
+
 with initial condition
+
 $$
 u(-1) = \begin{bmatrix}
 (e+e^{-1})^{-1} \\
@@ -127,14 +127,14 @@ s \\
 1
 \end{bmatrix}
 $$
+
 Once we solve this IVP, we get
+
 $$
 \phi(s) = u_1(1) - (e + e^{-1})^{-1}, \qquad \frac{d}{ds}\phi(s) = z_s(1) = u_3(1)
 $$
 
-+++
-
-## Now we start coding
+**Now we start coding**
 
 ```{code-cell} 
 import numpy as np
@@ -223,3 +223,4 @@ plt.grid(True)
 plt.title('Final solution')
 plt.legend(("Numerical","Exact"));
 ```
+:::