数值求解程序如下：
clear;clc;
% 参数设置
L = 10000; % 角动量
m = 4; % 质量
g = 9.81; % 重力加速度
alpha = pi/4; % 假设圆锥半顶角为45度
dt = 0.0001; % 时间步长
t_end = 100; % 结束时间
r0 = (L^2*tan(alpha)/(m^2*g))^(1/3); % 初始轨道半径由角动量决定
a = 100;
dr0 = L/(m*r0)/a; % v(0)要远小于角向初始速度L/(m*r0)
n = 3;
% 初始化变量
t = 0;
t_values = 0:dt:t_end;
r_values = zeros(1, length(t_values));%r(t)
dr_values = zeros(1, length(t_values));%v(t)
% 迭代求解
for i = 1:length(t_values)
if i == 1
% 第一个时间步长，使用初始条件
r_values(i) = r0;
dr_values(i) = dr0;
else
% 使用欧拉法更新r和dr
r_values(i) = r_values(i-1) + dt * dr_values(i-1);
dr_values(i) = dr_values(i-1) + dt * csc(alpha)^2*((L^2) / (m^2 * r_values(i-1)^n) - g * cot(alpha));
end
t = t + dt;
end
% 绘图
figure(1)
plot(t_values, r_values-r0, 'r-', 'LineWidth', 2);
hold on;
omega = sqrt(6*g/(r0*sin(2*alpha)));%微振动角频率
plot(t_values, max(r_values-r0)*sin(t_values*omega), 'b--', 'LineWidth', 2);
xlabel('Time (t)');
ylabel('r(t)');
legend('真实解','谐振子近似解');
title(['径向运动方程解，L = ',num2str(L),'; v0 = L/(m*r0)/',num2str(a)]);
figure(2)
plot(t_values, L./(m*r_values.^2), 'r-', 'LineWidth', 2);
xlabel('Time (t)');
ylabel('θdot(t)');
title(['角速度，L = ',num2str(L),'; v0 = L/(m*r0)/',num2str(a)]);
[~, locs] = findpeaks(r_values);
disp(['平衡位置R0为：', num2str(r0)]);
disp('实际周期为：');
disp((locs(2)-locs(1))*dt);
disp((locs(ceil(end/2))-locs(ceil(end/2)-1))*dt);
disp((locs(end)-locs(end-1))*dt);
disp('微振动理论值：');
disp(2*pi/omega);

另外，如果小球换成微观粒子，这个系统可以量子化吗？直接把哈密顿量里的物理量都换成对应的力学量算符就行了吗？能找到基态能级和基态波函数吗？

当几何约束不是圆锥，而是任意单增函数f(x)绕y轴旋转得到的柱面，运动方程会变成什么样？

python程序

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from scipy.signal import find_peaks
from math import exp, sqrt, sin, cos, pi
#曲面形状
R = 10
omega = 0.5
#f = lambda r: R*(1-cos(omega*r))
#df = lambda r: R*omega*sin(omega*r)
#ddf = lambda r: R*(omega**2)*cos(omega*r)
f = lambda r: R*(1-exp(-(r/R)**2))
df = lambda r: 2*(r/R)*exp(-(r/R)**2)
ddf = lambda r: (2/R)*(1-2*(r/R)**2)*exp(-(r/R)**2)
#f = lambda r: R-sqrt(R**2-r**2) # 定义符号表达式
#df = lambda r: r/sqrt(R**2 - r**2) # 一阶导数
#ddf = lambda r: 1/sqrt(R**2 - r**2) + r**2/(R**2 - r**2)**(3/2)
# 二阶导数
# 参数设置
r0 = 1 # 初始轨道半径
g = 9.81 # 重力加速度
m = 10 # 质量
a = 7 #微扰量
# 计算初始角动量 L
L = sqrt(m**2 * r0**3 * g * df(r0))
# 定义微分方程函数
def ode45_fun(t, y, m, g, L, df, ddf):
r, dr = y
dydt = [dr, (df(r) * ddf(r) * dr**2 + (L**2) / (m**2 * r**3) - g * df(r)) / (1 + df(r)**2)]
return dydt
# 初始条件
y0 = [r0, a * L / (m * r0)] # [r, dr]
# 使用 solve_ivp 求解
t_span = (0, 50)
t_eval = np.linspace(t_span[0], t_span[1], 10000)
sol = solve_ivp(ode45_fun, t_span, y0, args=(m, g, L, df, ddf), t_eval=t_eval)
peaks, _ = find_peaks(sol.y[0])
# 如果找到峰值，计算相邻峰值之间的时间差
if len(peaks) > 1:
# 计算相邻峰值之间的时间差
periods = np.diff(sol.t[peaks])
average_period = np.mean(periods) if len(periods) > 1 else periods[0]
print(f"找到的周期: {average_period}")
else:
print("没有找到足够的峰值来计算周期。")
# 绘图
plt.figure(1)
plt.plot(sol.t, sol.y[0] - r0, 'r-', linewidth=2)
plt.xlabel('Time (t)')
plt.ylabel('r(t) - r(0)')
plt.title(f'radical solution, L = {L}, vθ/vr = {a}')
plt.show()
vec_squared = [f(x) for x in sol.y[0]]
plt.figure(2)
plt.plot(sol.t, vec_squared , 'r-', linewidth=2)
plt.xlabel('Time (t)')
plt.ylabel('r(t) - r(0)')
plt.title(f'vertical solution, L = {L}, vθ/vr = {a}')
plt.show()
plt.figure(3)
plt.plot(sol.t, sol.y[1], 'r-', linewidth=2)
plt.xlabel('Time (t)')
plt.ylabel('r_dot(t)')
plt.title(f'radical velocity, L = {L}, vθ/vr = {a}')
plt.show()