用户注册 登录
珍珠湾全球网 返回首页

岳东晓 -- 珍珠湾全球网 ... https://ydx.zzwave.com [收藏] [复制] [分享] [RSS] 岳东晓 -- 珍珠湾全球网

日志

太阳系模拟(含哈雷彗星)-- ChatGPT 编程

已有 6203 次阅读2023-9-8 08:15 |个人分类:科普|系统分类:科普

solar-chat-gpt——使用Clipchamp制作.mp4

1677年,牛顿用拉丁语写了一本 《Philosophiæ Naturalis Principia Mathematica》,里面充斥着各种几何图形,没有用到任何微积分。为了解决一个现在看来简单的力学问题,牛顿必须证明各种几何定理、引理。100年后的1788 年,法国(意大利裔)数学家拉格朗日写了一本《Mécanique Analytique》(《分析力学》),他自豪地说,全书没有一张图。解决牛顿力学问题往往需要相当好的物理概念,要知道向量等概念,要进行受力分析之类。但使用分析力学方法,你只需要写下一个称为拉格朗日量的标量,剩下的就是计算。后人在其基础上又发展出哈密顿力学。你只要写出系统的能量 H,接下来就是两个一阶方程, q' = dH/dp, p' = -dH/dq,几乎不要动脑筋,硬算即可推导出系统的演化。

ChatGPT 自从推出以来,其能力水平大有提高,但我前天测试一下,它还是不能正确解决那个物体从斜面滑下、斜面后退、请计算加速度的中学力学题。但我突然想起何不让它用哈密顿力学算一下太阳系的运动。ChatGPT 看到我的提示,毫不犹豫的开始写代码。一次运行畅通无阻,但我发现地球的轨道不对,一年没有回到起点。我查看了一下,ChatGPT 在算星体之间的引力能时算了两次,A与B,然后是 B与A。于是我提示到 H 里 V 算重了,就这么个提示,ChatGPT 立刻知道了问题所在,而且立刻修改了代码。这不得不令人赞叹。(代码附在文末)

视频是太阳系运转图。从水星到火星,到太阳的距离比例是对的。但从木星开始,距离显示进行了比例缩减。那个白色的是哈雷彗星。




屏幕截图 2023-09-07 165053.png

import sympy as sp

# The Hamiltonian system definition as above...

class HamiltonianSystem:
    def __init__(self, N):
        self.xs = sp.symbols(f'x0:{N}')
        self.ys = sp.symbols(f'y0:{N}')
        self.pxs = sp.symbols(f'px0:{N}')
        self.pys = sp.symbols(f'py0:{N}')
        self.ms = sp.symbols(f'm0:{N}')
        
        T = sum(self.pxs[i]**2 / (2 * self.ms[i]) + self.pys[i]**2 / (2 * self.ms[i]) for i in range(N))
        V = sum(-self.ms[i] * self.ms[j] / sp.sqrt((self.xs[i] - self.xs[j])**2 + (self.ys[i] - self.ys[j])**2)
                for i in range(N) for j in range(i+1, N) if i != j)
    
        self.H = T + V
   
        self.H_numeric = lambdify((*self.xs, *self.ys, *self.pxs, *self.pys, *self.ms), self.H)

        
    def dHdx(self, i):
        return sp.diff(self.H, self.xs[i])
    
    def dHdy(self, i):
        return sp.diff(self.H, self.ys[i])
    
    def dHdpx(self, i):
        return sp.diff(self.H, self.pxs[i])
    
    def dHdpy(self, i):
        return sp.diff(self.H, self.pys[i])

    def H_value(self, states, masses):
        N = len(states)
        
        # Separate states for easy referencing
        x = [states[i][0] for i in range(N)]
        px = [states[i][1] for i in range(N)]
        y = [states[i][2] for i in range(N)]
        py = [states[i][3] for i in range(N)]
        
        # Use the lambdified function
        return self.H_numeric(*x, *y, *px, *py, *masses)


# The numerical simulation as above...

from sympy import lambdify

class NumericalSimulation:
    def __init__(self, N, masses):
        self.system = HamiltonianSystem(N)
        self.masses = masses
        
        self.dHdx_funcs = [lambdify((*self.system.xs, *self.system.ys, *self.system.pxs, *self.system.pys, *self.system.ms), 
                                    self.system.dHdx(i)) for i in range(N)]
        self.dHdy_funcs = [lambdify((*self.system.xs, *self.system.ys, *self.system.pxs, *self.system.pys, *self.system.ms), 
                                    self.system.dHdy(i)) for i in range(N)]
        self.dHdpx_funcs = [lambdify((*self.system.xs, *self.system.ys, *self.system.pxs, *self.system.pys, *self.system.ms), 
                                     self.system.dHdpx(i)) for i in range(N)]
        self.dHdpy_funcs = [lambdify((*self.system.xs, *self.system.ys, *self.system.pxs, *self.system.pys, *self.system.ms), 
                                     self.system.dHdpy(i)) for i in range(N)]

    def step(self, states, dt):
        N = len(states)
        
        # Separate states for easy referencing
        x = [states[i][0] for i in range(N)]
        px = [states[i][1] for i in range(N)]
        y = [states[i][2] for i in range(N)]
        py = [states[i][3] for i in range(N)]
        
        for i in range(N):
            px[i] -= 0.5 * dt * self.dHdx_funcs[i](*x, *y, *px, *py, *self.masses)
            py[i] -= 0.5 * dt * self.dHdy_funcs[i](*x, *y, *px, *py, *self.masses)
            
        for i in range(N):
            x[i] += dt * px[i] / self.masses[i]
            y[i] += dt * py[i] / self.masses[i]
            
        for i in range(N):
            px[i] -= 0.5 * dt * self.dHdx_funcs[i](*x, *y, *px, *py, *self.masses)
            py[i] -= 0.5 * dt * self.dHdy_funcs[i](*x, *y, *px, *py, *self.masses)
            
        # Repack the states
        new_states = [(x[i], px[i], y[i], py[i]) for i in range(N)]
        
        return new_states
    
    def energy(self, states):
        return self.system.H_value(states, self.masses)


import numpy as np

class SymplecticIntegrator:

    def __init__(self, dHdx_funcs, dHdy_funcs, masses):
        self.dHdx_funcs = dHdx_funcs
        self.dHdy_funcs = dHdy_funcs
        self.masses = masses
        
        # Coefficients for 2nd order method
        self.c2 = [0, 1]
        self.d2 = [0.5, 0.5]

        # Coefficients for 3rd order method
        self.c3 = [1, -2/3, 2/3]
        self.d3 = [-1/24, 3/4, 7/24]

        # Coefficients for 4th order Yoshida method
        t13 = np.power(2, 1/3)
        self.c4 = [1 /(2* (2-t13)),
                   (1-t13)/2/(2-t13),
                   (1-t13)/2/(2-t13),
                   1 /(2* (2-t13))]
        self.d4 = [1 / (2-t13),
                   -t13/(2-t13),
                   1 / (2-t13),
                   0]

    def position_step(self, x, px, y, py, dt, c):
        N = len(x)
        for i in range(N):
            x[i] += dt * c * px[i] / self.masses[i]
            y[i] += dt * c * py[i] / self.masses[i]
        return x, y

    def momentum_step(self, x, px, y, py, dt, d):
        N = len(x)
        for i in range(N):
            px[i] -= dt * d * self.dHdx_funcs[i](*x, *y, *px, *py, *self.masses)
            py[i] -= dt * d * self.dHdy_funcs[i](*x, *y, *px, *py, *self.masses)
        return px, py

    def step(self, states, dt, order=2):
        N = len(states)
        x = [states[i][0] for i in range(N)]
        px = [states[i][1] for i in range(N)]
        y = [states[i][2] for i in range(N)]
        py = [states[i][3] for i in range(N)]
        
        if order == 2:
            coeffs = zip(self.c2, self.d2)
        elif order == 3:
            coeffs = zip(self.c3, self.d3)
        elif order == 4:
            coeffs = zip(self.c4, self.d4)
        else:
            raise NotImplementedError(f"Order {order} not implemented.")
            
        for c, d in coeffs:
            x, y = self.position_step(x, px, y, py, dt, c)
            px, py = self.momentum_step(x, px, y, py, dt, d)
            
        new_states = [(x[i], px[i], y[i], py[i]) for i in range(N)]
        
        return new_states


import pygame
import sys
import math
import os

# Test code as requested...

v_earth = 29.78 #km/s

G = 6.67430e-11  # gravitational constant in m^3 kg^−1 s^−2
M = 1.989e30  # mass of the Sun in kg
r = 0.586 * 1.496e11  # distance from the Sun at perihelion in meters
a = 17.83 * 1.496e11  # semi-major axis in meters

v_haley = math.sqrt(G * M * (2/r - 1/a))/v_earth/1000

print(f"vh={v_haley}")

# Normalized masses based on the mass of the Sun set to 1
mass_sun = 1
mass_earth = 3e-6
mass_mercury = 1.66e-7
mass_venus = 2.45e-6
mass_mars = 3.3e-7
mass_jupiter = 9.55e-4
mass_saturn = 2.86e-4
mass_uranus = 4.37e-5
mass_neptune = 5.15e-5
mass_haley = 1e-8

masses = [mass_sun, mass_mercury, mass_venus,mass_earth,  mass_mars, mass_jupiter, mass_saturn, mass_uranus, mass_neptune, mass_haley]


# Normalized distances from the sun (in AU) and velocities with Earth's set to 1

earth = (1, 0, 0, mass_earth) # Earth's velocity is the reference (momentum in y-direction)
mercury = (0.39, 0, 0, mass_mercury * 1.61) # Orbital velocity normalized to Earth's
venus = (0.723, 0, 0, mass_venus * 1.18)
mars = (1.524, 0, 0, mass_mars * 0.81)
jupiter = (5.203, 0, 0, mass_jupiter * 0.44)
saturn = (9.537, 0, 0, mass_saturn * 0.33)
uranus = (19.191, 0, 0, mass_uranus * 0.23)
neptune = (30.069, 0, 0, mass_neptune * 0.18)
haley = (0.586, 0,0, v_haley*mass_haley)

sun_p= - sum(s[3] for s in [mercury, venus,earth, mars, jupiter, saturn, uranus, neptune, haley])

sun = (0, 0, 0, sun_p) # Sun stays stationary

states = [sun,  mercury, venus,earth, mars, jupiter, saturn, uranus, neptune, haley]


sim = NumericalSimulation(len(masses), masses)

stepper = SymplecticIntegrator(sim.dHdx_funcs, sim.dHdy_funcs,masses)



dt = 0.01
num_steps = 628  # For example, this evolves the system for 10 units of time

# for _ in range(num_steps):
#     state = sim.step(states, dt)
#     total_energy = sim.system.H_value(states, masses)
#     print(f"Total energy of the system: {total_energy}")
#     #print(f"Position of object 1: ({x[0]}, {y[0]})")
#     print(f"Position of object 2: ({state[1][0]}, {state[1][2]})")



# Initialize pygame
pygame.init()

# Screen settings
WIDTH, HEIGHT = 1800, 1800

screen = pygame.display.set_mode((WIDTH, HEIGHT))
pygame.display.set_caption("Solar System Visualization")
clock = pygame.time.Clock()
# Colors
WHITE = (255, 255, 255)
YELLOW = (255, 255, 0)  # Sun
COLORS = [
    (139, 69, 19),  # Mercury: brown
    (245, 123, 0),  # Venus: pink
    (0, 0, 255),  # Earth: blue
    (255, 0, 0),  # Mars: red
    (255, 165, 0),  # Jupiter: orange
    (210, 105, 30),  # Saturn: chocolate
    (64, 224, 208),  # Uranus: turquoise
    (0, 0, 128),  # Neptune: navy
    (125,125,125)
]

# Sizes
SUN_SIZE = 30
PLANET_SIZES = [8, 12, 15, 12, 25, 20, 18, 17,18]  # Adjust as needed to visualize comfortably
GRADIENT_RAYS=30

pygame.font.init()

font = pygame.font.SysFont(None, 20)

def display_energy(screen, steps, energy):
    # Initialize the font system (can be done once at the start of the program)
   

    # Render the text
    text_surface = font.render(f"Total Energy: {energy:.2f}", True, (255, 0, 255))  # White color
    # Display the text at the top left corner
    screen.blit(text_surface, (10, 10))
    text_surface = font.render(f"Years: {steps/628.0:0.2f}; Steps: {steps}", True, (255, 0, 255))  # White color
    # Display the text at the top left corner
    screen.blit(text_surface, (10, 32))
   


def draw_star(screen, color, center, radius, edge=1):
    pygame.draw.circle(screen, color, center, radius)

    edge_color = [int(c/2) for c in color]

    # Draw gradient rings
    for i in range(1, GRADIENT_RAYS + 1):
        alpha_value = 255 - i * (255 // GRADIENT_RAYS)
        color_with_alpha = tuple(edge_color) + (alpha_value,)

        # Create a temporary surface for this specific ring
        ring_radius = radius + i * edge
        ring_thickness = edge
        temp_surface = pygame.Surface((WIDTH, HEIGHT), pygame.SRCALPHA)
        
        # Draw the outer circle of the ring
        pygame.draw.circle(temp_surface, color_with_alpha, center, ring_radius)
        
        # Mask out the inner part to only leave the fading edge (ring) visible
        mask_color = (0, 0, 0, 0)  # fully transparent
        pygame.draw.circle(temp_surface, mask_color, center, ring_radius - ring_thickness)

        screen.blit(temp_surface, (0, 0))


# Positions (here we're just starting them at their initial distances, you'll want to update this dynamically based on your simulation)

LOG_BASE = 1.5

running = True

scale = WIDTH/5/2

steps=0
while running:
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            running = False

    screen.fill((0,0,0))

    # Drawing Sun
    states = stepper.step(states, dt,order=3)
    steps +=1
    positions = [(s[0], s[2]) for s in states]
    total_energy = sim.system.H_value(states, masses)
    #print(f"Total energy of the system: {total_energy}")
    draw_star(screen, YELLOW, (positions[0][0]*scale+WIDTH/2, positions[0][1]*scale+HEIGHT/2), SUN_SIZE)

    # Drawing planets
    for idx, position in enumerate(positions[1:]):
        x, y = position
        d = math.sqrt(x*x+y*y)
        #dl = math.log(d+1, LOG_BASE)
        dl=d
        if d>2:
            dl = 2+(d-2)/8 
       
        x = x/d *dl*scale+WIDTH/2
        y = y/d *dl*scale+HEIGHT/2
        pygame.draw.circle(screen, COLORS[idx], (x,y), PLANET_SIZES[idx])

    display_energy(screen, steps, total_energy)
    #imgpath = os.path.join("k:/tmp/solar/", f"{steps:05d}.png")
    #pygame.image.save(screen, imgpath)
    pygame.display.flip()
    clock.tick(100)

    # If you're running a simulation, update the positions here...

pygame.quit()
sys.exit()



路过

鸡蛋

鲜花

支持

雷人

难过

搞笑
 

评论 (0 个评论)

facelist

您需要登录后才可以评论 登录 | 用户注册

Archiver|手机版|珍珠湾全球网

GMT+8, 2024-12-5 18:55 , Processed in 0.010916 second(s), 8 queries , Memcache On.

Powered by Discuz! X2.5

回顶部