目錄

1.Python解微分方程數值解

2.驗證火箭發射模型

---

1.Python解微分方程數值解

Python解微分方程要用到幾個庫：numpy, matplotlib.pyplot, scipy.integrate，沒有的話就

pip install 相應的庫就行，本次用的python為3.6.8

我們先來看一下簡單的微分方程

?對于Python求解微分方程只需要跳相應的庫即可

from typing import Listimport numpy as np import matplotlib.pyplot as plt from scipy.integrate import odeint# 一階微分方程的例子 def diff_equation(y_list: List[float], x: float):y1, y2, y3 = y_listreturn np.array([y2 * y3, -y1 * y3, -0.51 * y1 * y2]) # 微分方程格式，左邊一定是dy/dx,返回右邊x = np.linspace(0, 12, 100) # 初始點是0,終點是12，其中有100個點 result = odeint(diff_equation, [0, 1, 1], x) # 中間那個是y初值 plt.plot(x, result[:, 0], label='y1') # result整個矩陣的第一列 plt.plot(x, result[:, 1], label='y2') plt.plot(x, result[:, 2], label='y3')plt.legend() plt.grid() # 網格 plt.show() # 這是y=f(x)的圖像

結果展示

2.驗證火箭發射模型

掌握了解方程之后我們就可以來驗證火箭發射模型了

高階微分方程，我們可以化為微分方程組來解，在這里書上義給出微分方程組，但只給出了0~t1階段即火箭燃料還有的階段，t1~t2為燃料耗盡的時間階段，其中r=18，t1~t2的微分方程組為 -g+(-k*x2*x2/m),這里的m是火箭凈重量。

為了驗證以下表格

代碼如下?

from typing import Listimport matplotlib.pyplot as plt import numpy as np from scipy.integrate import odeintdef get_acceleration(x: np.ndarray, result: np.ndarray, time: float, f: float, m: float, g: float,k: float) -> np.ndarray:"""獲取火箭發射的加速度:param x: 距離數組:param result::param time: 時間:param f: 外力:param m: 質量:param g :重力加速度:param k : 阻力系數:return: 加速度數組"""temp = []# 重力加速度result = result[:, 1]for i in range(np.size(result)):if x[i] < time:temp.append(-g + (f - k * pow(result[i], 2)) / (m - x[i] * 18))else:temp.append(-g + (- k * pow(result[i], 2)) / 520)return np.array(temp)def diff_equation1(x_list: List[float], x: float):x1, x2 = x_listtemp = ((27000 - 0.3 * x2 * x2) / (1600 - x * 18)) if x <= 60 else -0.3 * x2 * x2 / 520return np.array([x2, -9.8 + temp]) # 微分方程格式，左邊一定是dy/dx,返回右邊def diff_equation2(x_list: List[float], x: float):x1, x2 = x_listtemp = ((27000 * 2 - 0.3 * x2 * x2) / (1600 - x * 18)) if x <= 60 else -0.3 * x2 * x2 / 520return np.array([x2, -9.8 + temp]) # 微分方程格式，左邊一定是dy/dx,返回右邊def diff_equation3(x_list: List[float], x: float):x1, x2 = x_listtemp = ((27000 - 0.3 * x2 * x2) / (2680 - x * 18)) if x <= 120 else -0.3 * x2 * x2 / 520return np.array([x2, -9.8 + temp]) # 微分方程格式，左邊一定是dy/dx,返回右邊def diff_equation(x_list: List[float], x: float):x1, x2 = x_listtemp = ((27000 * 2 - 0.3 * x2 * x2) / (2680 - x * 18)) if x <= 120 else -0.3 * x2 * x2 / 520return np.array([x2, -9.8 + temp]) # 微分方程格式，左邊一定是dy/dx,返回右邊x = np.linspace(0, 140, 10000) # 初始點是0 x1 = np.linspace(0, 80, 10000) result = odeint(diff_equation, [0., 0.], x) # 中間那個是y0初值，即x=0時y=1 result1 = odeint(diff_equation1, [0., 0.], x1) result2 = odeint(diff_equation2, [0., 0.], x1) result3 = odeint(diff_equation3, [0., 0.], x)plt.plot(x, get_acceleration(x, result, 120, 27000 * 2, 2680, 9.8, 0.3), label='a3(t)') plt.plot(x1, get_acceleration(x1, result2, 60, 27000 * 2, 1600, 9.8, 0.3), label='a2(t)') plt.plot(x, get_acceleration(x, result3, 120, 27000, 2680, 9.8, 0.3), label='a1(t)') plt.plot(x1, get_acceleration(x1, result1, 60, 27000, 1600, 9.8, 0.3), label='a0(t)')plt.legend() plt.grid() # 網格 plt.show() # 這是微分方程的圖像

結果展示

x的關系圖

?v的關系圖

a的關系圖

?閑著沒事干寫的，寫的不是很好，有疑問可以發我郵箱liuzhi_wdq@foxmail.com

媽媽，他們拋棄了我 像歌唱一樣拋棄了我

媽媽，我是多么愛你 當你沉默的時候我愛你

只是那些猛烈的情緒

在睡不著的時候折磨著我

我那早已死去的父親

在沒有星星的夜晚看著你

媽媽，我會在夏天開放嗎

像你曾經的容顏那樣

媽媽，這種失落會持久嗎

這個世界會好嗎

忘記一些隱秘的委屈

在回頭觀望的時候迷失了自己

我的正在老去的身體

從某一天開始就在漸漸死去

媽媽我愛你

媽媽，我居然愛上了她

像歌唱一樣就愛上了她

媽媽，當你又回首一切

這個世界會好嗎

媽媽，我是多么恨你

當我歌唱的時候我恨你

總結

以上是生活随笔為你收集整理的Python解微分方程(验证数学建模第五版火箭发射模型)的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。