import sympy as smp from matplotlib import pyplot as plt import numpy as nmp from scipy.integrate import solve_ivp p = 1.225 # air density m = 0.003 # mass of the plane a = (9.3/180)*nmp.pi # angle of attack ar = 0.86 # aspect ratio s = 0.017 # wing area g = 9.807 # acceleration of free fall cla = (nmp.pi*ar)/1+nmp.sqrt(1+(ar/2)**2) # lift slope derivative cl = cla*a # lift coefficient e = 0.9 # oswald efficiency factor ε = 1/(nmp.pi*e*ar) # induced darg factor cd = 0.02 + ε*(cl**2) def airplane(t,f): v, γ, h, r = f return [-cd(p*v**2)*s/(2*m)-g*nmp.sin(γ), ((cl*(p*v**2/2))-g*nmp.cos(γ))/v, v*nmp.sin(γ), v*nmp.cos(γ)] solution = solve_ivp(airplane, (0,1000),[3.7,((-11/180)*nmp.pi),2.0,0.0]) print(solution)

时间: 2023-12-29 16:05:27 浏览: 29
这段代码是一个飞机在空气中运动的模拟。它使用了Sympy、Matplotlib、NumPy和SciPy库。 首先,定义了一些变量和常量,如空气密度(p)、飞机质量(m)、攻角(a)、展弦比(ar)、翼面积(s)、自由落体加速度(g)等。 然后,定义了一个名为"airplane"的函数,它描述了飞机在空中的运动方程。这个方程包括了飞机的速度(v)、爬升角(γ)、高度(h)和航向(r)随时间变化的关系。 最后,使用solve_ivp函数对"airplane"函数进行求解。该函数接受初始条件(速度、爬升角、高度和航向)和时间范围(0到1000秒),并返回飞机在这段时间内的运动轨迹。 最后一行代码打印出了求解的结果。 需要注意的是,代码中涉及到的数学模型和参数是根据特定的假设和公式进行计算的,具体的物理背景和准确性需要根据实际情况进行评估和调整。
相关问题

import matplotlib.pyplot as plt import matplotlib as mpl import numpy as np from scipy import integrate import sympy import mpmath

这段代码导入了matplotlib、numpy、scipy、sympy和mpmath这五个库。其中,matplotlib是用于画图的库,numpy是用于数值计算的库,scipy是用于科学计算的库,sympy是用于符号计算的库,mpmath是用于高精度计算的库。这些库都是Python中常用的科学计算库。

import sympy from scipy.integrate import quad y=x*sympy.atan(x)*sympy.ln(sympy.exp(1+x**2)) result,error=quad(y,0,1) print(result)出错怎么办

你的代码存在一些语法错误。可以尝试以下代码: ```python import sympy from scipy.integrate import quad x = sympy.symbols('x') y = x*sympy.atan(x)*sympy.ln(sympy.exp(1+x**2)) result, error = quad(sympy.lambdify(x, y), 0, 1) print(result) ``` 在这里,我们使用 `sympy.symbols` 函数来定义符号变量 `x`,并使用 `sympy.lambdify` 函数将 `y` 转换为可用于积分计算的函数。 另外,如果你想要输出误差值,可以将 `print` 语句改为以下代码: ```python print(f"Result: {result}, Error: {error}") ```

相关推荐

import numpy as np import matplotlib.pyplot as plt import sympy from scipy.interpolate import interp1d gamma = 1.2 R = 8.314 T0 = 500 Q = 50 * R * T0 a0 = np.sqrt(gamma * R * T0) M0 = 6.216 P_P0 = sympy.symbols('P_P0') num = 81 x0 = np.linspace(0,1,num) t_t0 = np.linspace(0,15,num) x = x0[1:] T_T0 = t_t0[1:] h0 = [] h1 = []#创建拉姆达为1的空数组 r = [] t = [] c = [] s = [] i = 0 for V_V0 in x: n1 = sympy.solve(1 / (gamma-1) * (P_P0 * V_V0 - 1) - 0.5 * (P_P0 + 1) * (1 - V_V0)- gamma * 0 * Q / a0 ** 2,P_P0)#lamuda=0的Hugoniot曲线方程 n2 = sympy.solve(1 / (gamma-1) * (P_P0 * V_V0 - 1) - 0.5 * (P_P0 + 1) * (1 - V_V0)- gamma * 1 * Q / a0 ** 2,P_P0)#lamuda=1的Hugoniot曲线方程 n3 = sympy.solve(-1 * P_P0 + 1 - gamma * M0 ** 2 * (V_V0 - 1),P_P0)#Reyleigh曲线方程 n4 = 12.014556 / V_V0#等温线 n5 = sympy.solve((P_P0 - 1 / (gamma+1) )* (V_V0-gamma / (gamma + 1)) - gamma / ((gamma + 1) ** 2),P_P0)#声速线 n6 = 10.6677 / np.power(V_V0,1.2)#等熵线 h0.append(n1) h1.append(n2) r.append(n3) t.append(n4) c.append(n5) s.append(n6) i = i+1 h0 = np.array(h0) h1 = np.array(h1) r = np.array(r) t = np.array(t) c = np.array(c) s = np.array(s) plt.plot(x,r,label='Rayleigh') plt.plot(x,t,color='purple',label='isothermal') plt.plot(x,s,color='skyblue',label='isentropic') a = np.where(h0 < 0) b = np.where(c < 0) h0 = np.delete(h0,np.where(h0 < 0)[0],axis = 0)#去除解小于0的值 h1 = np.delete(h1,np.where(h1 < 0)[0],axis = 0)#去除解小于0的值 c = np.delete(c,np.where(c < 0)[0],axis = 0)#去除解小于0的值 x0 = np.delete(x,a,axis = 0)#对应去除x轴上错误值的坐标 x1 = np.delete(x,b,axis = 0) plt.plot(x0,h0,label='Hugoniot(lambda=0)') plt.plot(x0,h1,label='Hugoniot(lambda=1)') plt.plot(x1,c,color='yellow',label='soniclocus') plt.ylim((0,50)) plt.legend() # 显示图例 plt.xlabel('V/V0') plt.ylabel('P/P0') f1 = interp1d(x1, c.T, kind='cubic') f2 = interp1d(x,r.T,kind='cubic') f3 = interp1d(x, t.T, kind='cubic') epsilon = 0.0001 x0 = 0.56 y0 = f1(x0) - f2(x0) while abs(y0) > epsilon: df = (f1(x0 + epsilon) - f2(x0 + epsilon) - y0) / epsilon x0 -= y0 / df y0 = f1(x0) - f2(x0) plt.scatter(x0, y0, 50, color ='red') plt.show()

#外点法(能运行出来) import math import sympy import numpy as np from matplotlib import pyplot as plt from mpl_toolkits.mplot3d import Axes3D plt.ion() fig = plt.figure() ax = Axes3D(fig) def draw(x,index,M): # F = f + MM * alpha # FF = sympy.lambdify((x1, x2), F, 'numpy') Z = FF(*(X, Y,M)) ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='rainbow',alpha=0.5) ax.scatter(x[0], x[1], FF(*(x[0],x[1],M)), c='r',s=80) ax.text(x[0], x[1], FF(*(x[0],x[1],M)), 'here:(%0.3f,%0.3f)' % (x[0], x[1])) ax.set_zlabel('F') # 坐标轴 ax.set_ylabel('X2') ax.set_xlabel('X1') plt.pause(0.1) # plt.show() # plt.savefig('./image/%03d' % index) plt.cla() C = 10 # 放大系数 M = 1 # 惩罚因子 epsilon = 1e-5 # 终止限 x1, x2 = sympy.symbols('x1:3') MM=sympy.symbols('MM') f = -x1 + x2 h = x1 + x2 - 1 # g=sympy.log(x2) if sympy.log(x2)<0 else 0 g = sympy.Piecewise((x2-1, x2 < 1), (0, x2 >= 1)) # u=lambda x: alpha = h ** 2 + g ** 2 F = f + MM * alpha # 梯度下降来最小化F def GD(x,M,n): # F = f + M * alpha # delta_x = 1e-11 # 数值求导 # t = 0.0001 # 步长 e = 0.001 # 极限 # my_print(e) np.array(x) for i in range(15): t = sympy.symbols('t') grad = np.asarray( [sympy.diff(F, x1).subs([(x1, x[0]), (x2, x[1]),(MM,M)]), sympy.diff(F, x2).subs([(x1, x[0]), (x2, x[1]),(MM,M)])]) # print('g',grad) # print((x-t*grad)) # print(F.subs([(x1,(x-t*grad)[0]),(x2,(x-t*grad)[1])])) t = sympy.solve(sympy.diff(F.subs([(x1, (x - t * grad)[0]), (x2, (x - t * grad)[1]),(MM,M)]), t), t) print('t',t) x = x - t * grad print('x', x) # print('mmm',M) draw(x,n*10+i,M) # my_print(np.linalg.norm(grad)) # print(type(grad)) if (abs(grad[0]) < e and abs(grad[1]) < e): # print(np.linalg.norm(grad)) print('g', grad) break return list(x) pass x = [-0.5, 0.2] X = np.arange(0, 4, 0.25) Y = np.arange(0, 4,

最新推荐

recommend-type

在树莓派4B安装 scipy 笔记,不需要删除numpy,不需要mkl

在树莓派4B安装 scipy 笔记,不需要删除numpy,不需要mkl 参考官网: 不要用sudo, 带上–user,否则有问题,官网 最好用pip安装 python -m pip install -i http://mirrors.aliyun.com/pypi/simple/ --user numpy ...
recommend-type

sympy-docs-pdf-1.6.2.pdf

sympy是python的一个用于科学计算的库,可解方程和方程组,计算极限,积分,微分,微积分,矩阵,级数等
recommend-type

基于网络的入侵检测系统源码+数据集+详细文档(高分毕业设计).zip

基于网络的入侵检测系统源码+数据集+详细文档(高分毕业设计).zip个人经导师指导并认可通过的高分毕业设计项目,评审分98分。主要针对计算机相关专业的正在做毕设的学生和需要项目实战练习的学习者,也可作为课程设计、期末大作业。 基于网络的入侵检测系统源码+数据集+详细文档(高分毕业设计).zip基于网络的入侵检测系统源码+数据集+详细文档(高分毕业设计).zip基于网络的入侵检测系统源码+数据集+详细文档(高分毕业设计).zip基于网络的入侵检测系统源码+数据集+详细文档(高分毕业设计).zip基于网络的入侵检测系统源码+数据集+详细文档(高分毕业设计).zip基于网络的入侵检测系统源码+数据集+详细文档(高分毕业设计).zip基于网络的入侵检测系统源码+数据集+详细文档(高分毕业设计).zip基于网络的入侵检测系统源码+数据集+详细文档(高分毕业设计).zip基于网络的入侵检测系统源码+数据集+详细文档(高分毕业设计).zip基于网络的入侵检测系统源码+数据集+详细文档(高分毕业设计).zip基于网络的入侵检测系统源码+数据集+详细文档(高分毕业设计).zip基于网络的入侵检测系统
recommend-type

本户型为2层独栋别墅D026-两层-13.14&12.84米-施工图.dwg

本户型为2层独栋别墅,建筑面积239平方米,占地面积155平米;一层建筑面积155平方米,设有客厅、餐厅、厨房、卧室3间、卫生间1间、杂物间;二层建筑面积84平方米,设有卧室2间、卫生间1间、储藏间、1个大露台。 本户型外观造型别致大方,采光通风良好,色彩明快,整体平面布局紧凑、功能分区合理,房间尺度设计适宜,豪华大气,富有时代气息。
recommend-type

Java_带有可选web的开源命令行RatioMaster.zip

Java_带有可选web的开源命令行RatioMaster
recommend-type

zigbee-cluster-library-specification

最新的zigbee-cluster-library-specification说明文档。
recommend-type

管理建模和仿真的文件

管理Boualem Benatallah引用此版本:布阿利姆·贝纳塔拉。管理建模和仿真。约瑟夫-傅立叶大学-格勒诺布尔第一大学,1996年。法语。NNT:电话:00345357HAL ID:电话:00345357https://theses.hal.science/tel-003453572008年12月9日提交HAL是一个多学科的开放存取档案馆,用于存放和传播科学研究论文,无论它们是否被公开。论文可以来自法国或国外的教学和研究机构,也可以来自公共或私人研究中心。L’archive ouverte pluridisciplinaire
recommend-type

【实战演练】MATLAB用遗传算法改进粒子群GA-PSO算法

![MATLAB智能算法合集](https://static.fuxi.netease.com/fuxi-official/web/20221101/83f465753fd49c41536a5640367d4340.jpg) # 2.1 遗传算法的原理和实现 遗传算法(GA)是一种受生物进化过程启发的优化算法。它通过模拟自然选择和遗传机制来搜索最优解。 **2.1.1 遗传算法的编码和解码** 编码是将问题空间中的解表示为二进制字符串或其他数据结构的过程。解码是将编码的解转换为问题空间中的实际解的过程。常见的编码方法包括二进制编码、实数编码和树形编码。 **2.1.2 遗传算法的交叉和
recommend-type

openstack的20种接口有哪些

以下是OpenStack的20种API接口: 1. Identity (Keystone) API 2. Compute (Nova) API 3. Networking (Neutron) API 4. Block Storage (Cinder) API 5. Object Storage (Swift) API 6. Image (Glance) API 7. Telemetry (Ceilometer) API 8. Orchestration (Heat) API 9. Database (Trove) API 10. Bare Metal (Ironic) API 11. DNS
recommend-type

JSBSim Reference Manual

JSBSim参考手册,其中包含JSBSim简介,JSBSim配置文件xml的编写语法,编程手册以及一些应用实例等。其中有部分内容还没有写完,估计有生之年很难看到完整版了,但是内容还是很有参考价值的。