Category: Python

Python 2.7.12 (v2.7.12:d33e0cf91556, Jun 27 2016, 15:19:22) [MSC v.1500 32 bit (Intel)] on win32
Type "copyright", "credits" or "license()" for more information.
>>> print(3+1)
4
>>> 3 + 1
4
>>> print(1 + 2 + 3 * 3)
12
>>> print((1 + 2 + 3) * 3)
18
>>> print(3**2)
9
>>> print(9 % 2)
1
>>> print(15 // 4)
3
>>> print(16 // 4)
4
>>> print(-5 // 4)
-2
>>> print((23 + 32 + 64)/ 3)
39

integer and floats

>>> print(3/4)
0
>>> print(16/4)
4
>>> 3 + 2.5
5.5
>>> 213.13
213.13
>>> 341.
341.0
>>> int(49.7)
49
>>> int(16/4)
4
>>> float(3520+3239)
6759.0
>>> print(0.1)
0.1
>>> print(0.1 + 0.1 + 0.1)
0.3

float e.g.
Length of a fish you caught, in metres
Length of time it took to catch the first fish, in hours

4.445e8 is equal to 4.445 * 10 ** 8 which is equal to 444500000.0.

# The current volume of a water reservoir (in cubic metres)
reservoir_volume = 4.445e8
# The amount of rainfall from a storm (in cubic metres)
rainfall = 5e6

# decrease the rainfall variable by 10% to account for runoff
rainfall *= 0.9
# add the rainfall variable to the reservoir_volume variable
reservoir_volume += rainfall
# increase reservoir_volume by 5% to account for stormwater that flows
reservoir_volume *= 1.05 
# into the reservoir in the days following the storm

# decrease reservoir_volume by 5% to account for evaporation
reservoir_volume *= 0.95
# subtract 2.5e5 cubic metres from reservoir_volume to account for water
# that's piped to arid regions.
reservoir_volume -= 2.5e5
# print the new value of the reservoir_volume variable
print(reservoir_volume)

Two types of friction
static: moves when F > μs・N
kinetic: F = μk * N

wheel slip
rolling, rocked

speed = 120 km/h
μ=1.0: 3.4s
μ=0.7: 4.9s

F = m*a
u*gravitational force(mg)
a = μ*g

t/ 120km/h = 1 / μg
t = 120km/h / μ*g

Computing the coefficient of Friction
s = 1 – w/v
F(s) = μ(s)*mqcG
V= -F(s)/mqc
F = ma
w = F(s)/mew – B

import math
from **** import *

h = 0.01
mass_quarter_car = 250.
mass_effective_wheel = 20.
g = 9.81

end_time = 5.
num_steps = int(end_time / h)

w = numpy.zeros(num_steps + 1)
v = numpy.zeros(num_steps + 1)
x = numpy.zeros(num_steps + 1)
times = h * numpy.array(range(num_steps + 1))

@show_plot(7, 7)
def plot_me():
	axes_x = matplotlib.pyplot.subplot(411)
	axes_v = matplotlib.pyplot.subplot(412)
	axes_w = matplotlib.pyplot.subplot(413)
	awes_s = matplotlib.pyplot.subplot(414)

	def friction_coeff(slip):
		return 1.1 * (1. - math.exp(-20. * slip)) - 0.4 * slip

	def wheel_slip():
		b_values = numpy.arange(70., 190.1, 30.)
		for b in b_values:
			x[0] = 0.

			for step in range(num_steps):
				if v[step] < 0.01:
					break
				s = max(0., 1. - w[step] / v[step])

				w[step + 1] = max(0., w[step + 1])

			axes_x.plot(times[:step], x[:step])
			axes_v.plot(times[:step], v[:step])
			axes_w.plot(times[:step], w[:step])
			axes_s.plot(times[:step], 1. - w[:step] / v[:step])
			p = int((0.35 + 0.4 * (b - b_values[0])/ (b_values[-1] - b_values[0])) * num_steps)
			axes_x.annotate(b, (times[p], x[p]),
					xytext = (-30, -30), textcoords = 'offset point',
					arrowprops = dict(arrowstyle = '-', connectionstyle = 'arc3, '))

			p = int((0.35 + 0.4 * (v - b_values[0]) / (b_values[-1] - b_values[0])) * num_steps)
			axes_x.annotate(b, (times[p], x[p]),
				xytext = (-30, -30), textcoords - 'offset points',
				arrowprops = dict(arrowstyle = '-', connectionstyle = 'arc3, rad = 0.2', shrinkB = 0.))

		return x, v, w

    axes_x.set_ylabel('Position\nin m', multialignment = 'center')
    axes_v.set_ylabel('Car velocity\nin m/s', multialignment = 'center')
    axes_w.set_ylabel('Wheel velocity\nin m/s', multialignment = 'center')
    axes_s.set_ylabel('Wheel\nslip', multialignment = 'center')
    axes_s.set_xlabel('Time in s')
    axes_s.set_ylim(0., 1.)

    return wheel_slip()

x, v, w = plot_me()

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

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

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