混沌台球模拟

Question

我来寻求一些数学和编程方面的帮助。

，我想做什么？，我试图实现一个混沌台球系统的模拟，遵循算法在这个节选中。

如何尝试?使用numpy和matplotlib的，实现了以下代码

def boundaryFunction(parameter):
return 1 + 0.1 * np.cos(parameter)

def boundaryDerivative(parameter):
    return -0.1 * np.sin(parameter)

def trajectoryFunction(parameter):
    aux = np.sin(beta - phi) / np.sin(beta - parameter)
    return boundaryFunction(phi) * aux

def difference(parameter):
    return trajectoryFunction(parameter) - boundaryFunction(parameter)

def integrand(parameter):
    rr = boundaryFunction(parameter)
    dd = boundaryDerivative (parameter)
    return np.sqrt(rr ** 2 + dd ** 2)

##### Main #####

length_vals = np.array([], dtype=np.float64)
alpha_vals = np.array([], dtype=np.float64)

# nof initial phi angles, alpha angles, and nof collisions for each.
n_phi, n_alpha, n_cols, count = 10, 10, 10, 0
# Length of the boundary
total_length, err = integrate.quad(integrand, 0, 2 * np.pi)

for phi in np.linspace(0, 2 * np.pi, n_phi):
    for alpha in np.linspace(0, 2 * np.pi, n_alpha):
        for n in np.arange(1, n_cols):

            nu = np.arctan(boundaryFunction(phi) / boundaryDerivative(phi))
            beta = np.pi + phi + alpha - nu

            # Determines next impact coordinate.
            bnds = (0, 2 * np.pi)
            phi_new = optimize.minimize_scalar(difference, bounds=bnds, method='bounded').x

            nu_new =  np.arctan(boundaryFunction(phi_new) / boundaryDerivative(phi_new))
            # Reflection angle with relation to tangent.
            alpha_new = phi_new - phi + nu - nu_new - alpha
            # Arc length for current phi value.
            arc_length, err = integrate.quad(integrand, 0, phi_new)

            # Append values to list
            length_vals = np.append(length_vals, arc_length / total_length)
            alpha_vals = np.append(alpha_vals, alpha)


        count += 1
    print  "{}%" .format(100 * count / (n_phi * n_alpha))

是什么问题？在计算phi_new时，方程有两个解(假设边界是凸的，即。)我必须强制phi_new是与phi不同的解决方案，但我不知道如何做到这一点。代码还有更多的问题吗？

输出应该是什么？是Sα，看上去像这样的相位图。

任何帮助都是非常感谢的！提前谢谢。

Answer 1

你可以尝试的一种方法是(考虑到只有两种解决方案)

epsilon = 1e-7 # tune this
delta = 1e-4 # tune this
# ...
bnds = (0, 2 * np.pi)
phi_new = optimize.minimize_scalar(difference, bounds=bnds, method='bounded').x
if abs(phi_new - phi) < epsilon:
    bnds_1 = (0, phi - delta)
    phi_new_1 = optimize.minimize_scalar(difference, bounds=bnds_1, method='bounded').x
    bnds_2 = (phi + delta, 2 * np.pi)
    phi_new_2 = optimize.minimize_scalar(difference, bounds=bnds_2, method='bounded').x
    if difference(phi_new_1) < difference(phi_new_2):
        phi_new = phi_new_1
    else:
        phi_new = phi_new_2

或者，您也可以引入一个惩罚术语，例如，带有epsilon和增量的适当选择的delta*exp(eps/(x-phi)^2)。