2018-04-22-175105.ipynb

Kernel: Python 3 (Ubuntu Linux)

%matplotlib -- inline 
import matplotlib.pyplot as plt
import numpy as np
import sympy as sp
import math
sp.init_printing()

עבודה על צמיגות

---

מיכל בגם ואמיר שפר

m=sp.symbols('m') #mass of the ball
eta=sp.symbols('eta') #Viscosity of liquid
raw=sp.symbols('rho') #Density of liquid
g=sp.symbols('g') #Gravitational acceleration
R=sp.symbols('R') #ball's Radius
pi=sp.symbols('pi') #pi
F=sp.symbols('F') #total Force of the ball
t=sp.var('t')
x=sp.Function('x') #place
v=sp.Function('v') #Velocity
a=sp.Function('a') #Acceleration

F_of_the_ball_1=sp.Eq(F,m*g-g*4*pi*R**3*raw/3-6*pi*eta*R*v(t))
F_of_the_ball_1

$$F = - \frac{4 R^{3} \pi \rho}{3} - 6 R \eta \pi v{\left (t \right )} + g m$$

F_of_the_ball_2=F_of_the_ball_1.subs(F,m*a(t))
F_of_the_ball_2

$$m a{\left (t \right )} = - \frac{4 R^{3} \pi \rho}{3} - 6 R \eta \pi v{\left (t \right )} + g m$$

F_of_the_ball_3=sp.Eq(a(t),g-(4*pi*R**3*raw/3)/m-(6*pi*eta*R*v(t))/m)
F_of_the_ball_3

$$a{\left (t \right )} = - \frac{4 R^{3} \pi \rho}{3 m} - \frac{6 R \eta \pi v{\left (t \right )}}{m} + g$$

F_of_the_ball_4=F_of_the_ball_3.subs([(a(t),v(t).diff())])
F_of_the_ball_4

$$\frac{d}{d t} v{\left (t \right )} = - \frac{4 R^{3} \pi \rho}{3 m} - \frac{6 R \eta \pi v{\left (t \right )}}{m} + g$$

F_of_the_ball_5=sp.dsolve(F_of_the_ball_4,v(t))
F_of_the_ball_5

$$v{\left (t \right )} = \frac{- 4 R^{2} \rho + \frac{3 g m}{R \pi} + \frac{e^{R \eta \pi \left(C_{1} - \frac{6 t}{m}\right)}}{R \pi}}{18 \eta}$$

F_of_the_ball_5_integral_x=sp.Eq(x(t),sp.integrate(F_of_the_ball_5.rhs,t))
F_of_the_ball_5_integral_x

$$x{\left (t \right )} = \begin{cases} - \frac{m e^{R \eta \pi \left(C_{1} - \frac{6 t}{m}\right)}}{108 R^{2} \eta^{2} \pi^{2}} & \text{for}\: 108 R^{2} \eta^{2} \pi^{2} \neq 0 \\t \left(- \frac{- 4 R^{3} \pi \rho + 3 g m}{18 R \eta \pi} + \frac{- 4 R^{3} \pi \rho + 3 g m + 1}{18 R \eta \pi}\right) & \text{otherwise} \end{cases} + \frac{t \left(- 4 R^{3} \pi \rho + 3 g m\right)}{18 R \eta \pi}$$

F_of_the_ball_4_2=F_of_the_ball_3.subs([(a(t),x(t).diff().diff()),(v(t),x(t).diff())])
F_of_the_ball_4_2

$$\frac{d^{2}}{d t^{2}} x{\left (t \right )} = - \frac{4 R^{3} \pi \rho}{3 m} - \frac{6 R \eta \pi \frac{d}{d t} x{\left (t \right )}}{m} + g$$

F_of_the_ball_5_2=sp.dsolve(F_of_the_ball_4_2,x(t))
F_of_the_ball_5_2

$$x{\left (t \right )} = C_{1} - \frac{2 R^{2} \rho t}{9 \eta} + \left(C_{2} - \begin{cases} - \frac{R m \rho e^{\frac{6 R \eta \pi t}{m}}}{27 \eta^{2} \pi} + \frac{g m^{2} e^{\frac{6 R \eta \pi t}{m}}}{36 R^{2} \eta^{2} \pi^{2}} & \text{for}\: 108 R^{2} \eta^{2} \pi^{2} \neq 0 \\- \frac{2 R^{2} \rho t}{9 \eta} + \frac{g m t}{6 R \eta \pi} & \text{otherwise} \end{cases}\right) e^{- \frac{6 R \eta \pi t}{m}} + \frac{g m t}{6 R \eta \pi}$$

####################################################################

#######################################################################

m_1=2.04/1000 #mass of the ball
#eta= #Viscosity of liquid
g_1=9.823 #Gravitational acceleration
R_1=0.005 #ball's Radius
pi_1=math.pi #pi

#######################################################################

#######################################################################

##### מטלות:
##### 1. לקחת את משוואת המהירות והמקום זמן ולהציב בהן את כל הנתונים הקבועים
##### 2. להציב גם את הצפיפות של אותה המדידה
##### להציב את הצמיגות שבודקים (או לולאת for)
##### 3. בעזרת V0 למצוא את הקבועה הראשון
##### 5. למצוא את הקבוע השני בעזרת x0
##### 6. ליצור רשימה ולהשוות לאמיתי

F_of_the_ball_3

$$a{\left (t \right )} = - \frac{4 R^{3} \pi \rho}{3 m} - \frac{6 R \eta \pi v{\left (t \right )}}{m} + g$$

F_of_the_ball_4

$$\frac{d}{d t} v{\left (t \right )} = - \frac{4 R^{3} \pi \rho}{3 m} - \frac{6 R \eta \pi v{\left (t \right )}}{m} + g$$

#1
F_of_the_ball_5

$$v{\left (t \right )} = \frac{- 4 R^{2} \rho + \frac{3 g m}{R \pi} + \frac{e^{R \eta \pi \left(C_{1} - \frac{6 t}{m}\right)}}{R \pi}}{18 \eta}$$

print(F_of_the_ball_5_integral_x)
F_of_the_ball_5_integral_x

Eq(x(t), Piecewise((-m*exp(R*eta*pi*(C1 - 6*t/m))/(108*R**2*eta**2*pi**2), Ne(108*R**2*eta**2*pi**2, 0)), (t*(-(-4*R**3*pi*rho + 3*g*m)/(18*R*eta*pi) + (-4*R**3*pi*rho + 3*g*m + 1)/(18*R*eta*pi)), True)) + t*(-4*R**3*pi*rho + 3*g*m)/(18*R*eta*pi))

$$x{\left (t \right )} = \begin{cases} - \frac{m e^{R \eta \pi \left(C_{1} - \frac{6 t}{m}\right)}}{108 R^{2} \eta^{2} \pi^{2}} & \text{for}\: 108 R^{2} \eta^{2} \pi^{2} \neq 0 \\t \left(- \frac{- 4 R^{3} \pi \rho + 3 g m}{18 R \eta \pi} + \frac{- 4 R^{3} \pi \rho + 3 g m + 1}{18 R \eta \pi}\right) & \text{otherwise} \end{cases} + \frac{t \left(- 4 R^{3} \pi \rho + 3 g m\right)}{18 R \eta \pi}$$

e=sp.symbols('e') #math.e
c1=sp.symbols('C1') #c1
c2=sp.symbols('C2') #c2
F_of_the_ball_x_last=sp.Eq(x(t),(-m*e**(R*eta*pi*(c1 - 6*t/m))/(108*R**2*eta**2*pi**2)+t*(-4*R**3*pi*raw + 3*g*m)/(18*R*eta*pi))+c2)
F_of_the_ball_x_last

$$x{\left (t \right )} = C_{2} + \frac{t \left(- 4 R^{3} \pi \rho + 3 g m\right)}{18 R \eta \pi} - \frac{e^{R \eta \pi \left(C_{1} - \frac{6 t}{m}\right)} m}{108 R^{2} \eta^{2} \pi^{2}}$$

F_of_the_ball_v_last=F_of_the_ball_5
F_of_the_ball_v_last

$$v{\left (t \right )} = \frac{- 4 R^{2} \rho + \frac{3 g m}{R \pi} + \frac{e^{R \eta \pi \left(C_{1} - \frac{6 t}{m}\right)}}{R \pi}}{18 \eta}$$

F_of_the_ball_v_last_sub1=F_of_the_ball_v_last.subs([(g,g_1),(m,m_1),(R,R_1),(pi,pi_1)])
F_of_the_ball_v_last_sub1

$$v{\left (t \right )} = \frac{- 0.0001 \rho + 63.6619772367581 e^{0.015707963267949 \eta \left(C_{1} - 2941.17647058824 t\right)} + 3.82715180666765}{18 \eta}$$

F_of_the_ball_x_last_sub1=F_of_the_ball_x_last.subs([(g,g_1),(m,m_1),(R,R_1),(pi,pi_1)])
F_of_the_ball_x_last_sub1

$$x{\left (t \right )} = C_{2} - \frac{0.076553783196433 e^{0.015707963267949 \eta \left(C_{1} - 2941.17647058824 t\right)}}{\eta^{2}} + \frac{3.53677651315323 t \left(- 1.5707963267949 \cdot 10^{-6} \rho + 0.06011676\right)}{\eta}$$

####### Example for water: eta=1.8

#2
raw_1=0.246/0.00025 #Density of liquid
#raw_1=998
F_of_the_ball_v_example_1=F_of_the_ball_v_last_sub1.subs([(raw,raw_1),(eta,2)])
F_of_the_ball_v_example_1

$$v{\left (t \right )} = 1.76838825657661 e^{0.0314159265358979 C_{1} - 92.3997839291116 t} + 0.103576439074101$$

#3
#C1 is related to v0:
#v0=0.1303
#(0,0.1303)
c1_example=sp.solve(F_of_the_ball_v_3.subs([(v(t),0.1303),(t,0)]),'C1')
c1_example=c1_example[0]
c1_example

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-41-6f7d2a97f343> in <module>()
      3 #v0=0.1303
      4 #(0,0.1303)
----> 5 c1_example=sp.solve(F_of_the_ball_v_3.subs([(v(t),0.1303),(t,0)]),'C1')
      6 c1_example=c1_example[0]
      7 c1_example
NameError: name 'F_of_the_ball_v_3' is not defined

F_of_the_ball_x_example_1=F_of_the_ball_x_last_sub1.subs([(raw,raw_1),(eta,2)])
F_of_the_ball_x_example_1

$$x{\left (t \right )} = C_{2} - 0.0191384457991082 e^{0.0314159265358979 C_{1} - 92.3997839291116 t} + 0.103576439074101 t$$

F_of_the_ball_x_example_2=F_of_the_ball_x_example_1.subs([(c1,c1_example)])
F_of_the_ball_x_example_2

#4
#C2 is related to x0:
#x0=0
#(0,0)
c2_example=sp.solve(F_of_the_ball_x_example_2.subs([(x(t),0),(t,0),(e,math.e)]),c2)
c2_example

F_of_the_ball_x_example_3=F_of_the_ball_x_example_2.subs([(c2,c2_example[0]),(e,math.e)])
F_of_the_ball_x_example_3

t_list=[]
x_list=np.array([])
for i in np.arange(0,61*0.05,0.05):
    t_list.append(i)
    x_list=np.append(x_list,[F_of_the_ball_x_example_3.rhs.subs([(t,i)])])
plt.plot(t_list,x_list)

x_list[-1]

t_list_example_real=[]
x_list_example_real=np.array([-0.001309,0.002946,0.007528,0.012111,0.016039,0.019967,0.023567,0.027168,0.030769,0.034042,0.037642,0.041243,0.044844,0.048117,0.05139,0.054336,0.057282,0.060555,0.063828,0.066774,0.06972,0.072994,0.076267,0.079213,0.085432,0.088705,0.091978,0.094924,0.09787,0.101143,0.104744,0.10769,0.110636,0.113909,0.117182,0.119474,0.122747,0.126348,0.130275,0.133221,0.136822,0.140423,0.143041,0.145332,0.148606,0.151552,0.155152,0.154825,0.157443,0.161044,0.163335,0.166609,0.169554,0.172173,0.175119,0.179702,0.182647,0.184611,0.188539,0.190831,0.193777])

for i in range(0,len(x_list_example_real)):
    t_list_example_real.append(i*(1/3))

plt.plot(t_list_example_real,x_list_example_real)

x_list_example_real[-1]

len(x_list)

error_present=(x_list_example_real-x_list)**2
sum(error_present)

##################################################################################

######################################################################################

#עכשיו "נגלה" מה צמיגות המים:

F_of_the_ball_3

F_of_the_ball_4

##בוא נציב ונמצא את משוואת המהירות זמן:

##נתונים קבועיים לכל המדידות
m_1=2.04/1000 #mass of the ball
g_1=9.823 #Gravitational acceleration
R_1=0.005 #ball's Radius
pi_1=math.pi #pi

#1
F_of_the_ball_v_1=F_of_the_ball_4.subs([(g,g_1),(m,m_1),(R,R_1),(pi,pi_1)])
F_of_the_ball_v_1

#!#
## נתון של המדידה הנוכחית:
raw_1=0.246/0.00025 # ץ(מסת החומר חלקי הנפח שלה)   צפיפות החומר    

F_of_the_ball_v_1_2=F_of_the_ball_v_1.subs([(raw,raw_1)])
F_of_the_ball_v_1_2  #חסר רק צמיגות      משוואה דיפרנציאלית עם פונקציית מהירות זמן מוצבת

print(F_of_the_ball_5_integral_x)
F_of_the_ball_5_integral_x

e=sp.symbols('e') #math.e
c1=sp.symbols('C1') #c1
c2=sp.symbols('C2') #c2
F_of_the_ball_x_last=sp.Eq(x(t),(-m*e**(R*eta*pi*(c1 - 6*t/m))/(108*R**2*eta**2*pi**2)+t*(-4*R**3*pi*raw + 3*g*m)/(18*R*eta*pi))+c2)
F_of_the_ball_x_last  ### פונקציית מקום זמן סופי עם פרמטרים

F_of_the_ball_x_last_sub1=F_of_the_ball_x_last.subs([(g,g_1),(m,m_1),(R,R_1),(pi,pi_1)])
F_of_the_ball_x_last_sub1 ## עדיין אין צפיפות וצמיגות       פונקציית מקום זמן עם הנתונים הקבועיים

F_of_the_ball_x_last_sub2=F_of_the_ball_x_last_sub1.subs([(raw,raw_1)])
F_of_the_ball_x_last_sub2 #חסר רק צמיגות ושני הקבועיים      פונקציית מקום זמן מוצבת

####3
F_of_the_ball_v_2=F_of_the_ball_v_1_2.subs([(eta,2)])
F_of_the_ball_v_2
##4
F_of_the_ball_v_3=sp.dsolve(F_of_the_ball_v_2)
F_of_the_ball_v_3

#!#
#בוא נמצא את הקבוע בעזרת המהירות ההתחלתית של המדידה הנוכחית
#C1 is related to v0:
#(0,0.1303)
v0=0.1303

c1_example=sp.solve(F_of_the_ball_v_3.subs([(v(t),v0),(t,0)]),'C1')
c1_example=c1_example[0]
c1_example
F_of_the_ball_x_last_sub3=F_of_the_ball_x_last_sub2.subs([(eta,2)])
F_of_the_ball_x_last_sub3 #צמיגות בו, רק שני הקבועים חסרים
F_of_the_ball_x_last_sub4=F_of_the_ball_x_last_sub3.subs([(c1,c1_example)])
F_of_the_ball_x_last_sub4 #
F_of_the_ball_x_last_sub4

#4
#C2 is related to x0:
#x0=0
#(0,0)
c2_example=sp.solve(F_of_the_ball_x_last_sub4.subs([(x(t),0),(t,0),(e,math.e)]),c2)
c2_example
F_of_the_ball_x_last_last=F_of_the_ball_x_last_sub4.subs([(c2,c2_example[0]),(e,math.e)])
F_of_the_ball_x_last_last

t_list=[]
x_list=np.array([])
for i in np.arange(0,61*0.05,0.05):
    t_list.append(i)
    x_list=np.append(x_list,[F_of_the_ball_x_last_last.rhs.subs([(t,i)])])
t_list_example_real=[]
x_list_example_real=np.array([-0.001309,0.002946,0.007528,0.012111,0.016039,0.019967,0.023567,0.027168,0.030769,0.034042,0.037642,0.041243,0.044844,0.048117,0.05139,0.054336,0.057282,0.060555,0.063828,0.066774,0.06972,0.072994,0.076267,0.079213,0.085432,0.088705,0.091978,0.094924,0.09787,0.101143,0.104744,0.10769,0.110636,0.113909,0.117182,0.119474,0.122747,0.126348,0.130275,0.133221,0.136822,0.140423,0.143041,0.145332,0.148606,0.151552,0.155152,0.154825,0.157443,0.161044,0.163335,0.166609,0.169554,0.172173,0.175119,0.179702,0.182647,0.184611,0.188539,0.190831,0.193777])

for i in range(0,len(x_list_example_real)):
    t_list_example_real.append(i*(1/3))
error_present=(x_list_example_real-x_list)**2
sum(error_present)

best_eta=0
min_error=1000
for eta_1 in np.arange(1.8,2.2,0.1):
    ####3
    F_of_the_ball_v_2=F_of_the_ball_v_1_2.subs([(eta,eta_1)])
    F_of_the_ball_v_2
    ##4
    F_of_the_ball_v_3=sp.dsolve(F_of_the_ball_v_2)
    F_of_the_ball_v_3
    #!#
    #בוא נמצא את הקבוע בעזרת המהירות ההתחלתית של המדידה הנוכחית
    #C1 is related to v0:
    #(0,0.1303)
    v0=0.1303
    c1_example=sp.solve(F_of_the_ball_v_3.subs([(v(t),v0),(t,0)]),'C1')
    c1_example=c1_example[0]
    c1_example
    F_of_the_ball_x_last_sub3=F_of_the_ball_x_last_sub2.subs([(eta,eta_1)])
    F_of_the_ball_x_last_sub3 #צמיגות בו, רק שני הקבועים חסרים
    F_of_the_ball_x_last_sub4=F_of_the_ball_x_last_sub3.subs([(c1,c1_example)])
    F_of_the_ball_x_last_sub4 #קבוע המהירות הוחלף
    #4
    #C2 is related to x0:
    #x0=0
    #(0,0)
    c2_example=sp.solve(F_of_the_ball_x_last_sub3.subs([(x(t),0),(t,0),(e,math.e)]),c2)
    c2_example
    F_of_the_ball_x_last_last=F_of_the_ball_x_last_sub4.subs([(c2,c2_example[0]),(e,math.e)])
    F_of_the_ball_x_last_last
    t_list=[]
    x_list=np.array([])
    for i in np.arange(0,61*0.05,0.05):
        t_list.append(i)
        x_list=np.append(x_list,[F_of_the_ball_x_last_last.rhs.subs([(t,i)])])
    t_list_example_real=[]
    x_list_example_real=np.array([-0.001309,0.002946,0.007528,0.012111,0.016039,0.019967,0.023567,0.027168,0.030769,0.034042,0.037642,0.041243,0.044844,0.048117,0.05139,0.054336,0.057282,0.060555,0.063828,0.066774,0.06972,0.072994,0.076267,0.079213,0.085432,0.088705,0.091978,0.094924,0.09787,0.101143,0.104744,0.10769,0.110636,0.113909,0.117182,0.119474,0.122747,0.126348,0.130275,0.133221,0.136822,0.140423,0.143041,0.145332,0.148606,0.151552,0.155152,0.154825,0.157443,0.161044,0.163335,0.166609,0.169554,0.172173,0.175119,0.179702,0.182647,0.184611,0.188539,0.190831,0.193777])

    for i in range(0,len(x_list_example_real)):
        t_list_example_real.append(i*(1/3))
    error_present=(x_list_example_real-x_list)**2
    final_error=sum(error_present)
    if(final_error<min_error):
        min_error=final_error
        best_eta=eta_1
best_eta

####################################################################################

F_of_the_ball_v_last

F_of_the_ball_v_lastF_subs1

F_of_the_ball_4

F_of_the_ball_4
F_of_the_ball_4_2=sp.dsolve(F_of_the_ball_4,v(t))

F_of_the_ball_4_2.subs([(g,g_1),(m,m_1),(R,R_1),(pi,pi_1),(raw,raw_1),(eta,2)])

(0,1.303)
c1=sp.solve(F_of_the_ball_v_3.subs([(v(t),1.303),(t,0)]),'C1')
c1

F_of_the_ball_v_4=F_of_the_ball_v_3.subs([('C1',c1[0])])
F_of_the_ball_v_4

#5
sp.integrate(F_of_the_ball_v_4,t)

F_of_the_ball_x_from_v_1=sp.Eq(x(t),sp.integrate(F_of_the_ball_v_4.rhs,t))
F_of_the_ball_x_from_v_1

t_list=[]
x_list=[]
for i in np.arange(0,5,0.05):
    t_list.append(i)
    x_list.append(F_of_the_ball_x_from_v_1.rhs.subs([(t,i)]))
plt.plot(t_list,x_list)

x_list[t_list.index(0.5)]

#The eq:
F_of_the_ball_v_1