T1 = 0.1/h(T2 – T1)
T2 = 0.1/h(T1 – T2)
Antilock Braking System
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()
Tangent
y=2^x, y=3^x
the slope of the tangent line at x < 1, x > 1
the exponential function
y = e^x -> slope of the tange line is exactly 1.
e = lim n->∞(1 + 1/N)^N = 2.7
e^π = (1+ 3.1416/10,000)^10,000
(1 + x/4)^4 = 1 + x + 3/8x^2 + 1/16x^3 + 1/256x^4
e^x = 1 + x + x^2/2 + x^3/6 + x^4/24 + …
e^x+h = e^x * e^h
dy(x)/dx = -3y(x)
y(0) = 5
自然対数の底
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