{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Introduction\n",
"\n",
"Your goal is to assess the trade-offs of stability, accuracy, and run-time of a least two methods applied to two ODE systems. \n",
"\n",
"The first ODE system is two proteins degrading. The ODEs are:\n",
"\n",
"\\begin{eqnarray*}\n",
"\\frac{dy_1}{dt}&=&-\\alpha y_1\\\\\n",
"\\frac{dy_2}{dt}&=&-\\beta y_2\\\\\n",
"\\end{eqnarray*}\n",
"\n",
"where $\\alpha=0.1$ and $\\beta=0.2$.\n",
"\n",
"The second ODE system is the Lotka-Volterra model of predator-prey populations and their interactions. You can check out the Wikipedia entry to learn more, but what you need to know here is that it is a 2-ODE system and that, with the right parameter values, it leads to oscillations. The ODEs are:\n",
"\n",
"\\begin{align*}\n",
"\\frac{dx}{dt} &= \\alpha x -\\beta x y\\\\\n",
"\\frac{dy}{dt} &= -\\gamma y + \\delta x y\\\\\n",
"\\end{align*}\n",
"\n",
"where $\\alpha=0.25$, $\\beta=0.01$, $\\gamma=1$, and $\\delta=0.01$.\n",
"\n",
"\n",
"## The Project\n",
"\n",
"Consider three methods -- the forward Euler method, the explicit trapezoidal method, and an adaptive step-size forward Euler method (that uses the trapezoidal method to estimate the local error).\n",
"\n",
"Compare the performance of each of your three methods on the degradation and predator-prey systems. In your analysis, you will want to consider both efficiency (how long does it take the method to run?) and accuracy (how close is the method to the \"true\" solution?). For the first two methods, be sure to consider the effect of step-size on accuracy and speed. For the third, be sure to analyze the step-sizes chosen and the effect of error tolerance on the step size. \n",
"\n",
"\n",
"## Extensions\n",
"\n",
"This project will be more interesting if you tackle one or more of these extensions. Feel free to develop your own extension. And, as always, a particularly thorough analysis is smiled upon.\n",
"\n",
"\n",
"- Analyze the van der Pol oscillator in addition to the models above. The equations for it are:\n",
"\\begin{align*}\n",
"\\frac{dy_1}{dt} &= y_2\\\\\n",
"\\frac{dy_2}{dt} &= \\mu (1-y_1^2)y_2 - y_1.\\\\\n",
"\\end{align*}\n",
"When $\\mu=50$, the period of oscillation is long (nearly 80) and the system is stiff. When $\\mu=1$, the period is short (approximately 6) and the system is not stiff. We expect implicit methods to perform better on the stiff system. Demonstrate that it is better to use the implicit method when $\\mu=50$, but that explicit methods are perfectly adequate when $\\mu=1$.\n",
"
- Implement additional explicit Runge-Kutta methods. In particular, it would be nice to see the results of a fourth order method.\n",
"
- Implement Adams-Bashforth method (this should be covered in Lecture 15).\n",
"
- Implement the Backward Euler method and compare its results to those of your other methods. (Note that backward euler is an implicit method, which means the solution at the new time step appears on both sides of the equation and, therefore, requires additional computation. For example, you will need to use Newton's method to approximately solve a nonlinear system of algebraic equations. If you chose to use Backward Euler, consult Stephanie before writing code.)\n",
"\n",
"
\n",
"\n",
"## Notes\n",
"\n",
"### Timing\n",
"To time your method, use the %timeit magic command. You use it either for a single line of code, e.g. \n",
"\n",
"\n",
"%timeit (t,y) = forwardEuler( models.protein_degradation, np.arange(0,50,0.01), np.array([1,1]), params)\n",
"\n",
"\n",
"or to time an entire cell, e.g.\n",
"\n",
"%%timeit\n",
"(t,y) = forwardEuler( models.protein_degradation, np.arange(0,50,0.01), np.array([1,1]), params) \n",
"\n",
"\n",
"If you use it in cell mode, then be sure to put just the solving code in the cell. Put any plotting code in a different cell, because you don't want to include that.\n",
"\n",
"### Accuracy\n",
"To compute the accuracy of your method, compare its output to the output of scipy.integrate.solve_ivp (which will stand in for the \"true solution\"). Set its relative tolerance so that it is \"tight\" (i.e. you will tolerate very little error; rtol=1e-8). Error calculations should be performed by comparing the solutions to $y$ at the same time steps, but you will likely be comparing the true solution to outputs from methods using lots of different time steps. The easiest way to handle this is to call scipy.integrate.solve_ivp for each unique set of timesteps you have. For example, when you use the adaptive step size function, you do not know ahead of time which timesteps you will have, so run that first, and use the output to set the t_eval parameter for solve_ivp:\n",
"\n",
"(t,y,steps) = forwardEulerAdaptiveStepSize( models.protein_degradation, 0, 50, \n",
" np.array([1,1]), params, 1e-4)\n",
"sol = scipy.integrate.solve_ivp(fun=models.protein_degradation, t_span=(0,50), y0=(1,1), \n",
" method='LSODA', t_eval=t, args=[params], rtol=1e-8, atol=1e-8 )\n",
"\n",
" \n",
"To compute the error, use one of the methods we described in Lecture 11 (first, compute the error at each time point, then summarize it over time)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
\n",
"\n",
"Begin your code and report here. Organize your report in the manner that best allows you to answer all the questions. Integrate your extensions with your main project.\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# Start code and write-up here."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
\n",
"\n",
"Acknowledgements\n",
"\n",
"*Replace this text with links to external resources and thanks to people you worked with.*"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.5"
}
},
"nbformat": 4,
"nbformat_minor": 4
}