空気抵抗がある場合の球の軌跡 微分方程式解析

球を投げます。

抵抗係数k 0.47 (kg/m)

質量使っていない

初速度 41 m/s

コード

    

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

g = 9.8
m = 0.145
k = 0.47

def ball(v,t):
    dvdt=np.zeros_like(v)
    dvdt[0] = -k*(v[0]**2+v[1]**2)*v[0]/(np.sqrt(v[0]**2+v[1]**2))
    dvdt[1] = -k*(v[0]**2+v[1]**2)*(v[1]/(np.sqrt(v[0]**2+v[1]**2)))-g
    dvdt[2] = v[0]
    dvdt[3] = v[1]
    return dvdt

t = np.arange(0, 10, 0.001)
kyori = [0]
kakudo = [0]

for theta in np.linspace(0,np.pi/2,101):
    

    v0 = [41*np.cos(theta), 41*np.sin(theta), 0, 0]

    y = odeint(ball, v0, t, args=())

    plt.plot(y[:, 2],y[:,3])
    

    for i in np.linspace(1, 10000, 10000):
        
        if float(y[int(i),3]) <= 0:
            y1 = y[int(i), 2]
            kyori.append(y1)
            kakudo.append(theta)
            break
    
    
plt.ylim(-1,7)
plt.xlabel('x')
plt.ylabel('y')
plt.show()

print(kyori)

plt.xlabel('theta')
plt.ylabel('distance')
plt.plot(kakudo, kyori)

    
maxx=max(kyori)
maxi=np.argmax(kyori)

print(maxx)
print(maxi)
print(180*kakudo[maxi]/np.pi)

    

––

少し解説

[Vx, Vy, x, y] って感じです。

数値解析  離散的なので  yがぴったり0になる点が 示されているか分からないので

 if y[int(i),3] <= 0:

yが0以下になった時の xの値の最大値を求めたい。

for文で回しています。

append で後ろに追加していきます。

--

10000個  4つ

-

max  最大値を取得

argmax 変数の最大値を取得

ーーー

軌跡を表したもの 100個

10000個出したものの内

30番目のものが 最大

投げる角度と飛距離の関係

0.455530934771

弧度法で表して

26.1°で  最長飛距離 5.92766280504 m

ーーー

デバッグのために

至る所で print(名)  していました。

……