%% 恒力始终垂直于速度的运动轨迹模拟
clear; clc; close all;
hold on;
grid on;
axis equal;
%% 参数设置
m = 2.0; % 质量 (kg)
v0 = 5.0; % 初始速度大小 (m/s)
f = 2.0; % 恒力大小 (N)
theta0 = 0; % 初始速度方向 (弧度,0表示+x方向)
t_total = 100; % 总时间 (s)
dt = 0.0001; % 时间步长 (s)% 计算理论值
r_theory = m * v0^2 / f; % 理论半径
omega_theory = f / (m * v0); % 理论角速度
T_theory = 2*pi / omega_theory; % 理论周期fprintf('理论计算值:\n');
fprintf(' 半径 r = %.3f m\n', r_theory);
fprintf(' 角速度 ω = %.3f rad/s\n', omega_theory);
fprintf(' 周期 T = %.3f s\n', T_theory);%% 数值模拟(欧拉法)
% 初始化数组
n_steps = round(t_total / dt);
t = linspace(0, t_total, n_steps);
x = zeros(1, n_steps);
y = zeros(1, n_steps);
vx = zeros(1, n_steps);
vy = zeros(1, n_steps);% 初始条件
x(1) = 0; % 初始位置
y(1) = 0;
vx(1) = v0 * cos(theta0); % 初始速度
vy(1) = v0 * sin(theta0);% 数值积分
for i = 1:n_steps-1% 当前速度大小和方向v = sqrt(vx(i)^2 + vy(i)^2);% 计算垂直于速度方向的单位向量% 垂直于(vx, vy)的向量可以是(-vy, vx)或(vy, -vx)% 这里选择使质点顺时针旋转的方向nx = -vy(i) / v;ny = vx(i) / v;% 计算加速度(力垂直于速度)ax = (f / m) * nx;ay = (f / m) * ny;% 更新速度和位置vx(i+1) = vx(i) + ax * dt;vy(i+1) = vy(i) + ay * dt;x(i+1) = x(i) + vx(i) * dt;y(i+1) = y(i) + vy(i) * dt;
endplot(x, y, 'b-', 'LineWidth', 1.5);