采用数值方法计算最速曲线

来源：互联网发布：java log4j日志 xml 编辑：程序博客网时间：2024/06/02 05:29

关于最速曲线的介绍有

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

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

最速曲线公式

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

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

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

Matlab代码如下

<span style="font-size:18px;">x0 = 10;y0 = 10;syms x;f = (1-cos(x)) / (x - sin(x)) - y0 / x0;x = vpasolve(f);       %数值求解器K = x0 / (x - sin(x));s = 0:0.01:x;x = K*(s - sin(s));y = K*(1 - cos(s));plot(x,y, 'r-');view([0,0,-1]);hold on;%-----------------------------------------上面是理论的结果-----------------------------------% 下面用数值方法求解最速曲线etof = 1e-2;         %误差，精度为0.1a0 = 0;        % 入射角的边值a1 = pi/30;    %dex = 0.1;        %最大步长dey = 0.1;while 1a_init = (a0 + a1) / 2;   % 入射角的初值a = a_init;dx = dex;   %起始步长dy = cot(a) * dx;while (dy > dey)dx = dx / 10;         %细化步长dy = cot(a) * dx;endx = dx;y = dy;while (x + dx <= x0)yt = y + cot(a) * dx;if (yt < 0) break;         %已久上升到初始高度，结束endif (sin(a) / sqrt(y) * sqrt(yt) >= 1)            % 考虑全反射a = -a; yt = y + cot(a) * dx;endat = asin(sin(a) / sqrt(y) * sqrt(yt));x = x + dx;y = yt;a = at;dy = cot(a) * dx;while (dy < dey / 10 && dx  < dex)             %适当放大步长dx = dx * 10;dy = cot(a) * dx;endendif (abs(y - y0) < etof && x + dx > x0)break;endif (x + dx <= x0 || y < y0)a1 = a_init;elsea0 = a_init;endenda = (a0 + a1) / 2;scatter(0, 0, 'go');dx = dex;dy = cot(a) * dx;while (dy > dey)dx = dx / 10;dy = cot(a) * dx;endscatter(dx, dy, 'go');x = dx;y = dy;while (x + dx<= x0)yt = y + cot(a) * dx;if (sin(a) / sqrt(y) * sqrt(yt) >= 1)a = -a;yt = y + cot(a) * dx;endat = asin(sin(a) / sqrt(y) * sqrt(yt));x = x + dx;y = yt;a = at;scatter(x, y, 'go');dy = cot(a) * dx;while (dy < dey / 10 && dx < dex)dx = dx * 10;dy = cot(a) * dx;endend</span>

效果如下

0 0