原文：https://blog.csdn.net/WhoisPo/article/details/47149395

关于最速曲线的介绍有

http://zhidao.baidu.com/s/daily/2014-04-21/1403015178.html

内容比较丰富，还比较好玩

最速曲线公式

理论解很久之前就已经有了，如下

我找了半天也没有找到这个理论解是如何求出来的方法，但是我找到了一篇怎样用数值方法求最速曲线的算法，这篇文章的题目是《应用斯涅尔公式求解最速下降曲线问题研究》，在百度文库中有。这个方法我觉得很有意思，就在matlab中实现了这个算法。

特别要注意的是文章中提到的，如果入射角大于90°时，要考虑全反射，也就是入射角从正变为负，下降变为上升。

# -*- coding: utf-8 -*-
# 下面用数值方法求解最速曲线

#导入必要的模块
import numpy as np
import matplotlib.pyplot as plt


import math

x0 = 10;
y0 = 10;
etof = 1e-2; # 误差，精度为0
.1
a0 = 0; # 入射角的边值
a1 = math.pi / 30;

dex = 0.1; # 最大步长
dey = 0.1;

data_x=[]
data_y=[]
while (1):
    a_init = (a0 + a1) / 2; # 入射角的初值

    a = a_init;
    dx = dex; # 起始步长
    dy = 1/math.tan(a) * dx;
    while (dy > dey):
        dx = dx / 10; # 细化步长
        dy = 1/math.tan(a) * dx;
    #end
    x = dx;
    y = dy;
    data_x.append(x)
    data_y.append(y)
    while (x + dx <= x0):
        yt = y +1/math.tan(a) * dx;
        if (yt < 0):
            break; # 已久上升到初始高度，结束
        #end
        if (math.sin(a) / math.sqrt(y) * math.sqrt(yt) >= 1): # 考虑全反射
            a = -a;
            yt = y + 1/math.tan(a) * dx;
        #end
        at = math.asin(math.sin(a) / math.sqrt(y) * math.sqrt(yt));
        x = x + dx;
        y = yt;
        # data_x.append(x) 这个不要
        # data_y.append(y)
        a = at;

        dy = 1/math.tan(a) * dx;
        while (dy < dey / 10 and dx < dex): # 适当放大步长
            dx = dx * 10;
            dy = 1/math.tan(a) * dx;
        #end
    #end

    if (abs(y - y0) < etof and x + dx > x0):
        break;
    #end
    if (x + dx <= x0 or y < y0):
        a1 = a_init;
    else:
        a0 = a_init;
    #end
#end

a = (a0 + a1) / 2;

# scatter(0, 0, 'go');
dx = dex;
dy = 1/math.tan(a) * dx;
while (dy > dey):
    dx = dx / 10;
    dy = 1/math.tan(a) * dx;
#end
# scatter(dx, dy, 'go');
x = dx;
y = dy;
data_x.append(x)
data_y.append(y)
while (x + dx <= x0):
    yt = y + 1/math.tan(a) * dx;
    if (math.sin(a) / math.sqrt(y) * math.sqrt(yt) >= 1):
        a = -a;
        yt = y + 1/math.tan(a) * dx;
    #end
    at = math.asin(math.sin(a) / math.sqrt(y) * math.sqrt(yt));
    x = x + dx;
    y = yt;
    data_x.append(x)
    data_y.append(y)
    a = at;
    # scatter(x, y, 'go');
    dy = 1/math.tan(a) * dx;
    while (dy < dey / 10 and dx < dex):
        dx = dx * 10;
        dy = 1/math.tan(a) * dx;

#产生测试数据

fig = plt.figure()
ax1 = fig.add_subplot(111)
#设置标题
ax1.set_title('Scatter Plot')
#设置X轴标签
plt.xlabel('X')
#设置Y轴标签
plt.ylabel('Y')
data_y = [20-c for c in data_y]
#画散点图
sValue = data_x
ax1.scatter(data_x,data_y,s=sValue,c='r',marker='x')
#设置图标
plt.legend('x1')
#显示所画的图
plt.show()