令$(x_k,y_k)$为第k次的反射点，$(a_k,b_k)$是第k次反射时入射线的方向向量。任意一点处的导数，亦即切线斜率可用隐函数求导法得
$$8x+2yy'=0,$$
即
$$y'=-\frac{4x}{y},$$
从而在反射点出切线的方向向量为$(y_k,-4x_k)$，法线的方向向量为$(4x_k,y_k)$。因此第k次反射时反射线的方向向量可以用正交投影计算：
\begin{align*}
\begin{pmatrix}a_{k+1} \\ b_{k+1}\end{pmatrix}&=\begin{pmatrix}a_{k} \\ b_{k}\end{pmatrix}-\frac{2}{16x_k^2+y_k^2}\langle\begin{pmatrix}a_{k} \\ b_{k}\end{pmatrix}, \begin{pmatrix}4x_k \\ y_k\end{pmatrix}\rangle\begin{pmatrix}4x_k \\ y_k\end{pmatrix} \\
&=\begin{pmatrix}a_{k} \\ b_{k}\end{pmatrix}-\frac{8a_kx_k+2b_ky_k}{16x_k^2+y_k^2}\begin{pmatrix}4x_k \\ y_k\end{pmatrix} \\
&=\frac{1}{16x_k^2+y_k^2}\begin{pmatrix}16a_{k}x_k^2+a_{k}y_k^2-32a_kx_k^2-8b_kx_ky_k \\ 16b_{k}x_k^2+b_{k}y_k^2-8a_kx_ky_k-2b_ky_k^2\end{pmatrix} \\
&=\frac{1}{16x_k^2+y_k^2}\begin{pmatrix}-16a_{k}x_k^2+a_{k}y_k^2-8b_kx_ky_k \\ 16b_{k}x_k^2-b_{k}y_k^2-8a_kx_ky_k\end{pmatrix}
\end{align*}
其中$\frac{1}{16x_k^2+y_k^2}\langle\begin{pmatrix}a_{k} \\ b_{k}\end{pmatrix}, \begin{pmatrix}4x_k \\ y_k\end{pmatrix}\rangle\begin{pmatrix}4x_k \\ y_k\end{pmatrix}$是入射方向$(a_k,b_k)$在法线方向$(4x_k,y_k)$上的正交投影，将其延长至2倍，再用原向量$(a_k,b_k)$减去之，就是反射后的向量，这是镜面反射的一般公式。
由此可得反射后的直线方程为
$$y-y_k=\frac{b_{k+1}}{a_{k+1}}(x-x_k),$$
将其代入椭圆方程$4x^2+y^2=100$得
$$4x^2+(\frac{b_{k+1}}{a_{k+1}}(x-x_k)+y_k)^2=100,$$
整理后可得
$$4a_{k+1}^2x^2+(b_{k+1}x-b_{k+1}x_k+a_{k+1}y_k)^2=100a_{k+1}^2,$$
即
$$(4a_{k+1}^2+b_{k+1}^2)x^2+2b_{k+1}(a_{k+1}y_k-b_{k+1}x_k)x+(a_{k+1}y_k-b_{k+1}x_k)^2-100a_{k+1}^2=0,$$
这个关于x的二次多项式的判别式为
\begin{align*}
\Delta&=4b_{k+1}^2(a_{k+1}y_k-b_{k+1}x_k)^2-4(4a_{k+1}^2+b_{k+1}^2)((a_{k+1}y_k-b_{k+1}x_k)^2-100a_{k+1}^2)\\
&=400a_{k+1}^2(4a_{k+1}^2+b_{k+1}^2)-16a_{k+1}^2(a_{k+1}y_k-b_{k+1}x_k)^2
\end{align*}
从而解得
$$x=\frac{-2b_{k+1}(a_{k+1}y_k-b_{k+1}x_k)\pm\sqrt{\Delta}}{2(4a_{k+1}^2+b_{k+1}^2)}.$$
这两个根分别就是$x_k$和$x_{k+1}$，将$x_{k+1}$代入反射线的方程可以求得$y_{k+1}$。
由题意，初值为$x_1=1.4$, $y_1=-9.6$, $a_1=1.4$, $b_1=-19.7$，重复上述过程直至$|x_{k+1}|\le 0.01$，此时激光可以离开“白腔”。
反射次数很多，舍入误差经多次累计后会影响精度，因此在计算中我们通过缩放以确保$a_k$和$b_k$的绝对值都较小（以下选择的是令其绝对值不超过$10$）。
结果是$354$。
注：本题涉及的导数知识是我国本科《高等数学》课程中的标准内容。正交投影和镜面反射公式在我国本科《线性代数》课程或有涉及，但一般在数学系开设的《解析几何》、《高等代数》、《泛函分析》等课程中都会讲授。
import math
x,y,a,b,k=1.4,-9.6,1.4,-19.7,0
while abs(x)>0.01 or y < 0:
a,b=-16*a*x*x+a*y*y-8*b*x*y,16*b*x*x-b*y*y-8*a*x*y
while abs(a)>10 or abs(b)>10:
a,b=a/10,b/10
s=b/a
xx=(4*x-s*s*x+2*s*y)/(-4-s*s)
yy=s*(xx-x)+y
x,y=xx,yy
k+=1
print k,x,y
|
