Ch21.3 Systems of ODEs

• We will apply Euler and RK4 methods to first order systems.
  
• Example of a \( 2 \times 2 \) first order linear system:

\[ \begin{aligned} y_1'(x) &= a_1(x)y_1(x) + b_1(x)y_2(x), \,\, y_1(x_0)=y_{10} \\ y_2'(x) &= a_2(x)y_1(x) + b_2(x)y_2(x), \,\, y_2(x_0)=y_{20} \end{aligned} \]

More simply:

\[ \begin{aligned} y_1' &= a_1y_1 + b_1y_2 = f_1(x,y_1,y_2) , \,\, y_1(x_0)=y_{10}\\ y_2' &= a_2y_2 + b_2y_2 = f_2(x,y_1,y_2), \,\, y_2(x_0)=y_{20} \end{aligned} \]

Vector form:

\[ \mathbf{y}' = \mathbf{f}(x,\mathbf{y}), \,\, \mathbf{y}(x_0)=\mathbf{y}_0 \]

It takes guts to be an organ donor.

• Math 466 Methods III: Modeling and Numerical Methods.
• The system below models logistic growth with interactions:

\[ \begin{aligned} y_1' &= a_1(1-y_1)y_1 + b_1y_1y_2, \,\, y_1(x_0)=y_{10}\\ y_2' &= a_2(1-y_2)y_2 + b_2y_1y_2, \,\, y_2(x_0)=y_{20} \end{aligned} \]

• Force-balance considerations lead to 2nd order IVP:

\[ y'' + y' + y = 0, \,\, y(0)= 1, y'(0) = 1 \]

• Analytical solution: Assume \( y(x) = e^{rx} \). Then

\[ \small{ \begin{aligned} r^2 + r + 1 & = 0 \Rightarrow r = -1/2 \pm \sqrt{3}/2 i \\ y(x) &= e^{-x/2}\cos\left(\sqrt{3}/2 \right)x + \sqrt{3} e^{-x/2}\sin\left(\sqrt{3}/2 \right)x \end{aligned} } \]

• We are interested in numerical method of solving

\[ y'' + y' + y = 0, \,\, y(0)= 1, y'(0) = 1 \]

• Let \( y_1 = y \) and \( y_2 = y' \). Then

\[ \begin{aligned} y_1' &= f_1(y_2) = y_2, & y_1(0) = 1 \\ y_2' &= f_2(y_1,y_2) = - y_2-y_1, & y_2(0)=1 \end{aligned} \]

• Transform a 2nd order IVP into system of two 1st order IVPs:

\[ \small{ y'' + a(x)y' + b(x)y = 0, \,\, y(x_0)= y_{10}, y'(x_0) = y_{20} } \]

• Let \( y_1 = y \) and \( y_2 = y' \). Then

\[ \begin{aligned} y_1' &= y_2, & y_1(x_0)=y_{10}\\ y_2' &= - ay_2 - by_1, & y_2(x_0)=y_{20} \end{aligned} \]

• Expressing right sides as \( f_1 \) and \( f_2 \) respectively:

\[ \begin{aligned} y_1' &= f_1(y_2), & y_1(x_0)=y_{10}\\ y_2' &= f_2(x,y_1,y_2), & y_2(x_0)=y_{20} \end{aligned} \]

• Transform 2nd order IVP into system of two 1st order IVPs:

\[ y'' -2y' -3y = 0, \,\, y(0)= 1, y'(0) = 4 \]

• Let \( y_1 = y \) and \( y_2 = y' \). Then

• Transform 2nd order IVP into system of two 1st order IVPs:

\[ y'' + y' + y = 0, \,\, y(0)= 1, y'(0) = 1 \]

• Let \( y_1 = y \) and \( y_2 = y' \). Then

import numpy as np
import matplotlib.pyplot as plt 
from scipy.linalg import norm

def eulersys(N):   #N = number of nodes in [a,b]
    f1 = lambda x: x              #f1 from ODE
    f2 = lambda x,y: -y-x         #f2 from ODE
    g = lambda x: np.exp(-x/2)*np.cos(np.sqrt(3)/2*x) + np.sqrt(3)*np.exp(-x/2)*np.sin(np.sqrt(3)/2*x) #Exact solution 

a = 0                #Left endpoint of [a,b]
b = 10               #right endpoint of [a,b] 
h = (b-a)/N          #Compute step size h  

x = np.zeros(N+1)    #Initialize x values 
y1 = np.zeros(N+1)   #Initialize y1 values
y2 = np.zeros(N+1)   #Initialize y2 values
gn = np.zeros(N+1)   #Initialize gn values 
en = np.zeros(N+1)   #Initialize en values

x[0] = 0             #Initial x value 
y1[0] = 1            #Initial y1 = y(x0) value 
y2[0] = 1            #Initial y2 = y'(x0) value 

#Euler Method
for n in range(0,N): 
    x[n+1] = x[n] + h    
    y1[n+1] = y1[n] + h*f1(y2[n])       
    y2[n+1] = y2[n] + h*f2(y1[n],y2[n])  
    gn[n+1] = g(x[n+1])                  
    en[n+1] = g(x[n+1]) - y1[n+1]          

NormE = norm(en)        
print('Norm of Error Vector = ', NormE )
PRE = NormE/norm(gn)*100  
print('Percent Relative Error = ', PRE )

#Plot the numerical and analytical solutions
L = N*h  #[x0,x0+L] =x-axis interval for solution
xvalues = np.linspace(x0,x0+L,1000) 
plt.figure(figsize=(12,5.5))
plt.plot(xvalues,g(xvalues),'r', linewidth=3) 
plt.plot(x,y1,'b--', linewidth=3) 
plt.xlabel('x')
plt.legend(('Exact','Euler'),loc = 0)
plt.show()

#Runge-Kutta Method
for n in range(0,N): 
    x[n+1] = x[n] + h        
    a1 = h*f1(y2[n])   
    a2 = h*f2(y1[n],y2[n])   
    b1 = h*f1(y2[n]+0.5*a2)    
    b2 = h*f2(y1[n]+0.5*a1,y2[n]+0.5*a2)    
    c1 = #Homework     
    c2 = #Homework     
    d1 = #Homework        
    d2 = #Homework         
    y1[n+1] = y1[n] + (a1+2*b1+2*c1+d1)/6  
    y2[n+1] = y2[n] + (a2+2*b2+2*c2+d2)/6  
    gn[n+1] = g(x[n+1]) 
    en[n+1] = g(x[n+1]) - y1[n+1]          

• For both the Euler and RK4 methods, the error norm \( || \mathbf{E} || \) and percent relative error decreases as the number of nodes in \( [0,10] \) is increased from \( N=100 \) to \( N=1000 \).
• The Euler solution graph for \( N=100 \) follows the trend of exact solution, but with visible error.
  
• The Euler graph for \( N=1000 \) follows solution closely.
• The RK4 graph closely fits solution in both cases.
• The PRE values for Euler are substantially larger than RK4; see error ratios below.

  PRE Ratios   RK4/Euler
0    N = 100     62500.0
1   N = 1000  60000000.0

The following ODE arises in Math 260/360 (Ch5):

\[ \small{ y'' - xy = 0 } \]

Method of series \( y = \sum a_n x^n \) is used to find analytic solution.

\[ \scriptsize{ \begin{aligned} y_1(x) &= y_{10}\left[ 1 + \frac{1}{3!}x^3 + \frac{4}{6!}x^6 + \cdots \right] + y_{20}\left[ x + \frac{2}{4!}x^4 + \frac{10}{7!}x^8 + \cdots \right] \\ y_2(x) &= y_{10}\left[\frac{3}{3!}x^2 + \frac{24}{6!}x^5 + \cdots \right] + y_{20}\left[1 + \frac{8}{4!}x^3 + \frac{80}{7!}x^7 + \cdots \right] \end{aligned}} \]

We next seek to find numerical solution to Airy's IVP:

\[ y'' - xy = 0, \,\, y(0)= y_{10}, \,\, y'(0) = y_{20} \]

From the book example,

\[ \small{ y_1(0) = y_{10} = \frac{1}{3^{2/3}\Gamma(2/3)}, \,\, y_2(0)= y_{20} = -\frac{1}{3^{1/3}\Gamma(1/3)} } \]

Let \( y_1 = y \) and \( y_2 = y' \). Then

\[ \begin{aligned} y_1' &= f_1(y_2) = y_2, & y_1(0) = y_{10} \\ y_2' &= f_2(x,y_1) = x y_1, & y_2(0)= y_{20} \end{aligned} \]

import numpy as np
import matplotlib.pyplot as plt 
from scipy.linalg import norm
from scipy import special   
def rk4sys(N):   #N = number of nodes in [a,b]
    f1 = lambda x: x              #f1 from ODE
    f2 = lambda x,y: x*y         #f2 from ODE
    x0 = 0;
    y10 = 1/(3**(2/3)*special.gamma(2/3));
    y20 = -1/(3**(1/3)*special.gamma(1/3));

• The error norm \( \lVert \mathbf{E} \rVert \) and percent relative error values are very small, on the order of \( 10^{-10} \) and \( 10^{-8} \), respectively.
• The Desmos graph below shows the graphs of the \( 6^{th} \) degree partial sums for the analytical solutions \( y_1 \) and \( y_2 \). The plots of the numerical solution matches these graphs, and appears to do better near the endpoints than the partial sums.