每日算法 — 使用java实现滚球:刚体滚动摩擦与能量守恒路径规划

引言

滚球游戏要求玩家控制小球在复杂地形中滚动,避开障碍物并到达目标点。小球的运动受重力、摩擦力和地形坡度共同影响。本文用Java实现一个二维滚球物理模拟系统,核心算法包括刚体滚动动力学方程数值求解、A*路径搜索在重力场中寻找最优路线,以及能量守恒约束下的运动规划。

核心算法思路

1. 刚体滚动动力学

假设小球为均匀实心球体,质量 m,半径 r,转动惯量 I = 2/5 * m * r²。在无滑动滚动条件下:
– 平动方程:m * a = mg * sin(θ) - f
– 转动方程:I * α = f * r
– 无滑动约束:a = α * r

联立得加速度:a = (5/7) * g * sin(θ)

其中 θ 为坡度角。该结果表明滚动小球的等效加速度小于滑动情况(g * sin(θ)),因为部分重力势能转化为转动动能。

2. 地形表示与碰撞

地形用二维高度函数 h(x) 表示。小球位置 (x, y) 处:
– 地形高度 groundY = h(x)
– 坡度 tan(θ) = dh/dx

小球与地面接触时受到法向支持力和切向摩擦力。若小球腾空(y > h(x) + r),则只受重力做抛体运动。

3. A*路径规划

将连续地形离散化为网格图,每个格子 (i, j) 为一个节点。边权考虑:
– 距离代价:欧几里得距离
– 能量代价:上坡需要额外能量,下坡可加速
– 障碍代价:不可通行区域设为无穷大

启发函数使用欧几里得距离加上高度差估计,保证可采纳性。

4. 能量守恒验证

系统总能量 E = K_trans + K_rot + U = (1/2)mv² + (1/2)Iω² + mgh
无滑动滚动时 v = ωr,故 E = (7/10)mv² + mgh

在理想无摩擦情况下,总能量守恒,可用于验证数值积分精度。

完整Java实现

import java.util.*;

/**
 * 滚球游戏物理模拟与路径规划系统
 * 核心算法:刚体滚动动力学 + A*路径搜索 + 能量守恒验证
 */
public class RollingBallSimulator {

    static final double GRAVITY = 9.8;
    static final double DT = 0.01; // 积分步长

    // 小球状态
    static class BallState {
        double x, y;      // 位置
        double vx, vy;    // 线速度
        double omega;     // 角速度(标量,顺时针为正)
        double mass;
        double radius;
        double momentOfInertia;

        BallState(double x, double y, double mass, double radius) {
            this.x = x;
            this.y = y;
            this.mass = mass;
            this.radius = radius;
            this.momentOfInertia = 0.4 * mass * radius * radius; // 实心球 I = 2/5 mr²
            this.vx = 0;
            this.vy = 0;
            this.omega = 0;
        }

        // 总动能(平动+转动)
        double kineticEnergy() {
            double v2 = vx * vx + vy * vy;
            return 0.5 * mass * v2 + 0.5 * momentOfInertia * omega * omega;
        }

        // 势能(以y=0为参考面)
        double potentialEnergy() {
            return mass * GRAVITY * (y - radius);
        }

        // 总机械能
        double totalEnergy() {
            return kineticEnergy() + potentialEnergy();
        }
    }

    // 地形接口
    interface Terrain {
        double height(double x);
        double slope(double x); // 导数 dh/dx
    }

    // 正弦波地形
    static class SineTerrain implements Terrain {
        double amplitude;
        double frequency;

        SineTerrain(double amplitude, double frequency) {
            this.amplitude = amplitude;
            this.frequency = frequency;
        }

        @Override
        public double height(double x) {
            return amplitude * Math.sin(frequency * x);
        }

        @Override
        public double slope(double x) {
            return amplitude * frequency * Math.cos(frequency * x);
        }
    }

    // 分段线性地形(带平台与斜坡)
    static class PiecewiseTerrain implements Terrain {
        private List<double[]> segments = new ArrayList<>(); // {x_start, x_end, h_start, slope}

        void addSegment(double xStart, double xEnd, double hStart, double slope) {
            segments.add(new double[]{xStart, xEnd, hStart, slope});
        }

        @Override
        public double height(double x) {
            for (double[] seg : segments) {
                if (x >= seg[0] && x <= seg[1]) {
                    return seg[2] + seg[3] * (x - seg[0]);
                }
            }
            return 0;
        }

        @Override
        public double slope(double x) {
            for (double[] seg : segments) {
                if (x >= seg[0] && x <= seg[1]) {
                    return seg[3];
                }
            }
            return 0;
        }
    }

    private Terrain terrain;
    private double frictionCoeff; // 滚动摩擦系数
    private double restitution;   // 碰撞恢复系数

    public RollingBallSimulator(Terrain terrain, double frictionCoeff, double restitution) {
        this.terrain = terrain;
        this.frictionCoeff = frictionCoeff;
        this.restitution = restitution;
    }

    /**
     * 单步物理积分:使用半隐式欧拉法
     */
    public void step(BallState ball) {
        double groundH = terrain.height(ball.x);
        double slope = terrain.slope(ball.x);
        double theta = Math.atan(slope);

        // 判断接触状态
        boolean onGround = ball.y <= groundH + ball.radius + 0.001;

        if (onGround) {
            // 地面接触:法向力约束
            ball.y = groundH + ball.radius;

            // 切向加速度(沿坡面方向)
            // 无滑动条件:a = alpha * r
            // 平动:m*a = mg*sin(theta) - f
            // 转动:I*alpha = f*r
            // 联立得:a = (g*sin(theta)) / (1 + I/(m*r²)) = 5/7 * g * sin(theta)
            double effectiveAcc = (GRAVITY * Math.sin(theta)) /
                (1.0 + ball.momentOfInertia / (ball.mass * ball.radius * ball.radius));

            // 摩擦力方向与运动方向相反(阻碍滑动)
            double tangentVel = ball.vx * Math.cos(theta) + ball.vy * Math.sin(theta);
            double frictionDecel = frictionCoeff * GRAVITY * Math.cos(theta);

            // 沿坡面方向的加速度
            double alongSlope = effectiveAcc;
            if (Math.abs(tangentVel) > 0.01) {
                alongSlope -= Math.signum(tangentVel) * frictionDecel;
            }

            // 分解为x,y方向
            double ax = alongSlope * Math.cos(theta);
            double ay = alongSlope * Math.sin(theta);

            // 半隐式欧拉积分
            ball.vx += ax * DT;
            ball.vy += ay * DT;
            ball.x += ball.vx * DT;
            ball.y += ball.vy * DT;

            // 角速度更新(无滑动条件)
            double vTangent = ball.vx * Math.cos(theta) + ball.vy * Math.sin(theta);
            ball.omega = vTangent / ball.radius;
        } else {
            // 空中:自由落体
            ball.vy -= GRAVITY * DT;
            ball.x += ball.vx * DT;
            ball.y += ball.vy * DT;
        }
    }

    /**
     * 模拟小球从起点滚到终点的过程
     */
    public List<double[]> simulate(BallState ball, double maxTime, double targetX) {
        List<double[]> trajectory = new ArrayList<>();
        trajectory.add(new double[]{ball.x, ball.y, 0});

        double t = 0;
        double initialEnergy = ball.totalEnergy();

        while (t < maxTime && ball.x < targetX + 5) {
            step(ball);
            t += DT;

            if ((int)(t / DT) % 10 == 0) { // 每10步记录一次
                trajectory.add(new double[]{ball.x, ball.y, t});
            }

            // 能量检查(用于调试,允许小幅耗散)
            double currentE = ball.totalEnergy();
            if (Math.abs(currentE - initialEnergy) > initialEnergy * 0.1) {
                // 能量变化过大,可能数值不稳定
            }
        }

        return trajectory;
    }

    // A*搜索节点
    static class AStarNode implements Comparable<AStarNode> {
        int ix, iy; // 网格索引
        double g;   // 从起点代价
        double h;   // 启发式估计
        AStarNode parent;

        AStarNode(int ix, int iy) {
            this.ix = ix;
            this.iy = iy;
        }

        double f() { return g + h; }

        @Override
        public int compareTo(AStarNode o) {
            return Double.compare(this.f(), o.f());
        }

        @Override
        public boolean equals(Object o) {
            if (this == o) return true;
            if (!(o instanceof AStarNode)) return false;
            AStarNode other = (AStarNode) o;
            return ix == other.ix && iy == other.iy;
        }

        @Override
        public int hashCode() {
            return ix * 10000 + iy;
        }
    }

    /**
     * A*路径规划:在离散网格上寻找从起点到终点的最优路径
     * @param startX, startY 起点
     * @param goalX, goalY 终点
     * @param gridSize 网格尺寸
     * @param cellSize 每个格子的实际大小
     */
    public List<double[]> planPath(double startX, double startY,
                                    double goalX, double goalY,
                                    int gridSize, double cellSize) {
        // 初始化网格
        boolean[][] blocked = new boolean[gridSize][gridSize];
        double[][] heightMap = new double[gridSize][gridSize];

        for (int i = 0; i < gridSize; i++) {
            for (int j = 0; j < gridSize; j++) {
                double x = i * cellSize;
                double y = j * cellSize;
                heightMap[i][j] = terrain.height(x);
                // 过陡区域标记为障碍(坡度>60度)
                double s = Math.abs(terrain.slope(x));
                if (s > Math.tan(Math.toRadians(60))) {
                    blocked[i][j] = true;
                }
            }
        }

        int si = (int)(startX / cellSize);
        int sj = (int)(startY / cellSize);
        int gi = (int)(goalX / cellSize);
        int gj = (int)(goalY / cellSize);

        PriorityQueue<AStarNode> open = new PriorityQueue<>();
        Set<AStarNode> closed = new HashSet<>();
        Map<String, AStarNode> openMap = new HashMap<>();

        AStarNode start = new AStarNode(si, sj);
        start.g = 0;
        start.h = heuristic(si, sj, gi, gj, cellSize, heightMap);
        open.offer(start);
        openMap.put(si + "," + sj, start);

        int[][] dirs = {{-1,0},{1,0},{0,-1},{0,1},{-1,-1},{-1,1},{1,-1},{1,1}};

        while (!open.isEmpty()) {
            AStarNode current = open.poll();
            openMap.remove(current.ix + "," + current.iy);

            if (current.ix == gi && current.iy == gj) {
                return reconstructPath(current, cellSize);
            }

            closed.add(current);

            for (int[] d : dirs) {
                int ni = current.ix + d[0];
                int nj = current.iy + d[1];

                if (ni < 0 || ni >= gridSize || nj < 0 || nj >= gridSize) continue;
                if (blocked[ni][nj]) continue;

                AStarNode neighbor = new AStarNode(ni, nj);
                if (closed.contains(neighbor)) continue;

                // 移动代价:距离 + 上坡惩罚(能量代价)
                double dx = d[0] * cellSize;
                double dy = d[1] * cellSize;
                double dist = Math.hypot(dx, dy);
                double heightDiff = heightMap[ni][nj] - heightMap[current.ix][current.iy];
                double moveCost = dist + Math.max(0, heightDiff * 2); // 上坡更 costly

                double tentativeG = current.g + moveCost;

                String key = ni + "," + nj;
                AStarNode existing = openMap.get(key);
                if (existing == null || tentativeG < existing.g) {
                    if (existing == null) {
                        existing = neighbor;
                        openMap.put(key, existing);
                        open.offer(existing);
                    }
                    existing.g = tentativeG;
                    existing.h = heuristic(ni, nj, gi, gj, cellSize, heightMap);
                    existing.parent = current;
                }
            }
        }

        return new ArrayList<>(); // 无路径
    }

    private double heuristic(int i, int j, int gi, int gj, double cellSize, double[][] heightMap) {
        double dx = (gi - i) * cellSize;
        double dy = (gj - j) * cellSize;
        double dist = Math.hypot(dx, dy);
        // 加入高度差的启发估计
        double hDiff = Math.abs(heightMap[gi][gj] - heightMap[i][j]);
        return dist + hDiff;
    }

    private List<double[]> reconstructPath(AStarNode node, double cellSize) {
        List<double[]> path = new ArrayList<>();
        while (node != null) {
            path.add(new double[]{node.ix * cellSize, node.iy * cellSize});
            node = node.parent;
        }
        Collections.reverse(path);
        return path;
    }

    public static void main(String[] args) {
        // 创建起伏地形
        PiecewiseTerrain terrain = new PiecewiseTerrain();
        terrain.addSegment(0, 3, 0, 0.1);      // 缓坡
        terrain.addSegment(3, 5, 0.3, -0.3);   // 下坡
        terrain.addSegment(5, 7, -0.3, 0);     // 平台
        terrain.addSegment(7, 10, -0.3, 0.2);  // 上坡到终点

        RollingBallSimulator sim = new RollingBallSimulator(terrain, 0.05, 0.7);

        // 模拟小球滚动
        BallState ball = new BallState(0, terrain.height(0) + 0.1, 0.5, 0.1);
        ball.vx = 1.0; // 初始推动

        System.out.println("=== 滚球物理模拟 ===");
        double e0 = ball.totalEnergy();
        List<double[]> traj = sim.simulate(ball, 10, 10);

        System.out.printf("初始总能量: %.3f J\n", e0);
        System.out.printf("最终总能量: %.3f J\n", ball.totalEnergy());
        System.out.printf("最终位置: (%.2f, %.2f)\n", ball.x, ball.y);
        System.out.printf("最终速度: %.2f m/s\n", Math.hypot(ball.vx, ball.vy));
        System.out.println("轨迹点数量: " + traj.size());

        // A*路径规划
        System.out.println("\n=== A*路径规划 ===");
        List<double[]> path = sim.planPath(0, 0.5, 9, 1.5, 50, 0.25);
        System.out.println("规划路径点数: " + path.size());
        for (int i = 0; i < Math.min(path.size(), 10); i++) {
            double[] p = path.get(i);
            System.out.printf("  路径点%d: (%.2f, %.2f)\n", i, p[0], p[1]);
        }
    }
}

算法复杂度分析

步骤 时间复杂度 空间复杂度 说明
单步物理积分 O(1) O(1) 半隐式欧拉法
轨迹模拟 O(T/DT) O(T/DT) T为总时间,DT为步长
A*搜索 O(N log N) O(N) N为网格节点数
路径重建 O(L) O(L) L为路径长度

扩展方向

  1. 三维扩展:将地形扩展为 h(x,z),小球在三维空间滚动,用四元数表示姿态
  2. 关节体连接:多个小球用杆连接,形成蛇形滚动物体,用约束动力学求解
  3. 最优控制:用庞特里亚金极大值原理求解最少时间到达终点的控制策略
  4. 神经网络控制器:训练深度网络直接输出施加的力/扭矩,用进化策略优化

总结

本文从刚体动力学出发,推导了无滑动滚动条件下小球的等效加速度公式,通过半隐式欧拉积分实现稳定高效的物理模拟。同时利用A*搜索在考虑地形坡度和能量代价的图中规划最优路径,并用能量守恒验证数值积分的物理正确性。该框架可扩展至更复杂的三维滚球游戏和机器人运动规划场景。