```
%pylab inline
import numpy as np
import matplotlib.pyplot as plt
import time
```

In this problem we'll study the efficiency of various approaches to solving the initial value problem  
$$\begin{align} u_1'(t) & = -K_1 u_1 \\\\  
u_2'(t) & = K_1 u_1 -K_2 u_2 \\\\
u_3'(t) & = K_2 u_2\end{align}$$  
$$u(0) = [1,0,0].$$  
This corresponds to the decay process  
$$A\to^{K_1} B \to^{K_2} C $$  
To begin, we define a Python function that evaluates the right-hand-side of this system.  


```
def f(t,u):
    K1=1
    K2=2
    du=np.zeros([3])
    du[0]=-K1*u[0]
    du[1]= K1*u[0] - K2*u[1]
    du[2]=           K2*u[1]
    return du
```

The scipy.integrate.ode package contains interfaces to various numerical ODE solvers, similar to those available in MATLAB. Take a moment to read the help on this package. The scipy.integrate.odeint package contains similar functionality.  


```
from scipy.integrate import ode
```

We'll use the 'dopri5' integrator, which is an embedded Runge-Kutta pair (the same one used by MATLAB's ode45 command). The code below sets up and solves the problem. Note that the parameter dt_output below, does not determine the step size used by the integrator; it is merely the interval between outputs that we'll use for plotting. The step size is chosen to achieve specified error tolerances; in this case, we have specified $10^{-10}$ for both absolute and relative errors.  


```
t0 = 0.                  # Initial time
u0 = np.array([1.,0.,0.])# Initial values
tfinal = 4.              # Final time
r = ode(f).set_integrator('dopri5',atol=1.e-10,rtol=1.e-10) # Dormand-Prince RK4(5) method
r.set_initial_value(u0, t0)
dt_output=0.02           # Interval between output for plotting
N=tfinal/dt_output       # Number of output times
tt=np.zeros(N+1); tt[0]=t0  # Output times
uu=np.zeros([3,N+1]);       # Output values
uu[:,0]=u0
i=0
while r.successful() and r.t<tfinal:
    i=i+1
    uu[:,i]=r.integrate(r.t+dt_output)
    tt[i]=r.t
```

Finally, we plot the solution. Take some time to observe what is happening and make sure that it makes sense to you.


```
plt.clf()
plt.plot(tt,uu[0,:],tt,uu[1,:],tt,uu[2,:]);
plt.legend(['A','B','C']);
```

#Exercise 1
(a) Now consider the decay process  
$$A\to^{K_1} D \to^{K_3} B \to^{K_2} C.$$  
Modify the function f() above to return the right-hand-side for this system by adding $u_4=[D]$ with the reaction rate $K_3=3$ and then solve the system using the initial value $u_4=0$. Plot all 4 components of the result and make sure your results are reasonable.  


```

```

(b) Suppose $K_3$ is much larger than $K_1$ and $K_2$. Then as $A$ is converted to $D$, it decays almost instantly into $B$. In this case we would expect that $u_4(t)$ will always be very small (though nonzero for $t&gt;0$), while $u_j(t)$ for $j=1,2,3$ will be nearly identical to what would be obtained by solving the first system of 3 components (at the top of the worksheet) with the same reaction rates $K_1$ and $K_2$. Test this by using $K_3=1000$ and solving the system with 4 components.  


```

```

(c) Try solving the system with $K_3=1000$ and different tolerances, from $10^{-2}$ down to $10^{-10}$. Compare the time required to compute the solution (you may find it useful to use the python package time and the function time.time()). Compare this with the time taken to solve the original system (of three equations) with the same tolerances. You should observe that it is much slower to solve the system that includes $D$, even though the solution is essentially the same. You should also see that the cost of computing the solution doesn't increase much as the tolerance is reduced. Explain these two observations. <span style="color: #ff0000;">Hint: in order to get useful timing figures, set dt_output=tfinal. That way no time will be spent copying the output and you will measure just the time to compute the solution.</span>


```

```

#Exercise 2

(continuation of exercise 1)  
(a) Plot the computed solution from part (c) of the previous exercise, using the tolerances $10^{-2}$ and $10^{-4}$, and comment on what you observe.  
(b) Now set the tolerances to $10^{-6}$ and vary $K_3$ from 500 up to 2000. You should observe that the time to compute the solution grows linearly with $K_3$. Explain why you would expect this to be true (rather than being roughly constant, or growing at some other rate, such as quadratically in $K_3$). How long would it take to compute a solution for $K_3=10^7$?  
(c) Repeat part (b) but use the scipy.integrate.odeint function, which has an adaptive BDF integrator, in place of dopri5. Explain why the computational time is much smaller and now roughly constant for large $K_3$. Also try $K_3=10^7$.


```

```

4b00e501d3ba6dfee1ad8c4e937eb3d0b7f488ad
ipynb
Jupyter Notebook
ipython/HW_4_part_1.ipynb
bfdeandrade/Finite-Difference
97e93f394051a70b6aa2c26fa266952ae821deb2
CC-BY-2.0
2015-02-05T23:13:07.000Z
2022-01-05T03:09:04.000Z
ipython/HW_4_part_1.ipynb
MuriloHMoreira/finite-difference-course
97e93f394051a70b6aa2c26fa266952ae821deb2
CC-BY-2.0
ipython/HW_4_part_1.ipynb
MuriloHMoreira/finite-difference-course
97e93f394051a70b6aa2c26fa266952ae821deb2
CC-BY-2.0
2015-02-16T17:36:38.000Z
2021-12-11T07:13:03.000Z
Qwen/Qwen-72B
1. YES
2. YES