龙格-库塔法里的能量漂移

问题描述

我们用最常用的数值求解方法—龙格库塔法求解了平面三体问题的运动方程；但发现由于误差的积累效应，整个体系不保持能量守恒。

问题状态

已经解决

进一步解释

我们以 Pythogras 三体问题为例，用 RK4 求解，然后观察体系的能量。



In [1]:

    
import numpy as np
import wq.core.physics.unit.au as au

from math import sqrt
from wq.core.math.ode import rk4 as solver
from wq.core.physics.nbody.body3p import derivativeOf

deriv = derivativeOf(au, 5.0, 3.0, 4.0)
step = solver(deriv)

time = 0
phase = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 3.0, 0.0, 0.0, 0.0])

K = 0.0
U = - au.G * (5.0 * 4.0 / 3.0 + 3.0 * 5.0 / 4.0 + 3.0 * 4.0 / 5.0)
E = K + U

def evolve(tao):
    global time, phase, K, U, E

    time, phase = step(time, phase, tao)
    x1  = phase[0]
    y1  = phase[1]
    vx1 = phase[2]
    vy1 = phase[3]
    x2  = phase[4]
    y2  = phase[5]
    vx2 = phase[6]
    vy2 = phase[7]
    x3  = phase[8]
    y3  = phase[9]
    vx3 = phase[10]
    vy3 = phase[11]

    r12 = sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2))
    r13 = sqrt((x1 - x3) * (x1 - x3) + (y1 - y3) * (y1 - y3))
    r23 = sqrt((x2 - x3) * (x2 - x3) + (y2 - y3) * (y2 - y3))

    v1 = sqrt(vx1 * vx1 + vy1 * vy1)
    v2 = sqrt(vx2 * vx2 + vy2 * vy2)
    v3 = sqrt(vx3 * vx3 + vy3 * vy3)

    K = (5.0 * (vx1 * vx1 + vy1 * vy1) + 3.0 * (vx2 * vx2 + vy2 * vy2) + 4.0 * (vx3 * vx3 + vy3 * vy3)) /2
    U = - au.G * (5.0 * 3.0 / r12 + 5.0 * 4.0 / r13 + 3.0 * 4.0 / r23)
    E = K + U

    return min(r12, r13, r23), max(v1, v2, v3)

演化程序做好了，让我们做一个图表看看星体的轨道演进和能量的变化趋势。



In [2]:

    
import plotly.plotly as py
from plotly.graph_objs import *

serial_time = [time]
serial_x1 = [0]
serial_y1 = [0]
serial_x2 = [0]
serial_y2 = [4]
serial_x3 = [3]
serial_y3 = [0]
serial_energy = [E]

dt = 0.1
for i in range(5000):
    minDist, maxVelo = evolve(dt)
    dt = min(minDist / maxVelo / 100, 0.1)

    x1 = phase[0]
    y1 = phase[1]
    x2 = phase[4]
    y2 = phase[5]
    x3 = phase[8]
    y3 = phase[9]

    serial_time.append(time)
    serial_x1.append(x1)
    serial_y1.append(y1)
    serial_x2.append(x2)
    serial_y2.append(y2)
    serial_x3.append(x3)
    serial_y3.append(y3)
    serial_energy.append(E)

trace = Scatter(x=serial_time, y=serial_energy)
data = Data([trace])
unique_url = py.plot(data, filename = 'energy-drift-by-rk4')

trace1 = Scatter(x=serial_x1, y=serial_y1)
trace2 = Scatter(x=serial_x2, y=serial_y2)
trace3 = Scatter(x=serial_x3, y=serial_y3)

data = Data([trace1, trace2, trace3])
unique_url = py.plot(data, filename = 'pythogras-3-body-rk4')









    



High five! You successfuly sent some data to your account on plotly. View your plot in your browser at https://plot.ly/~mountain/0 or inside your plot.ly account where it is named 'energy-drift-by-rk4'

龙格-库塔法里的能量漂移

问题描述

问题状态

进一步解释

相关问题