How to solve a differential system using Runge-Kutta 2nd order method?

Question

I have question to solve using RK2. Your help is highly appreciated.

In a reversible chemical reaction concentrations $A$ and $B$ given in $\rm \text{g-mol}/L$ of two compound, change at rates given by $$ \begin{align} \frac{dA}{dt} &= -k_1 A^2 + k_2 B\\ \frac{dB}{dt} &= 2k_1 A^2 - k_2 B^2 \end{align} $$ with $A(0) = 1$, $B(0) = 0$, and $k_1 = 0.03$, $k_2 = 0.02$. Use RK2 method to find B at $t= 1,2$. Show how to apply an alternate method for the same purpose.

Answer 1

We have the system:

$$ \begin{align} \frac{dA}{dt} &= -k_1 A^2 + k_2 B = -0.03 A^2 + 0.02 B\\ \frac{dB}{dt} &= 2k_1 A^2 - k_2 B^2 = 2 (0.03) A^2 - 0.02 B^2 \end{align} $$

$$A(0) = 1, B(0) = 0$$

Use the Runge-Kutta Order 2, RK2, method to find $B$ at $t = 1,2$.

We arrive at the iterates:

• $t = 0, (A, B) = (1., 0.)$
• $t = 0.1, (A, B) = (0.997015, 0.005982)$
• $t = 0.2, (A, B) = (0.99406, 0.0119283)$
• $\ldots$
• $t = 1.0, (A, B) = (0.971451, 0.0582517)$
• $\ldots$
• $t = 2.0, (A, B) = (0.945618 , 0.113204)$

You can fill in the details for RK2 and also use a second method to verify these results.

Answer 2

The method can be implemented in general fashion (in python using numpy) as

def RK2(f,u,t):
     uout = np.zeros((len(t),)+u.shape)
     uout[0] = u;
     for k in range(len(t)-1):
         h = t[k+1]-t[k]
         k1 = f(uout[k],t[k])*h
         k2 = f(uout[k]+0.5*k1, t[k]+0.5*h)*h
         uout[k+1]=uout[k]+k2
     return uout

and called for the given problem as

def solve(k1,k2,A,B,M):
     def derivs(u,t): 
         A,B = u; 
         return np.array([ 
             -k1*A**2+k2*B, 
             2*k1*A**2-k2*B**2 
         ])

     u0 = np.array([A, B])
     t  = np.arange(0,2.05,1./M); 
     u  = RK2(derivs, u0, t)
     for k in range(1,3):
         print "   %5.2f  |  %s"%(t[k*M],u[k*M]);

Calling it with 3 different step sizes $h=1, 0.1, 0.01$ via

solve(0.03, 0.02, 1.0, 0.0, 1)
solve(0.03, 0.02, 1.0, 0.0, 10)
solve(0.03, 0.02, 1.0, 0.0, 100)

gives the results

1.00  |  [ 0.97149325  0.0581955 ]
2.00  |  [ 0.94569454  0.11310238]

1.00  |  [ 0.97145105  0.05825173]
2.00  |  [ 0.9456183   0.1132042 ]

1.00  |  [ 0.97145066  0.05825226]
2.00  |  [ 0.94561759  0.11320516]

indicating good results already with step size $h=1$.