Simple Differential Equations

In this section we solve simple differential equations. These will serve as an example of the concept of an initial value problem and the methods of solution, as well as provide a fist look at the accuracy of the numerical solutions.

Radioactive Decay

Arguably the simplest differential equation is the one describing Radioactive Decay.

This equation is

$$\frac{ d n(t) }{ dt } = - \lambda n(t)$$

The function $n(t)$ gives the number of radioactive nuclei remaining in a sample at time $t$. It is derived via a probability argument. $\lambda$ is the probability per unit time for a nucleus to decay. This constant is fitted to decay data of the nucleus being described. The model describes the number of radioactive nuclei remaining in the sample at all time given $n(0)=No$, the number of radioactive nuclei initially at time $t=0$. The half life of a radioactive species is the time at which only 1/2 of the radioactive nuclei remain in the sample.

The exact solution of this equation is:

$$n(t) = N_0 e^{- \lambda t }$$

Here is how Maple exactly solve this equation: Radioactive_decay.mw

Exercise #1: Write a C code to solve this equation via the Euler and the Mid-Point methods. Compare these numerical solution to the exact solution. Here are the steps you may want to follow:

1. write a skeleton C code which establishes a time grid ( use $t_{min} = 0.0$, $t_{max} = 12$, and $dt = 0.2$ )
2. add the model definition and initial value of $n(t)$ ( use $\lambda = 0.2$ and $No = 1200$ )
3. add the statements to calculate and output the function at the first time grid point
4. add the statements to calculate and output the result the Euler step inside the time loop
5. plot the Euler solution
6. add statements to the code to calculate the exact solution
7. plot the Euler solution and the exact solution
8. add statements to the code to calculate and output the Mid-Point solution
9. plot this solution with the exact and the Euler solution

The code radiactive_decay.c compares both the Euler and Mid-Point solutions to the exact solution of this problem.

Aids Epidemics

Under epidemics conditions, the number of sick individuals, $n(t)$, in a population sample containing a total of $Np$ individuals is given as a crude first approximation by the equation

$$\frac{d n(t)} { dt } = \lambda n(t) ( N_p - n(t) )$$

where $\lambda$ is the probability per unit time that a healthy individual become sick when in contact with a sick individual.

Exercise #2: Solve this equation for $n(t)$ given $\lambda=0.0005$, $Np=1200$, and $n(0)=No=5$. Use a time range $t=[0,15]$ and $dt=0.3$. Solve the differential equation via the Euler and the Mid-Point methods. Plot these solutions simultaneously.

Exercise #3: The solution of a differential equation can be checked via direct replacement of the solution in the equation and checking whether the LHS equals the RHS. Do this for the Euler and Mid-Point solutions above. Here are the steps to follow:

1. store the Euler and Mid-Point solutions in arrays as well as the time grid (or read in the data from the piped output file in a separate program)
2. compute and output the LHS of the differential equation via the symmetric difference formula for both solutions (within a new loop)
3. compute and output the RHS of the differential equation
4. plot the results
5. what do you conclude about the accuracy of the solutions

The exact solution of this differential equation can be obtained. Here is a Aids_epidemic.mw illustrating how.

Fixed Point

The following differential equation has an interesting solution.

$$\frac{ d \theta(t) } {d t } = sin( \sqrt{2} \theta(t) ) - cos( \sqrt{3} \theta(t) )$$

Exercise #4: Numerically solve this differential equation

1. Use the Mid-Point method to get your solution using $dt=0.1$ and a time interval $t=[0,12]$
2. Run your code with $\theta(0)=0.55$ and $\theta(0)=0.45$, plotting $theta(t)$ as a function of time in the given interval
3. Explain the behavior of your solution