e = 2.71828 18284 59045 23536 02874 71352
Euler’s number & the exponential function
5^1 = 5, 5^0 = 1, 5^-2= 1/25, 5^1/2= √5, 5^3*5^4=5^7, (5^3)^4=5^12
Infinity Time
x(t) = (x(t)^2 + 1)/2
with * (0) = 3
from *** import * end_time = 42. h = 0.0001 num_steps = int(end_time / h) times = h * numpy.array(range(num_steps) + 1) x = numpy.zeros(num_steps + 1) x[0] = 3. def forward_euler(): for step in range(num_steps): return x x = forward_euler() @show_plot def plot_me(): matplotlib.pyplot.plot(times, x) matplotlib.pyplot.show() plot_me()
Sources of Errors
-model, parameters, initial data, numerics
Logistic growth, Constant Harvesting
Exp growth, Amount of fish
from xxx import *
harvest_rates = [2e4, 5e4, 1e5, 2e5]
data = []
def logistic_growth():
maximum_growth_rate = 0.5
carrying_capacity = 2e6
end_time = 10.
h = 0.1
num_steps = in(end_time / h)
times = h * numpy.array(range(num_steps + 1))
fish = numpy.zeros(num_steps + 1)
fish[0] = 2e5
for harvest_rate in harvest_rates:
data.append(([time for time in times], [f for f in fish],
str(harvest_rate)))
return fish
fish = logistic_growth()
@show_plot
def plot_me():
fish_plots = []
for (times, fish, rate_label) in data:
fish_plots.append(matplotlib.pyplot.plot(times, fish, label=rate_label))
matplotlib.pyplot.legend([str(rate) for rate in harvest_rates], loc='upper right')
axes = matplotlib.pyplot.gca()
axes.set_xlabel('Time in years')
axes.set_ylabel('Amount of fish in tons')
plot_me()
Sir Backwards
R(t) = 1 / 5days I(t)
R2 = R1 + h 1/5daysI2
S2 = S1 + h(- 5 * 10^-9/ day*Person I1,S1)
Implicit solvers: use it or not?
Trapezoidal Stability
Cases – Moderate value of h, hk/2 = 0.1
hk/2 = 1
High value of h, e.g. hk/2 = g
Models with flows between different compartments
from xxx import *
h = 0.5
end_time = 60.
num_steps = int(end_time / h)
times = h * numpy.array(range(num_steps + 1))
def waning():
transmission_coeff = 5e-9
infectious_time = 5.
waning_time
s = numpy.zeros(num_steps + 1)
i = numpy.zeros(num_steps + 1)
r = numpy.zeros(num_steps + 1)
s[0] = 1e8 - 1e6 - 1e5
i[0] = 1e5
r[0] = 1e6
for step in range(num_steps):
s2i = h * transmission_coeff * s[step] * i[step]
i2r = h / infectious_time * i[step]
s[step + 1] = s[step] - s2i
i[step + 1] = i[step] + s2i - i2r
r[step + 1] = r[step] + i2r
return s, i, r
s, i, r = waning()
@show_plot
def plot_me():
s_plot = matplotlib.pyplot.plot(times, s, label = 'S')
i_plot = matplotlib.pyplot.plot(times, i, label = 'I')
r_plot = matplotlib.pyplot.plot(times, r, label = 'R')
matplotlib.pyplot.legend(('S', 'I', 'R'), loc = 'upper right')
axes = matplotlib.pyplot.gca()
axes.set_xlabel('Time in days')
axes.set_ylabel('Number of persons')
matplotlib.pyplot.xlim(xmin = 0.)
matplotlib.pyplot.ylim(ymin = 0.)
plot_me()
Contagion
A Contagious Disease
infectious, suspicious, Recovered
SIR Model
I(t)
R(t)= 1/5days * I(t)
A model for infection
S(t):number of suspicious person
I(t):number of infectious person
I(t) = 5 * 10^-9 / day*person I(t)*S(t)
SIR Model
S(t)= -5*10^-9/day*person I(t)S(t)
I(t) = 5 * 10^-9 / day*person I(t)*S(t)
R(t)= 1/5days * I(t)
System, Explicit
from xxx import *
h = 0.5
transmission_coeff = 5e-9
latency_time = 1.
infectious_time = 5.
end_time = 60.0
num_steps = int(end_time / h)
times = h * numpy.array(range(num_steps + 1))
def seir_model():
s = numpy.zeros(num_steps + 1)
e = numpy.zeros(num_steps + 1)
i = numpy.zeros(num_steps + 1)
r = numpy.zeros(num_steps + 1)
s[0] = 1e8 - 1e6 - 1e5
e[0] = 0.
i[0] = 1e5
r[0] = 1e6
for step in range(num_steps):
return s, e, i, r
s, e, i, r = seir_model()
@show_plot
def plot_me():
s_plot = matplotlib.pyplot.plot(times, s, label = 'S')
e_plot = matplotlib.pyplot.plot(times, e, label = 'E')
i_plot = matplotlib.pyplot.plot(times, i, label = 'I')
r_plot = matplotlib.pyplot.plot(times, r, label = 'R')
matplotlib.pyplot.legend(('S', 'E', 'I', 'R'), loc = 'upper right')
axes = matplotlib.pyplot.gca()
axes.set_xlabel('Time in days')
axes.set_ylabel('Number of persons')
matplotlib.pyplot.xlim(xmin = 0.)
matplotlib.pyplot.ylim(ymin = 0.)
x(t)=-k*(t), k=const>0
x(h)=x(0)+h(-k*(h))
Area in Phase Space
For simplicity: ID
x = v
v = F(x)/m
Errors on Different Time Scales
LTE, GTE
Long-term
Kinetic energy: 1/2ms ||Vs||^2
Potential enery: -G mems/||Xs||
def total_enery(): num_steps = 20000 x = numpy.zeros([num_steps + 1, 2]) v = numpy.zeros([num_steps + 1, 2]) enery = numpy.zeros(num_steps + 1) x[0, 0] = 15e6 x[0, 1] = 1e6 v[0, 0] = 2e3 v[0, 1] = 4e3 for step in range(num_steps): x[step + 1] = x[step] + h * v[step] v[step + 1] = v[step] + h * acceleration(x[step]) return x, energy
Circular Orbit
Orbital period:T
Vs
speed = ||Vs|| = 2πr/T
Acceleration = ||as|| = (2π)^2 r/t^2
Geostationary Orbit
T = 24h
Acceleration => (2π)^2 r/T^2 = 3√GmeT^2/(2π)^2 = 42,000 km
GTE
Error/h1 = E/(h2-h1) = Error = E h1/(h2-h1)
“Order” of a Solves
GTE = C*h order
x = sin(t)
Ch^2new = tolerance
Differential equation
xs = G * mE * -xs / |xs| + G * mm * Xm-xs/ |xm-xs|^3
non-linear
x(t)=sin(t)
x(t)=sin(x(t))
x(t)=cos(y(t))
t=(y(t))^2
differential equation that govern the motion
simple solution method
import math
from xxx import *
moon_distance = 384e6
def orbit():
num_steps = 50
x = numpy.zeros([num_steps + 1, 2])
for i in range(num_steps + 1):
angle = 2.*path.pi * i / num_steps
x[i, 0] = moon_distance * math.cos(angle)
x[i, 1] = moon_distance * math.sin(angle)
return x
x = orbit()
@show_plot
def plot_me():
matplotlib.pyplot.axis('equal')
matplotlib.pyplot.plot(x[:, 0], x[:, 1])
axes = matplotlib.pyplot.gca()
axes.set_xlabel('Longitudinal position in m')
axes.set_ylabel('Lateral position in m')
plot_me()
import math
from xxx import *
h = 0.1
g = 9.81
acceleration = numpy.array
initial_speed = 20. # m / s
@show_plot
def trajectory():
angles = numpy.linspace(20., 70., 6)
num_steps = 30
x = numpy.zeros([num_steps + 1, 2])
v = numpy.zeros([num_steps + 1, 2])
for angle in angles:
matplotlib.pyplot.plot(x[:, 0], x[:, 1])
matplotlib.pyplot.axis('equal')
axes = matplotlib.pyplot.gca()
axes.set_xlabel('Horizontal position in m')
axes.set_ylabel('Vertical position in m')
return x, v
trajectory()
Vector
vector is arrows
a-> = (4 3)
b-> = (2 1)
c-> = (1 -3)
b + c = sum of vector
||a|| = √4^2 + 3^2 = 5
Newton’s low of gravity
F = G m1*m2 / d^2
Gravity = 6.67*10^-11 Nm^2/kg^2
Mass of the moon
m2 = 1.62 m/s^2
R = 1.737*10^6m
A2 = 4A, A1
xM, dES, dMS
import numpy earth_mass = 5.97e24 moon_mass = 7.35e22 gravitational_constant = 6.67e-11 def acceleration(moon_position, spaceship_position): vector_to_moon = moon_position - spaceship_position vector_to_earth = - spaceship_position nonsense = numpy.linalg.norm(moon_position) * spaceship_position return nonsense