Python数值分析案例01--------四阶龙格库塔法解抛体运动

在考虑到空气阻力（ $f=-kv^n$ ）和因地球自转引起的科里奥利力后，列出运动方程

$\begin{bmatrix} 1 & & \\ &1 & \\ & &1 \end{bmatrix}\begin{bmatrix} x''\\ y''\\ z'' \end{bmatrix}=\begin{bmatrix} -2\omega (cos\varphi *z'-sin\varphi *y')-k(x')^n\\ -2\omega sin\varphi *x'-k(y')^n\\ -g+2\omega cos\varphi *x'-k(z')^n \end{bmatrix}$

$\varphi$ ---->抛出点的北纬读数

$\omega$ ---->表示地球自转角速度

算法:

四阶龙格库塔

代码

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from math import *def runge_kutta4(df,a,b,h,y0):num = len(y0)t = np.arange(a,b+h,h)w = np.zeros((t.size,num))w[0,:] = y0for i in range(t.size - 1):s0 = df(t[i], w[i,:])s1 = df(t[i] + h/2.,w[i,:] + h * s0 / 2.)s2 = df(t[i] + h/2.,w[i,:] + h * s1 / 2.)s3 = df(t[i+1], w[i,:] + h * s2)w[i+1,:] = w[i,:] + h * (s0 + 2*(s1+s2) + s3) / 6.if w[i+1,:][2] <= 0:q = ixl,yl,zl =  w[i+1,:][0],w[i+1,:][1],max(w[:,2])breakreturn t,w,q,xl,yl,zldef df(t, variables):x,y,z,x1,y1,z1 = variablesA = np.zeros((3,3))b = np.zeros(3)A[0,0] = 1A[0,1] = 0A[0,2] = 0A[1,0] = 0A[1,1] = 1A[1,2] = 0A[2,0] = 0A[2,1] = 0A[2,2] = 1phi = 30 * pi/180k = 0.00001n = 5g = 9.8b[0] = -2 * pi * 15/180/3600 * (cos(phi) * z1 - sin(phi) * y1) - k * (x1)**nb[1] = -2 * pi * 15/180/3600 * sin(phi) * x1 - k * (y1)**nb[2] = -g + pi * 15/180/3600 * 2 * cos(phi) * x1 - k * (z1)**nx2,y2,z2 = np.linalg.solve(A, b)return np.array([x1,y1,z1,x2,y2,z2])a,b = 0.0,20.0
h = 0.01x = 0
y = 0
z = 0
x1 = 10
y1 = 0
z1 = 70y0 = np.array([x,y,z,x1,y1,z1])t,w,q,xl,yl,zl = runge_kutta4(df, a, b, h, y0)x = w[:,0]
y = w[:,1]
z = w[:,2]
x1 = w[:,3]
y1 = w[:,4]
z1 = w[:,5]if 1:fig = plt.figure()ax = fig.add_subplot(projection='3d')ax.scatter(x,y,z,label="projectile",s=0.5)ax.set_xlabel("X"+" covers "+str(round(xl,2)) +" meters in x")ax.set_ylabel("Y"+" covers "+str(round(yl,2)) +"meters in y")ax.set_zlabel("Z"+" the highest is "+str(round(zl,2)) +" meters in z")ax.set_title("3D projectile draw"+"  end in  "+str(q*h)+"  second")ax.legend()else:#plt.plot(x,z)plt.plot(x,y)
plt.show()

效果图

Vx=10,Vy=10,Vz=20,k=0.000001,n=5

Vx,y=0.5Vz=10 ，k=0.05 n=2

Python数值分析案例01--------四阶龙格库塔法解抛体运动相关推荐

二阶偏微分方程组龙格库塔法_1、经典四阶龙格库塔法解一阶微分方程组
1.经典四阶龙格库塔法解一阶微分方程组陕西科技大学数值计算课程设计任务书理学院信息与计算科学/应用数学专业信息08/数学08 班级学生: 题目:典型数值算法的C++语言程序设计课程 ...
四阶龙格库塔法的基本思想_经典四阶龙格库塔法解一阶微分方程组讲义.doc
1.经典四阶龙格库塔法解一阶微分方程组 1.1运用四阶龙格库塔法解一阶微分方程组算法分析 , 经过循环计算由推得 -- 每个龙格-库塔方法都是由一个合适的泰勒方法推导而来,使得其最终全局误差为,一种 ...
四阶龙格库塔法的基本思想_请问用四阶龙格库塔法解二阶微分方程的思想是什么？...
默认y的单位是弧度 k=1000; t=0:0.001:1; Y=[]; err=1 K=[]; Ymax=[]; xishu=1.01; while err X=[0 0]; k=xishu*k; ...
Python实战案例01
文章目录 1.打印购物小票 2.打印蚂蚁森林植树证书 3.绝对温标 4.身体质量指标 5.计算器 6.猜数字 7.逢 7 拍手游戏 8.打印五子棋棋盘 9.房贷计算器 1.打印购物小票购物小票又称购 ...
python解常微分方程龙格库_excel实现四阶龙格库塔法runge-kutta解二阶常微分方程范例.xls...
excel实现四阶龙格库塔法runge-kutta解二阶常微分方程范例,rungekutta,四阶rungekutta法,rungekuttamatlab,四阶rungekutta,rungekutt ...
四阶龙格库塔法求解一次常微分方程组(python实现)
四阶龙格库塔法求解一次常微分方程组一.前言二.RK4求解方程组的要点 1. 将方程组转化为RK4求解要求的标准形式 2. 注意区分每个方程的独立性三.python实现RK4求解一次常微分方程组 ...
代码检查规则：Python语言案例详解
在之前的文章中代码检查规则:Java语言案例详解学习了Java的检查规则.我们今天将学习<代码检查规则:Python语言案例详解>,内容主要分为两个部分:Python的代码检查规则和Pyt ...
AidLux“换脸”案例源码详解（Python）
"换脸"案例源码详解 (Python) faceswap_gui.py用于换脸,可与facemovie_gui.py身体互换源码(上一篇文章)对照观看打开faceswap_gui ...
AidLux“人像抠图”案例源码详解（Python）
"人像抠图"案例源码详解 (Python) seg_gui_meet.py用于人像抠图导入基础包作用详解构建程序图形化类初始化处理函数(人体抠图应用启动时首先被调用) 程序入 ...
python编程入门与案例详解-Python程序设计案例课堂
Python程序设计案例课堂第Ⅰ篇基础知识 1 揭开Python神秘面纱 1.1 什么是Python 1.2 Python的优点和特性 1.2.1 Python的优点 1.2.2 Python的特 ...

Python数值分析案例01--------四阶龙格库塔法解抛体运动

Python数值分析案例01--------四阶龙格库塔法解抛体运动相关推荐

最新文章

热门文章