Assignment 5: Solve ODE and store in NetCDF file

After an ice storm, the doors of one apartment building are frozen shut, trapping everyone inside.  To make matters worse, some of the tenants suddenly become zombies.  The zombies go rampaging through the building, turning other apartment dwellers into zombies.

Fortunately, some people know how to kill zombies and are able to teach the other people in the apartment building.  The change of surviving a zombie encounter differs between zombie killers and non-zombie killers, but if a person does not survive the encounter, they become a zombie.

This process can be described by the Modified Zombie Apocalypse equations:

\(\displaystyle\frac{dx}{dt} = f(x,y,z) = -\beta x z - \delta x y\)

\(\displaystyle\frac{dy}{dt} = g(x,y,z) = -\gamma y z + \delta x y\)

\(\displaystyle\frac{dz}{dt}= h(x,y,z) = \beta xz + \gamma y z - \alpha yz\)

where

• \(x(t)\) is the uninfected population without zombie killing skills ("regular people")
• \(y(t)\) is the uninfected population with zombie killing skills ("zombie killers")
• \(z(t)\) is the zombie population

These are fractions of the original population of the apartment building, so that \(0\leq x\leq1\), \(0\leq y\leq1\), \(0\leq z\leq1\), \(x(0)+y(0)+z(0) = 1\).

The parameters of the model are constants which can be interpreted as follows:

• \(\alpha\) is the rate at which zombies are killed (by \(y\))
• \(\beta\) is the rate at which regular people are turned into zombies  
• \(\gamma\) is the rate at which zombie killers are turned into zombies
• \(\delta\) is the rate at which zombie killers teach regular people how to kill zombies.

(This model was inspired by Munz et al., Infectious Disease Modeling Research Progress, 2009.)

Your assignment is write a program 'mzasolve' that uses GSL's odeiv2 sublibrary to solve this system of equations and writes out the time series to a NetCDF file.

Details:

• Assume the values \(\alpha=3, \beta=2,\gamma=1,\delta=3/2\).
• Let initial conditions for y and z be given as command line parameters.  
• Let the filename of the output netcdf file also be given as a command line parameter.
• I.e., "./mzasolve 0.018 0.004 timeseries.nc" should solve the equations with $y(0)=0.018$ and $z(0)=0.004$ and write the time series to the file "timeseries.nc".
• Use an adaptive time step with rk8pd as the step type with a relative error of \(10^{-6}\)and an initial time step of \(10^{-3}\).
• The time series should be saved in a netcdf file as four 1d arrays containing the time values, x values, y values and z values.
• The time evolution should run until a steady state has been reached.  As a criterion to see if the system is in a steady state, you can use

$$ \frac{\sqrt{f^2+g^2+h^2}}{x+y+z} < 10^{-6}$$

While there is little opportunity for modularization in this assignment, we do expect the code to contain a number of functions, and the compilation of the code and at least one run case to be automated in a Makefile. We also expect the code to be under git version control with several meaningful commits, and that your code is properly commented.

The git repository only has to contain the code for the assignment and its Makefile, not the netcdf file. Upload a zip file with your git repo (as in previous assignments) by 11:59 PM on Friday March 3, 2023.