物理老师教你学Python语言(下篇)

off999 2025-05-25 14:51 74 浏览 0 评论

下篇：物理建模与综合项目

核心目标：掌握微分方程数值解、面向对象编程和交互式可视化，构建可扩展的物理仿真系统

第7章动态系统模拟

7.1 数值解法与经典力学

案例1：弹簧振子动力学（欧拉法）

import numpy as np  
import matplotlib.pyplot as plt  

# 系统参数  
m = 0.5   # 质量 (kg)  
k = 20    # 劲度系数 (N/m)  
x0 = 0.1  # 初始位移 (m)  
v0 = 0    # 初速度  
dt = 0.01 # 时间步长  

# 初始化数组  
t = np.arange(0, 10, dt)  
x = np.zeros_like(t)  
v = np.zeros_like(t)  
x[0], v[0] = x0, v0  

# 迭代计算  
for i in range(1, len(t)):  
    a = -k * x[i-1] / m  # 牛顿第二定律  
    v[i] = v[i-1] + a * dt  
    x[i] = x[i-1] + v[i] * dt  

# 绘制相空间图 (位移-速度关系)  
plt.plot(x, v, label='数值解')  
plt.xlabel('位移 (m)')  
plt.ylabel('速度 (m/s)')

物理讨论：

对比解析解 x(t) = x0 * cos(√(k/m) * t)
讨论能量守恒：数值解的振幅为何会逐渐增大？（欧拉法的误差积累）

案例2：热传导模拟（有限差分法）

# 金属棒初始温度分布  
L = 1.0       # 棒长 (m)  
N = 100       # 空间分段数  
alpha = 1e-4  # 热扩散系数  
dx = L/N  
dt = 0.1  

# 初始条件：左端100°C，右端0°C  
T = np.zeros(N)  
T[0] = 100  

# 温度场更新  
for _ in range(500):  
    T[1:-1] += alpha * dt / dx**2 * (T[2:] - 2*T[1:-1] + T[:-2])  

plt.plot(np.linspace(0, L, N), T, label=f't={500*dt}s')

7.2 波动与场方程

案例：一维波动方程模拟

# 参数设置  
c = 300  # 波速 (m/s)  
x = np.linspace(0, 10, 200)  
dx = x[1] - x[0]  
dt = dx / (2*c)  # 满足稳定性条件  

# 初始条件：高斯脉冲  
u_prev = np.exp(-(x-3)**2 / 0.2)  
u_current = u_prev.copy()  
u_next = np.zeros_like(x)  

# 时间演化  
for step in range(500):  
    u_next[1:-1] = 2*u_current[1:-1] - u_prev[1:-1] + (c*dt/dx)**2 * (  
        u_current[2:] - 2*u_current[1:-1] + u_current[:-2])  
    u_prev, u_current = u_current, u_next  

# 生成波动传播动画  
from matplotlib.animation import FuncAnimation  
fig, ax = plt.subplots()  
line, = ax.plot(x, u_current)  

def update(frame):  
    # 更新代码...  
    return line,  

ani = FuncAnimation(fig, update, frames=100, interval=50)  
ani.save('wave.gif')  # 用于课堂展示

第8章面向对象物理建模

8.1 粒子系统构建

案例：带电粒子在磁场中的运动

class ChargedParticle:  
    def __init__(self, m, q, pos, vel):  
        self.mass = m          # 质量  
        self.charge = q        # 电荷量  
        self.position = pos    # 位置矢量 [x,y,z]  
        self.velocity = vel    # 速度矢量  

    def lorentz_force(self, E, B):  
        F = self.charge * (E + np.cross(self.velocity, B))  
        return F  

    def update(self, F, dt):  
        a = F / self.mass  
        self.velocity += a * dt  
        self.position += self.velocity * dt  

# 使用示例  
particle = ChargedParticle(m=1e-26, q=1.6e-19,  
                          pos=[0,0,0], vel=[1e4,0,0])  
B_field = np.array([0, 0, 0.5])  # 0.5 T 沿z轴  
E_field = np.array([0,0,0])  

# 模拟轨迹  
trajectory = []  
for _ in range(1000):  
    F = particile.lorentz_force(E_field, B_field)  
    particile.update(F, 1e-8)  
    trajectory.append(particle.position.copy())

物理验证：

绘制轨迹应显示螺旋运动（回旋半径 r = \frac{mv}{qB}）
对比理论周期 T = \frac{2\pi m}{qB}

8.2 多体系统模拟

案例：验证开普勒第三定律

class CelestialBody:  
    G = 6.67e-11  # 引力常量  

    def __init__(self, mass, pos, vel):  
        self.mass = mass  
        self.pos = np.array(pos)  
        self.vel = np.array(vel)  

    def gravity_force(self, other):  
        r_vec = other.pos - self.pos  
        r = np.linalg.norm(r_vec)  
        F = self.G * self.mass * other.mass / r**3 * r_vec  
        return F  

# 创建太阳和地球  
sun = CelestialBody(1.989e30, [0,0], [0,0])  
earth = CelestialBody(5.972e24, [1.496e11,0], [0,2.978e4])  

# 模拟一年的轨道（简化版）  
positions = []  
for _ in range(365):  
    F = sun.gravity_force(earth)  
    earth.vel += F / earth.mass * 86400  # 每天更新一次  
    earth.pos += earth.vel * 86400  
    positions.append(earth.pos.copy())  

# 绘制椭圆轨道并计算周期  
plt.plot(*zip(*positions))  
plt.scatter([0], [0], c='orange', s=100)  # 太阳位置

第9章交互式教学工具开发

9.1 动态参数调节

案例：凸透镜成像规律演示（使用ipywidgets）

from ipywidgets import interact, FloatSlider  

@interact(  
    f=FloatSlider(15, min=10, max=30, step=1, description='焦距 f/cm'),  
    u=FloatSlider(30, min=15, max=60, step=1, description='物距 u/cm')  
)  
def lens_simulation(f, u):  
    v = 1 / (1/f - 1/u) if u > f else '无法成像'  # 透镜公式  
    plt.figure(figsize=(8,4))  
    # 绘制光路图代码...  
    plt.show()

输出效果：

实时拖动滑块观察成像变化（虚像/实像、放大/缩小）
自动标注关键点：焦点、二倍焦距点

9.2 电路仿真工具

案例：RC电路充放电过程

from IPython.display import display  
import ipywidgets as widgets  

R_slider = widgets.FloatSlider(1000, min=500, max=5000, description='R(Ω)')  
C_slider = widgets.FloatSlider(1e-6, min=1e-7, max=1e-5, description='C(F)')  

def update_plot(R, C):  
    t = np.linspace(0, 5*R*C, 100)  
    Vc = 5 * (1 - np.exp(-t/(R*C)))  # 充电曲线  
    plt.plot(t*1e3, Vc, label=f'τ={R*C*1e3:.1f}ms')  

widgets.interact(update_plot, R=R_slider, C=C_slider)