Моделирование лаба 1

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
def RK4(f):
return lambda t, y, dt: (
lambda dy1: (
lambda dy2: (
lambda dy3: (
lambda dy4: (dy1 + 2*dy2 + 2*dy3 + dy4)/6
)( dt * f( t + dt , y + dy3 ) )
)( dt * f( t + dt/2, y + dy2/2 ) )
)( dt * f( t + dt/2, y + dy1/2 ) )
)( dt * f( t , y ) )
#velocity (m per s)
v = 5
#acseleration (m per s^2)
g = 9.8
#mass kg
m = .1
#angle
from math import pi
alpha = pi/4
#air paramethr
beta = .025
from math import copysign
from math import tan
from math import sin
from math import cos
def haliley(v, alpha):
if alpha <= 0 or alpha >= pi/2:
print("Bad angle")
return .0
return (tan(alpha)*2*(v**2)*cos(alpha)**2)/g
def newton(v, alpha, beta, m):
if alpha <= 0 or alpha >= pi/2:
print("Bad angle")
return .0
def fx(t, v):
return -(beta*v**2)/m
def fy(t, v):
return -copysign(1, v)*beta*v**2/m-g
du, dw = RK4(fx), RK4(fy)
dt = 0.0001
u, w = v*cos(alpha), v*sin(alpha)
t = dt
x_old = y_old = 0
x, y = u*dt, w*dt
while (y > 0):
t += dt
u_new, w_new = u + du(t, u, dt), w + dw(t, w, dt)
x_old = x
x += (u_new + u) * dt / 2
y_old = y
y += (w_new + w) * dt / 2
u, w = u_new, w_new
if y_old == 0:
return x_old
return (x_old*(y_old - y) - y_old*(x_old-x))/(y_old - y)
print("distance haliley = %2.6f " % haliley(v, alpha))
print("distance newton = %2.6f " % newton(v, alpha, beta, m))
print("error = %2.6f " % (haliley(v, alpha) - newton(v, alpha, beta, m)))
a = 0.25*0.1
print a