2018-10-06-152433.ipynb

Kernel: Python 3 (Ubuntu Linux)

import sympy as sp
import matplotlib.pyplot as plt
import numpy as np
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_of_the_ball_2=F_of_the_ball_1.subs(F,m*a(t))
F_of_the_ball_2

F_of_the_ball_3=sp.Eq(a(t),g-(4*g*pi*R**3*raw/3)/m-(6*pi*eta*R*v(t))/m)
F_of_the_ball_4=F_of_the_ball_3.subs([(a(t),v(t).diff())])
F_of_the_ball_4 #משוואה דיפרנציאלית

#נתונים קבועים לכל המדידות
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
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  #חסר רק צמיגות      משוואה דיפרנציאלית עם פונקציית מהירות זמן מוצבת

#x_list
t_list_example_real=[]
x_list_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])

min_error=1000
min_eta=-1
x_list_min=[]
for eta_1 in np.arange(0.8,1.8,0.05):
    v0=0.1303
    x0=0
    F_of_the_ball_v_1_3=F_of_the_ball_v_1_2.subs([(eta,eta_1)])
    F_of_the_ball_v_solve_1=sp.dsolve(F_of_the_ball_v_1_3,v(t))
    c1_solved=sp.solve(F_of_the_ball_v_solve_1.subs([(v(t),v0),(t,0)]),'C1')
    c1_solved=c1_solved[0]
    F_of_the_ball_v_solve_2=F_of_the_ball_v_solve_1.subs('C1',c1_solved)
    c2=sp.symbols('C2')
    F_of_the_ball_x_1=sp.integrate(F_of_the_ball_v_solve_2.rhs,t)+c2
    F_of_the_ball_x=F_of_the_ball_x_1.subs(c2,x0)
    x_list=[]
    n=len(x_list_real)
    t_list=np.arange(0,n*1/30,1/30)
    for i in t_list:
        x_list=np.append(x_list,[F_of_the_ball_x.subs([(t,i)])])
    error=(1/(2*n)) * sum((x_list-x_list_real)**2)
    if error<min_error:
        min_error=error
        min_eta=eta_1
        x_list_min=x_list
print(min_eta)
plt.plot(t_list,x_list_real ,)
plt.plot(t_list,x_list_min ,)
plt.xlabel("t")
plt.ylabel('x')

#in 1 code:
def find_Viscosity(x_list_real,density_of_liquid):
    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_2=F_of_the_ball_1.subs(F,m*a(t))
    F_of_the_ball_3=sp.Eq(a(t),g-(4*g*pi*R**3*raw/3)/m-(6*pi*eta*R*v(t))/m)
    F_of_the_ball_4=F_of_the_ball_3.subs([(a(t),v(t).diff())])
    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
    F_of_the_ball_v_1=F_of_the_ball_4.subs([(g,g_1),(m,m_1),(R,R_1),(pi,pi_1)])
    raw_1=density_of_liquid
    F_of_the_ball_v_1_2=F_of_the_ball_v_1.subs([(raw,raw_1)])
    min_error=1000
    min_eta=-1
    x_list_min=[]
    for eta_1 in np.arange(0,10,0.05):
        v0=0.1303
        x0=0
        F_of_the_ball_v_1_3=F_of_the_ball_v_1_2.subs([(eta,eta_1)])
        F_of_the_ball_v_solve_1=sp.dsolve(F_of_the_ball_v_1_3,v(t))
        c1_solved=sp.solve(F_of_the_ball_v_solve_1.subs([(v(t),v0),(t,0)]),'C1')
        c1_solved=c1_solved[0]
        F_of_the_ball_v_solve_2=F_of_the_ball_v_solve_1.subs('C1',c1_solved)
        c2=sp.symbols('C2')
        F_of_the_ball_x_1=sp.integrate(F_of_the_ball_v_solve_2.rhs,t)+c2
        F_of_the_ball_x=F_of_the_ball_x_1.subs(c2,x0)
        x_list=[]
        n=len(x_list_real)
        t_list=np.arange(0,n*0.05,0.05)
        for i in t_list:
            x_list=np.append(x_list,[F_of_the_ball_x.subs([(t,i)])])
        error=(1/(2*n)) * sum((x_list-x_list_real)**2)
        if error<min_error:
            min_error=error
            min_eta=eta_1
            x_list_min=x_list
    return(min_eta)

find_Viscosity(x_list_real,1000)

$$2.35$$