A story about calculating the speed of a small ball falling while receiving air resistance with Python and Sympy

Introduction

I remember playing with Mathematica while I was in a university lab. Over the decades since then, differential equations have appeared in front of me many times in my research life, and I wish I had computer algebra software, but when it came to numerical integration. I'm an old engineer who has just reached the end of the calendar with the corruption. Nice to meet you tonight. Recently, however, a mathematician named Mr. Tanaka, a metamorphosis companion, introduced me to Sympy, and while I was using it, my memories of my youth revived, and I went to Tobita Shinchi. I decided to try to solve the problem.

environment

Windows10 pro Anaconda + Python 3.7.6 (64bit) + Sympy 1.5.1 Visual Studio Code 1.44.2

Basic formula

In high school, an old engineer first learned the power of differential equations when he encountered the problem of falling during air resistance.

A small ball with a mass of $ m $ (kg) is falling vertically downward. The globules start falling at the initial speed of $ 0 $ at time $ 0 $ (s), and during the fall, the globules have an air resistance $ kv $ (N) proportional to the speed $ v $ (m / s). Find the speed of the globules at time $ t $ (s), assuming that they act in the opposite direction of. </ b>

A long time ago, there was a great physics teacher named Mr. Sakuma in Sundai, who wrote the differential equations in a quick letter. I was in the first year of high school and was diving in the morning, so I flinch for a moment. However, let's first formulate the equation of motion while taking pride in moss-sucking, as the days of copying it into a notebook somehow continue and it is sunny and active in science class 1. Let the speed of the small ball at time $ t $ (s) be $ v (t) $ (m / s)

Equation of motion </ b>

m \\frac{dv\\left(t \\right)}{d t} {} = mg-kv{\\left(t \\right)}

Solve this $v{\\left(t \\right)} = \\frac{mg}{k} + \\frac{e^{\\frac{C_{1} k}{m}} }{k}e^{- \\frac{k t}{m}}$ Constant $A=\\frac{e^{\\frac{C_{1} k}{m}} }{k}$ To summarize $v{\\left(t \\right)} = \\frac{mg}{k} + A e^{- \\frac{k t}{m}}$ It becomes a general solution. When v {\ left (t \ right)} = 0 $ is inserted under the initial condition $ t = 0 (s) to obtain $ A , $A=-\frac{mg}{k}$ And $v{\left(t \right)} = \frac{mg}{k} (1- e^{- \frac{k t}{m}})$$ Can be obtained. Speed $ v {\ left (t \ right)} $ t → \ infinity $ limit, that is, terminal velocity $\lim_{t \to \infty} v(t)=\\frac{mg}{k}$ It will be.

Solve this with Sympy and draw a graph

python


#%%
import sympy as sym
from IPython.display import Math, display
from sympy.solvers import ode
from sympy import oo
import matplotlib.pyplot as plt

#Define variables
t, v, k, m, g = sym.symbols("t v k m g")
v = sym.Function('v')

#Equation of motion (ordinary differential equation)
eq = sym.Eq(m * v(t).diff(t,1) - m * g + k * v(t), 0)

#General solution
ans = sym.expand(sym.dsolve(eq))
display(Math(f"General solution: {sym.latex(ans)}"))

#Special solution(Character expression)
ans2 = sym.expand(sym.dsolve(eq, ics={v(0):0}))
display(Math(f"Special solution(Character expression): {sym.latex(ans2)}"))

#Special solution (numerical value)
m1 = 0.4
k1 = 0.2
g1 = 9.8
eq2 = eq.subs([(m, m1), (k, k1), (g, g1)])
ans3 = sym.dsolve(eq2, ics={v(0):0})
p1 = sym.plot(ans3.rhs, (t, 0, 10), show=False) #rhs is the right side
display(Math(f"Special solution(m={m1}, k={k1}, g={g1}): {sym.latex(ans3)}"))

#Terminal velocity
voo = sym.limit(ans3.rhs, t, oo)
display(Math(r"Terminal velocity:\lim_{t \to \infty} v(t)" \
    f"={voo:.2f}"))
p2 = sym.plot(voo, (t, 0, 10), show=False)
p1.extend(p2)
p1[1].line_color = "0.8"
p1.show()
fig = p1._backend.fig
fig.legend([f"f(t)={ans3.rhs}", f"f(t)={voo:.1f}"], loc='center')
fig

When you do this, image.png image.png

The good thing about Sympy is that you can ask for a special solution. Furthermore, if you pass the formula of the special solution to sympy.plot, you can draw a graph. You can draw an asymptote for the limit.

Recommended Posts