logo资料库

四阶龙格库塔法模拟二维三体运动.docx

发布时间：2022-06-15 发布人：admin 分类：说明书 资料大小：0.22M 资料格式：docx 举报 版权申诉
第1页 / 共10页
第2页 / 共10页
第3页 / 共10页
第4页 / 共10页
第5页 / 共10页
第6页 / 共10页
第7页 / 共10页
第8页 / 共10页
资料共10页，剩余部分请下载后查看
    数蕴物语之数据可视化大赛
    目 录
    (一)、呈现效果图
    (二)、实现方式
    1.使用软件：
    2.使用代码：见附录
    (三)、数据来源
    (四)、数据处理方式
    (五)、数据内涵
    (六)、对本次实验的心得与体会
    (七)、对本次活动的建议
    (八)、附录
    数蕴物语之数据可视化大赛 队名： 队长学号： 队长学院： 联系方式：
    玫瑰线的设计与实现 目 录 (一)、 呈现效果图............................................................................... 3 (二)、 实现方式....................................................................................4 1. 使用软件：..........................................................................................4 2. 使用代码：..........................................................................................4 (三)、 数据来源....................................................................................4 (四)、 数据处理方式........................................................................... 4 (五)、 数据内涵....................................................................................4 (六)、 对本次实验的心得与体会...................................................... 4 (七)、 对本次活动的建议...................................................................4 (八)、 附录............................................................................................5 2 / 10
    玫瑰线的设计与实现 (一)、 呈现效果图 不受限制混沌的三体运动 星体 1（红线）坐标（3,0）星体 2（绿线）（0,0）星体 3（蓝线）（-3,0） 质量 1 速度（0,4） 质量 3 速度（0,0） 质量 1 速度（0，-4） 受限制三体运动的一种较稳定的解 3 / 10
    玫瑰线的设计与实现 (二)、 实现方式 1. 使用软件： Python3.6 2. 使用代码：见附录 (三)、 数据来源 由牛顿第二定律 F=ma 和万有引力定律 F= 1m mG 2 2 r 得出三个质点相互作用下的三个二阶常 微分运动方程组（约束在平面），通过初始条件的改变，求解得到不同的运动轨迹。 (四)、 数据处理方式 通过引入速度 v，将一个二阶常微分方程（x’’=a）化为含两个一阶常微分方程的方程组（x’=v， v’=a），将每个质点的运动分为 x，y 两个方向从而得到一个含 12 个一阶常微分方程的方程 组。 编程实现 4 阶龙格-库塔方法（RK4）求解得到的一阶常微分方程组，得到质点组在 h=Δt 之 后的位置，把这些点连起来就得到三个质点在相互作用下的运动轨迹。由于万有引力常数 G=6.67*10^-11,模拟时不可能使用这么小的数计算，所以本次模拟设定 G=10,也能得到比较 准确的天体运动模拟图像。 (五)、 数据内涵 设定三个质点的质量和所在的位置和速度即初始条件，可以通过编程计算得到从 t=0 开始， 每隔一小段时间 h=Δt 时刻三个质点位置，从而通过控制 h 足够小来实现三体运动的轨迹绘 制。 (六)、 对本次实验的心得与体会 通过使用 Python 绘制该图形，让我更加熟悉了数学软件的基本运用。以后可以籍此在别的 方面得到一定的帮助。同时了解了三体问题虽无解析解，但一定程度上可以模拟其运动轨迹。 (七)、 对本次活动的建议 1. 建议将编程的文件和代码一并上交。 2. 建议去掉“对本次活动的建议”一项。 3. 增加赛前交流会。 4 / 10
    玫瑰线的设计与实现 (八)、 附录 import numpy as np import matplotlib.pyplot as plt def m(i, c, n, p): #i 为时间间隔，c 为初值条件，n 为分割次数，p 为质量，龙格库塔方法近 似 h = (i[1] - i[0])/n x1 = c[0] vx1 = c[1] x2 = c[2] vx2 = c[3] x3 = c[4] vx3 = c[5] y1 = c[6] vy1 = c[7] y2 = c[8] vy2 = c[9] y3 = c[10] vy3 = c[11] m1 = p[0] m2 = p[1] m3 = p[2] plx1 = np.empty((n,), dtype=np.float) ply1 = np.empty((n,), dtype=np.float) plx2 = np.empty((n,), dtype=np.float) ply2 = np.empty((n,), dtype=np.float) plx3 = np.empty((n,), dtype=np.float) ply3 = np.empty((n,), dtype=np.float) g = 10 for t in range(n):#梯形法则 x1k1 = vx1 y1k1 = vy1 x2k1 = vx2 y2k1 = vy2 x3k1 = vx3 y3k1 = vy3 x1l1 = g*m2*(x2-x1)/pow(pow((x2-x1), 2)+pow((y2-y1), 1.5)+g*m3*(x3-x1)/pow(pow(x3-x1, 2)+pow(y3-y1, 2), 1.5) y1l1 = g*m2*(y2-y1)/pow(pow((x2-x1), 2)+pow((y2-y1), 1.5)+g*m3*(y3-y1)/pow(pow(x3-x1, 2)+pow(y3-y1, 2), 1.5) x2l1 = g*m3*(x3-x2)/pow(pow((x3-x2), 2)+pow((y3-y2), 1.5)+g*m1*(x1-x2)/pow(pow(x1-x2, 2)+pow(y1-y2, 2), 1.5) 2), 2), 2), 5 / 10
    玫瑰线的设计与实现 y2l1 = g*m3*(y3-y2)/pow(pow((x3-x2), 2)+pow((y3-y2), 1.5)+g*m1*(y1-y2)/pow(pow(x1-x2, 2)+pow(y1-y2, 2), 1.5) x3l1 = g*m1*(x1-x3)/pow(pow((x1-x3), 2)+pow((y1-y3), 1.5)+g*m2*(x2-x3)/pow(pow(x2-x3, 2)+pow(y1-y3, 2), 1.5) y3l1 = g*m1*(y1-y3)/pow(pow((x1-x3), 2)+pow((y1-y3), 2), 2), 2), 1.5)+g*m2*(y2-y3)/pow(pow(x2-x3, 2)+pow(y2-y3, 2), 1.5) x1k2 = x1k1 + h*x1l1/2 y1k2 = y1k1 + h*y1l1/2 x2k2 = x2k1 + h*x2l1/2 y2k2 = y2k1 + h*y2l1/2 x3k2 = x3k1 + h*x3l1/2 y3k2 = y3k1 + h*y3l1/2 x1l2 = g*m2*(x2+x2k1*h/2-x1-x1k1*h/2)/pow(pow((x2+x2k1*h/2-x1-x1k1*h/2), 2)+pow((y2+y2k1*h/2-y1-y1k1*h/2), 2), 1.5)+g*m3*(x3+x3k1*h/2-x1-x1k1*h/2)/pow(pow(x3+x3k1*h/2-x1-x1k1*h/2, 2)+pow(y3+y3k1*h/2-y1-y1k1 * h / 2, 2), 1.5) y1l2 = g*m2*(y2+y2k1*h/2-y1-y1k1*h/2)/pow(pow((x2+x2k1*h/2-x1-x1k1*h/2), 2)+pow((y2+y2k1*h/2-y1-y1k1*h/2), 2), 1.5)+g*m3*(y3+y3k1*h/2-y1-y1k1*h/2)/pow(pow(x3+x3k1*h/2-x1-x1k1*h/2, 2)+pow(y3+y3k1*h/2-y1-y1k1 * h / 2, 2), 1.5) x2l2 = g*m3*(x3+x3k1*h/2-x2-x2k1*h/2)/pow(pow((x3+x3k1*h/2-x2-x2k1*h/2), 2)+pow((y3+y3k1*h/2-y2-y2k1*h/2), 2), 1.5)+g*m1*(x1+x1k1*h/2-x2-x2k1*h/2)/pow(pow(x1+x1k1*h/2-x2-x2k1*h/2, 2)+pow(y1+y1k1*h/2-y2-y2k1 * h / 2, 2), 1.5) y2l2 = g*m3*(y3+y3k1*h/2-y2-y2k1*h/2)/pow(pow((x3+x3k1*h/2-x2-x2k1*h/2), 2)+pow((y3+y3k1*h/2-y2-y2k1*h/2), 2), 1.5)+g*m1*(y1+y1k1*h/2-y2-y2k1*h/2)/pow(pow(x1+x1k1*h/2-x2-x2k1*h/2, 2)+pow(y1+y1k1*h/2-y2-y2k1 * h / 2, 2), 1.5) 6 / 10
    玫瑰线的设计与实现 x3l2 = g*m1*(x1+x1k1*h/2-x3-x3k1*h/2)/pow(pow((x1+x1k1*h/2-x3-x3k1*h/2), 2)+pow((y1+y1k1*h/2-y3-y3k1*h/2), 2), 1.5)+g*m2*(x2+x2k1*h/2-x3-x3k1*h/2)/pow(pow(x2+x2k1*h/2-x3-x3k1*h/2, 2)+pow(y1+y1k1*h/2-y3-y3k1 * h / 2, 2), 1.5) y3l2 = g*m1*(y1+y1k1*h/2-y3-y3k1*h/2)/pow(pow((x1+x1k1*h/2-x3-x3k1*h/2), 2)+pow((y1+y1k1*h/2-y3-y3k1*h/2), 2), 1.5)+g*m2*(y2+y2k1*h/2-y3-y3k1*h/2)/pow(pow(x2+x2k1*h/2-x3-x3k1*h/2, 2)+pow(y2+y2k1*h/2-y3-y3k1 * h / 2, 2), 1.5) x1k3 = x1k1 + h*x1l2/2 y1k3 = y1k1 + h*y1l2/2 x2k3 = x2k1 + h*x2l2/2 y2k3 = y2k1 + h*y2l2/2 x3k3 = x3k1 + h*x3l2/2 y3k3 = y3k1 + h*y3l2/2 x1l3 = g*m2*(x2+x2k2*h/2-x1-x1k2*h/2)/pow(pow((x2+x2k2*h/2-x1-x1k2*h/2), 2)+pow((y2+y2k2*h/2-y1-y1k2*h/2), 2), 1.5)+g*m3*(x3+x3k2*h/2-x1-x1k2*h/2)/pow(pow(x3+x3k2*h/2-x1-x1k2*h/2, 2)+pow(y3+y3k2*h/2-y1-y1k2 * h / 2, 2), 1.5) y1l3 = g*m2*(y2+y2k2*h/2-y1-y1k2*h/2)/pow(pow((x2+x2k2*h/2-x1-x1k2*h/2), 2)+pow((y2+y2k2*h/2-y1-y1k2*h/2), 2), 1.5)+g*m3*(y3+y3k2*h/2-y1-y1k2*h/2)/pow(pow(x3+x3k2*h/2-x1-x1k2*h/2, 2)+pow(y3+y3k2*h/2-y1-y1k2 * h / 2, 2), 1.5) x2l3 = g*m3*(x3+x3k2*h/2-x2-x2k2*h/2)/pow(pow((x3+x3k2*h/2-x2-x2k2*h/2), 2)+pow((y3+y3k2*h/2-y2-y2k2*h/2), 2), 1.5)+g*m1*(x1+x1k2*h/2-x2-x2k2*h/2)/pow(pow(x1+x1k2*h/2-x2-x2k2*h/2, 2)+pow(y1+y1k2*h/2-y2-y2k2 7 / 10
    玫瑰线的设计与实现 * h / 2, 2), 1.5) y2l3 = g*m3*(y3+y3k2*h/2-y2-y2k2*h/2)/pow(pow((x3+x3k2*h/2-x2-x2k2*h/2), 2)+pow((y3+y3k2*h/2-y2-y2k2*h/2), 2), 1.5)+g*m1*(y1+y1k2*h/2-y2-y2k2*h/2)/pow(pow(x1+x1k2*h/2-x2-x2k2*h/2, 2)+pow(y1+y1k2*h/2-y2-y2k2 * h / 2, 2), 1.5) x3l3 = g*m1*(x1+x1k2*h/2-x3-x3k2*h/2)/pow(pow((x1+x1k2*h/2-x3-x3k2*h/2), 2)+pow((y1+y1k2*h/2-y3-y3k2*h/2), 2), 1.5)+g*m2*(x2+x2k2*h/2-x3-x3k2*h/2)/pow(pow(x2+x2k2*h/2-x3-x3k2*h/2, 2)+pow(y1+y1k2*h/2-y3-y3k2 * h / 2, 2), 1.5) y3l3 = g*m1*(y1+y1k2*h/2-y3-y3k2*h/2)/pow(pow((x1+x1k2*h/2-x3-x3k2*h/2), 2)+pow((y1+y1k2*h/2-y3-y3k2*h/2), 2), 1.5)+g*m2*(y2+y2k2*h/2-y3-y3k2*h/2)/pow(pow(x2+x2k2*h/2-x3-x3k2*h/2, 2)+pow(y2+y2k2*h/2-y3-y3k2 * h / 2, 2), 1.5) x1k4 = x1k1 + h*x1l3 y1k4 = y1k1 + h*y1l3 x2k4 = x2k1 + h*x2l3 y2k4 = y2k1 + h*y2l3 x3k4 = x3k1 + h*x3l3 y3k4 = y3k1 + h*y3l3 x1l4 = g*m2*(x2+x2k3*h-x1-x1k3*h)/pow(pow((x2+x2k3*h-x1-x1k3*h), 2)+pow((y2+y2k3*h-y1-y1k3*h), 2), 1.5)+g*m3*(x3 + x3k3*h-x1-x1k3*h)/pow(pow(x3+x3k3*h-x1-x1k3*h, 2)+pow(y3+y3k3*h-y1-y1k3 * h, 2), 1.5) y1l4 = g*m2*(y2+y2k3*h-y1-y1k3*h)/pow(pow((x2+x2k3*h-x1-x1k3*h), 2)+pow((y2+y2k3*h-y1-y1k3*h), 2), 1.5)+g*m3*(y3 + y3k3*h-y1-y1k3*h)/pow(pow(x3+x3k3*h-x1-x1k3*h, 2)+pow(y3+y3k3*h-y1-y1k3 * h, 2), 1.5) x2l4 = g*m3*(x3+x3k3*h-x2-x2k3*h)/pow(pow((x3+x3k3*h-x2-x2k3*h), 2)+pow((y3+y3k3*h-y2-y2k3*h), 2), 1.5)+g*m1*(x1 + x1k3*h-x2-x2k3*h)/pow(pow(x1+x1k3*h-x2-x2k3*h, 2)+pow(y1+y1k3*h-y2-y2k3 * h, 2), 1.5) 8 / 10
    分享到：
    赞 收藏

    相关推荐

    资料库

    后端

    共收录2142份资料，累计10个分类，关注成员有19位，主要包括：Perl,Docker,Delphi,VB,PHP,Java,C,C#,Python,C++

    热门标签

    Perl Docker Delphi VB PHP Java C C# Python C++

    最新资料