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

关于最速曲线的介绍有

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.1
a0 = 0;        % 入射角的边值
a1 = pi/30;    %

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

while 1
	a_init = (a0 + a1) / 2;   % 入射角的初值
	
	a = a_init;
	dx = dex;   %起始步长
	dy = cot(a) * dx;
	while (dy > dey)
		dx = dx / 10;         %细化步长
		dy = cot(a) * dx;
	end
	x = dx;
	y = dy;
	while (x + dx <= x0)
		yt = y + cot(a) * dx;
		if (yt < 0) 
			break;         %已久上升到初始高度，结束
		end
		if (sin(a) / sqrt(y) * sqrt(yt) >= 1)            % 考虑全反射
			a = -a; 
			yt = y + cot(a) * dx;
		end
		at = 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;
		end
	end

	if (abs(y - y0) < etof && x + dx > x0)
		break;
	end
	if (x + dx <= x0 || y < y0)
		a1 = a_init;
	else
		a0 = a_init;
	end	
end


a = (a0 + a1) / 2;

scatter(0, 0, 'go');
dx = dex;
dy = cot(a) * dx;
while (dy > dey)
	dx = dx / 10;
	dy = cot(a) * dx;
end
scatter(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;
	end
	at = 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;
	end
end</span>

效果如下