Category: Python

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()

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))

For simplicity: ID
x = v
v = F(x)/m

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

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

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 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

F = m * a(9.81m/s^2)

x(t) = v(t)
v(t) = a(t) = F/m

small time steps of size h
x(h) = x(0) + hv(0)
v(h) = v(0) + h F/m

from xxxplots import

def forward_euler():
	h = 0.1
	g = 9.81

	num_steps = 50

	t = numpy.zeros(num_steps + 1)
	x = numpy.zeros(num_steps + 1)
	v = numpy.zeros(num_steps + 1)

	for step in range(num_steps):
		t[step + 1] = h * (step + 1)
		x[step + 1] = x[step] + h * v[step]
		v[step + 1] = v[step] - h * g
	return t, x, v

t, x, v = forward_euler()

@show_plot
def plot_me():
	axes_height = matplotlib.pyplot.subplot(211)
    matplotlib.pyplot.plot(t, x)
    axes_velocity = matplotlib.pyplot.subplot(212)
    matplotlib.pyplot.plot(t, v)
    axes_height.set_ylabel('Height in m')
    axes_velocity.set_ylabel('Velocity in m/s')
    axes_velocity.set_xlabel('Time in s')
    # Uncomment the line below when running locally.
    # matplotlib.pyplot.show() 

plot_me()

import numpy
import matplotlib.pyplot

def forward_euler():
	h = 0.1
	g = 9.81

	num_steps = 50

	t = numpy.zeros(num_steps + 1)
	x = numpy.zeros(num_steps + 1)
	v = numpy.zeros(num_steps + 1)

	for step in range(num_steps):
		t[step + 1] = h * (step + 1)
		x[step + 1] = x[step] + h * v[step]
		v[step + 1] = v[step] - h * g
	return t, x, v

t, x, v = froward_euler()

k = 130
m = 7655

if k > 42:
    print 'Hi'
elif m > 43:
    print 'Hello'
else:
    print 'Bye'

======================== RESTART: C:/Python27/test.py ========================
Hi
>>> 

import numpy

n = numpy.zeros(5)
n[0] = 20.
n[4] = 99.

print n
print n[-1]

import numpy

p = numpy.zeros([2, 3])

p[0, 1] = 7.
p[1, 2] = 8.

print p
p = 1. + 2. * p

print p

def three_times(u):
    r = 2 * u
    return r + u

print three_times(1), three_times(7)

======================== RESTART: C:/Python27/test.py ========================
3 21
>>> 

import math
from xxx import *

def sin_cos():
    num_points = 50

    x = numpy.zeros(num_points)
    sin_x = numpy.zeros(num_points)
    cos_x = numpy.zeros(num_points)

    for i in range(num_points):
        x[i] = 2. * math.pi * i / (num_points - 1.)
        sin_x[i] = math.sin(x[i])
        cos_x[i] = math.cos(x[i])

    return x, sin_x, cos_x

x, sin_x, cos_x = sin_cos()

@show_plot
def plot_me():
    matplotlib_pyplot.plot(x, sin_x)
    matplotlib_pyplot.plot(x, cos_x)
plot_me()

import math

f = math.sqrt(4)
g = radius_of_earth * math.sin(math.pi / 180.0 * 23.4)
g = radius_of_earth * math.sin(math.pi / 180.0 * 23.4)

print f, g

======================== RESTART: C:/Python27/test.py ========================
2.0 2530229.21123
>>> 

for h in range(7):
    i = h**2
    print i

======================== RESTART: C:/Python27/test.py ========================
0
1
4
9
16
25
36

k = 42
m = 43

while k < 130 or m == 140: k += 1 m += k print k, m [/python]

======================== RESTART: C:/Python27/test.py ========================
130 7655
>>>