In [1]:
# RUN ME FIRST!!!
from IPython.core.display import HTML
def css_styling():
    styles = open("custom.css", "r").read()
    return HTML(styles)
css_styling()
Out[1]:
/* From Lorena Barba's AeroPython course: https://github.com/barbagroup/AeroPython */
In [2]:
# we'll turn off inline plotting this time so that plots open in separate windows
# then you can easily save them by clicking the save button...
%matplotlib
Using matplotlib backend: MacOSX

Homework 1

Problems:

  1. For $\Delta t =1$, $y_0=1$, and $y' = y + 1$, use Euler's method by hand to find $y_{4}$.

  2. For $\Delta t =1$, $y_0=1$, and $y' = y + 1$, use the Trapezoidal method by hand to find $y_{4}$.

  3. For this problem you should use/modify the code below in the next cell block.

    1. Modify the code to solve $y' = 2y - y^{2}$ for Nt=10000 and tf=10. For y0=1 and y0=4 numerically solve the equation and plot the solutions. Why do the two solutions behave differently? (Make sure to turn in your two plots with your homework.)
    2. For y0=1 find the true solution to the preceding ODE. (Hint, separate and integrate with a partial fraction decomposition.)
    3. For y0=1, solve the ODE by Euler's method for Nt = 128, 256, 512, 1024, 2048, and 4094. For each numerical solution calculate the maximum error with the true solution you found in part b) as $$ E_{Nt} = \max_{n=0\dots Nt} \left| y_{n} - y(t_{n}) \right|. $$ Plot $E_{Nt}$ vs Nt on a loglog scale. What is the approximate slope of the curve? (Hint, you might find it useful in Python to store two vectors:
      Nt = np.array([128,256,512,1024,2048,4096])
      and a vector err of the corresponding values of $E_{Nt}$. You can estimate the slope of the best fit line using the np.polyfit command on np.log(Nt) and np.log(err).)
    4. What would you expect the slope to be? Why? (This part requires material we will cover 9/26 in class.)
In [1]:
"""
Euler's method routines for MA579 homework 1, problem 3.

I RECOMMEND YOU COPY THIS CODE INTO THE BLANK CODE BLOCK BELOW AND WORK THERE!

Created 9/16/15 by Prof. Samuel Isaacson (isaacson@math.bu.edu)

"""

import numpy as np
from matplotlib import pyplot as plt
import matplotlib as mpl


def eulersMethodFunction( y , t ):
    """
    Returns f(y,t) for use in the ODE y' = f(y,t).
    
    Parameters
    ----------
    y : scalar 
        Solution at time t.
    
    t: scalar
        Time to evaluate function at. 
    
    Returns
    -------
    f : scalar
        Value of the function f(y,t), the rhs of the ODE.       
    """
    
    f =  -y*y
    return f


def eulersMethod( f, y0, Nt, tf ):
    """
    Runs Euler's method.
    
    Parameters
    ----------
    f  : function
        Evaluates f(y,t) from the ODE y' = f(y,t).
    y0 : scalar
        Initial value of y at time t0.
    Nt : int
        Number of timesteps to take.
    tf : float
        Final time to solve for y(t) at.
    
    Returns
    -------
    y  : array 
        Solution values at the times t. 
    t  : array
        Times to solve for y(t) at.        
    """
    
    dt     = tf / Nt
    t      = np.linspace(0, Nt, Nt+1) * dt
    y      = np.zeros( Nt+1 )
    y[0]   = y0
    
    
    for  n in range(0, Nt):
        y[n+1] = y[n] + dt * f(y[n], t[n])

    return y,t
        

# This remainder of the code actually runs everything:

# PARAMETERS:
y0 = 1.                        # initial condition
Nt = 10000                     # number of time steps to take
tf = 10.                       # final time to run to

#f = lambda y, t: -y*y        # alternative lambda function way to define f, try it out
                              # by uncommenting it and then commenting the line below!
f = eulersMethodFunction      # standard function way to define f (as above)


# RUN EULER'S METHOD, saving y and t
y,t = eulersMethod( f, y0, Nt, tf )


# PLOTTING:

# pretty up plots
mpl.rcParams.update({'font.size': 12})  # this sets the default font size in plots
mpl.rcParams['axes.labelsize']  = 14    # make the xlabel/ylabel text size a bit bigger
mpl.rcParams['lines.linewidth'] = 2     # change the default linewidth in plots
    
# plot the solution
plt.figure()
plt.plot(t, y, 'k')
plt.xlabel('$t$')
plt.ylabel('$y(t)$')

plt.show()
In [2]:
# CODE UP YOUR SOLUTION TO PROBLEM 3 HERE

Continuing problems:

  1. In this problem we will modify our HIV model to incorporate latently infected cells. Let $L(t)$ denote the amount of latently infected T-cells at time $t$. We modify our pre-treatment model to \begin{align*} \frac{d T^{*}}{dt} &= k V T + aL - \delta T^{*},\\ \frac{d L}{dt} &= f k V T - \mu L, \\ \frac{dV}{dt} &= N \delta T^{*} - c V. \end{align*}

    Here $a$ denotes the rate latent cells reactivate into regular infected cells, $f$ denotes the fraction for how much smaller the rate of latent cell production is than of infected T-cell production, and $\mu$ is the combined rate of loss of latent cells (incorporating both cell death and conversion to $T^{*}$ cells).
    1. Find the new steady-state relations for this model. (You may again assume that $T$ is a constant.)
    2. Write down a new set of ODEs for the post treatment model that now incorporate the latent cells, $L(t)$.
    3. Using the new steady-state relations, modify the code below to incorporate latent cells. (You may take $\mu = .111$, $a = .098$, $f = .1$.) Make sure to add an equation for the latent cell steady state value. Solve both the old model and the new model for up to 40 days (`tf = 40`). Turn in your semilog plots of the total viral RNA amount in each model along with your new code. What is the effect of incorporating latent cells in the model on treatment? Why?
In [3]:
"""
  eulerHIVModel.py (modify for problem 4)
  Euler's method based solver for Perelson HIV model with treatment.

  This code should essentially be what you implemented in lab01, i.e.
  the model without latent cells, but with treatment.
  
  I RECOMMEND YOU COPY THIS INTO A NEW CODE BLOCK BELOW AND THEN UPDATE FOR PROBLEM 3!
  
  Created by Prof. Samuel Isaacson (isaacson@math.bu.edu) 9/15/15
"""

import numpy as np
from matplotlib import pyplot as plt
import matplotlib as mpl

# model parameters (patient 104):
T     =  2000.    # constant level of uninfected cells
VSS   = 2*52000.  # steady-state plasma RNA copies before treatment = 2 * number of virions
                  # we multiply by two as each HIV viral capsid actually contains two RNAs...
k     = 10.       # reaction-rate for  V + T -> TI
delta = .5        # death rate for infected cells (TI)
c     = 3.68      # virus clearance rate

# numerical parameters
Tf    = 40.      # time to run simulation to (days)
Nt    = 160000    # number of time steps


####### Everything below here is determined by the parameters above #######

# addition parameters
N     = c / (k*T)          # number of virions produced by the death of an infected cell,
                           # determined from SS relations before treatment
TISS  = c*VSS/(N*delta)    # steady-state level of infected cells before treatment,
                           # determined from SS relations before treatment
TI    = np.zeros(Nt+1)     # TI(n+1)  = amount of infected cells at time t_n
VI    = np.zeros(Nt+1)     # VI(n+1)  = amount of infective virus particles at time t_n
VNI   = np.zeros(Nt+1)     # VNI(n+1) = amount of non-infective virus particles at time t_n
dt    = Tf / Nt            # number of time steps
t     = np.arange(0,Nt+1) * dt        # vector of times, t_n
                           
# initial conditions
TI[0]  = TISS
VI[0]  = VSS
VNI[0] = 0.                # no non-infective virus particles at time zero


# Euler's method
for n in range(0, Nt):
    TI[n+1]  = TI[n]  + dt * ( k*VI[n]*T - delta*TI[n] )
    VI[n+1]  = VI[n]  - dt * c * VI[n]
    VNI[n+1] = VNI[n] + dt * ( N*delta*TI[n] - c * VNI[n] )


# plot solutions

# pretty up plots
mpl.rcParams.update({'font.size': 12})  # this sets the default font size in plots
mpl.rcParams['axes.labelsize']  = 14    # make the xlabel/ylabel text size a bit bigger
mpl.rcParams['lines.linewidth'] = 2     # change the default linewidth in plots

# black = TI, blue = V = VI + VNI
plt.figure()           # create a new figure window
plt.plot(t, TI, 'k')   # plot TI
plt.xlabel('time (days)')
plt.ylabel('number of infected cells per ml, $T^{*}(t)$')

plt.figure()
plt.plot(t, VI+VNI, 'b')  # plot V
plt.xlabel('time (days)')
plt.ylabel('number of virus RNA copies per ml, $V(t)$')

# black = TI, blue = V = VI + VNI
plt.figure()           # create a new figure window
plt.semilogy(t, TI, 'k')   # plot TI
plt.xlabel('time (days)')
plt.ylabel('number of infected cells per ml, $T^{*}(t)$')

plt.figure()
plt.semilogy(t, VI+VNI, 'b')  # plot V
plt.xlabel('time (days)')
plt.ylabel('number of virus RNA copies per ml, $V(t)$')

plt.show()
In [31]:
# CODE THE LATENT CELL MODEL HERE USING THE CODE ABOVE AS A TEMPLATE

Continuing problems:

  1. Calculate the remainder $R_2$ to the Taylor Series expansion about $x=0$ of $\exp(x)$, $\exp(x) = 1 + x + \frac{x^2}{2!} +R_2$.
  2. Calculate the remainder $R_4$ to the Taylor Series expansion about $x=0$ of $\cos(x)$, $\cos(x) =1- \frac{x^2}{2!} + \frac{x^4}{4!} + R_4$.
In [ ]: