# RUN ME FIRST!!!
from IPython.core.display import HTML
def css_styling():
styles = open("custom.css", "r").read()
return HTML(styles)
css_styling()
# 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
For $\Delta t =1$, $y_0=1$, and $y' = y + 1$, use Euler's method by hand to find $y_{4}$.
For $\Delta t =1$, $y_0=1$, and $y' = y + 1$, use the Trapezoidal method by hand to find $y_{4}$.
For this problem you should use/modify the code below in the next cell block.
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.)y0=1
find the true solution to the preceding ODE. (Hint, separate and integrate with a partial fraction decomposition.)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])
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)
.) """
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()
# CODE UP YOUR SOLUTION TO PROBLEM 3 HERE
"""
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()
# CODE THE LATENT CELL MODEL HERE USING THE CODE ABOVE AS A TEMPLATE