概述

之前我一直用Matlab计算微分方程，这次尝试下用python进行计算，目标是一个二阶的微分方程组，方程组的物理意义是气体动力学建立的乒乓球某一时间段的轨迹。

题目

题目描述如下：
目标方程为一组参数方程$x(t),y(t)$，有如下定义:

$V_x=\frac{dx(t)}{dt} \\ V_y=\frac{dy(t)}{dt} \\ \theta = arctan\frac{V_x}{V_y} \\ V=(V_x^2+V_y^2)^{\frac{1}{2}}$

已知的条件包括微分方程组和方程初始值。微分方程组形式如下：

$\frac{d^2}{dt^2}x(t) = \frac{1}{m}(F_3 \cdot V \cdot sin\theta - F_2 \cdot V^2 \cdot cos \theta) \\ \frac{d^2}{dt^2}y(t) = - \frac{1}{m}(G - F_1 + F_3\cdot V \cdot cos\theta + F_2 \cdot V \cdot sin \theta )$

其中参数值为：$m = 0.0025$，$G = 0.0245$，$F_1 = 3.6406E-4$，$F_2=1.4664E-4$，$F_3=0.008$。
已知的初始值为：$x(0) = 0$，$y(0) = 0.3$，$V_x(0)=42$，$V_y(0)=0.3$。
要求画出以$x(t)$为横坐标，以$y(t)$为纵坐标，且满足$y(t)>0,t>0$的轨迹图像。

思路

我利用了scipy.integrate库提供的数值积分和常微分方程组求解方法odeint。实现参考了官方文档，但是官方文档中仅仅给出了一阶微分方程组的示例。
我引入了两个原始方程的一阶微分方程：$x_1(t)和y_1(t)$，这样需要求解的方程就包含了四个，即两个原始方程以及它们的一阶微分方程分别为$x_0(x),x_1(t),y_0(t),y_1(t)$，将上述微分方程进行化简，得到如下方程组：

$\frac{d}{dt}x\_0(t) = x\_1(t) \\ \frac{d}{dt}y\_0(t) = y\_1(t) \\ \frac{d}{dt}x\_1(t) = \frac{1}{m}(F_3 \cdot y\_1(t) - F_2 \cdot x\_1(t) \cdot V ) \\ \frac{d}{dt}y\_1(t) = - \frac{1}{m}(G - F_1 + F_3\cdot x\_1(t) + F_2 \cdot y\_1(t))$

再得到以上化简结果后利用odeint方法求解并做图即可。

代码

关于odeint的参数，查看官方文档即可，完成程序如下：

# -*- coding: utf-8 -*-
# By YanWang
# 2018 4 16 

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

# coefficient
coe_m = 0.0025
coe_G = 0.0245
coe_F1 = 3.6406E-4
coe_F2 = 1.4664E-4
coe_F3 = 0.008

# equation
def equation(w,t,m,G,F1,F2,F3):
    X_0,X_1,Y_0,Y_1 = w

    V = np.power(np.power(X_1,2.0) + np.power(Y_1,2.0),0.5)
    diff_X0 = X_1
    diff_X1 = (F3*Y_1 - F2*X_1*V)/m
    diff_Y0 = Y_1
    diff_Y1 = -(G - F1 + F3*X_1 + F2*Y_1)/m

    return np.array([diff_X0,diff_X1,diff_Y0,diff_Y1])

# the time which we want to get the data 
time = np.arange(0,0.7,0.0001)

# use "odeint" to count the data
track = odeint(equation,(0,42,0.3,0.3),time,args=(coe_m,coe_G,coe_F1,coe_F2,coe_F3))

# print the data
print track[:,0]
print track[:,2]

# get the photo
plt.plot(track[:,0], track[:,2])
plt.show()

需要注意的是，odeint返回的是四个方程的时间区间内的值，根据我代码中的求解方程的顺序，找到$x(t),y(y)$对应的序号为第一组解和第三组解，所以输出为track[:,0]和track[:,2]。

结果与分析

一开始我赋值时间为:

1	time = np.array(0,30,0.01)

然后得到结果如下：
$photo1$
很明显，乒乓球不会运动的这么放飞自我，然后我输出了部分数据，发现数据的前一部分基本复合预定轨迹，所以问题的原因是时间区间选择不合适，时间过过久导致模型因超出定义域而失去物理意义。根据输出数据值调整时间区间，最后得到正常轨迹图如下：
$photo2$
至此问题解决。

小结

我觉得相比Matlab，pyhton作为科学计算的工具可能更和我的胃口，以上。

参考

[1] SciPy官方文档
[2] Treysure的博客