计算方法实验1：圆形镜面成像问题

本文主要是介绍计算方法实验1：圆形镜面成像问题，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

二维平面上的镜面反射示意图
在这里插入图片描述

Algorithm Description

设 $T(cos\theta,sin\theta)$ ,则有
$PT+QT=\sqrt{(P_x-cos\theta)^2+sin^2\theta}+\sqrt{(Q_x-cos\theta)^2+(Q_y-sin\theta)^2}$
$PT+QT=\sqrt{P_x^2-2P_xcos\theta+1}+\sqrt{Q_x^2+Q_y^2+1-2Q_xcos\theta-2Q_ysin\theta}$
由费马原理，光线沿 $PT + QT$ 最短的路径传播，因此只需对上式求导求极小值点。关于 $\theta$ 求导得
$\frac{P_xsin\theta}{\sqrt{P_x^2-2P_xcos\theta+1}}+\frac{Q_xsin\theta-Q_ycos\theta}{\sqrt{Q_x^2+Q_y^2+1-2Q_xcos\theta-2Q_ysin\theta}}$
故只需用二分法解非线性方程
$P_xsin\theta\sqrt{Q_x^2+Q_y^2+1-2Q_xcos\theta-2Q_ysin\theta}+(Q_xsin\theta-Q_ycos\theta)\sqrt{P_x^2-2P_xcos\theta+1}=0$
$T_x=cos\theta$
$T_y=sin\theta$
由对称性易知
$Rx=\frac{2Q_y-Qxtan\theta -kQ_x +k(2Tx - Qx) - 2Ty}{k-tan\theta}$
$Rx)\theta$

Code

#include<iostream>
#include<cmath>
#include <iomanip>
using namespace std;int main(int argc, char *argv[])
{if(argc != 4) {cout << "Usage: " << argv[0] << " <Px> <Qx> <Qy>\n";return 1;}long double Px = stold(argv[1]);long double Qx = stold(argv[2]);long double Qy = stold(argv[3]);if(Px >= -1 || Qx >= 0 || Qy <= 0 || Qx * Qx + Qy * Qy <= 1) {cout << "输入错误，Px应该小于-1，Qx应该小于0，Qy应该大于0，Qx^2+Qy^2应该大于1\n";return 1;}long double Tx = 0;long double Ty = 0;long double theta = 0;long double low = 3.14159265358979323 , high = 1.570796326794896;long double k = 0;while(1){theta = (low + high) / 2;long double eq1 = Px * sin(theta);long double eq2 = sqrt(Qx * Qx + Qy * Qy +1 - 2 * Qx * cos(theta) - 2 * Qy * sin(theta));long double eq3 = (Qx * sin(theta) - Qy * cos(theta));long double eq4 = sqrt(Px * Px -2 * Px * cos(theta) + 1);long double res = eq1 * eq2 + eq3 * eq4;long double absres = abs(res);if(absres <= 1e-7){Tx = cos(theta);Ty = sin(theta);k = 1 / tan(3.14159265358979323 - theta);break;}else if(res > 0){low = theta;}else{high = theta;}}long double eq5 = 2 * Qy - Qx * tan(theta) + k * (2 * Tx - Qx) - 2 * Ty;long double eq6 = k - tan(theta); long double Rx = eq5 / eq6;long double Ry = Qy - tan(theta) * (Qx - Rx);cout << "T = (" << fixed << setprecision(6) << Tx << " , " << fixed << setprecision(6) << Ty << ") , R = (" << fixed << setprecision(6) << Rx << " , " << fixed << setprecision(6) << Ry << ")";return 0;
}

Results

$P = (- 2, 0), Q = (- 1, 1) : T = (- 0.885670, 0.464316), R = (- 0.380057, 0.674993)$

$P = (- 10, 0), Q = (- 2, 1) : T = (- 0.959312, 0.282350), R = (0.304214, 0.321811)$

$P = (- 1.000001, 0), Q = (- 2, 2) : T = (- 1.000000, 0.000002), R = (0.000007, 1.999996)$

$P = (- 2, 0), Q = (- 1, 0.000001) : T = (- 1.000000, 0.000001), R = (- 1.000000, 0.000001)$

$P = (- 2.33, 0), Q = (- 3, 1) : T = (- 0.989279, 0.146038), R = (1.182424, 0.382590)$

$P = (- 3, 0), Q = (- 1, 0.5) : T = (- 0.922615, 0.385721), R = (- 0.786920, 0.410917)$

$P = (- 3, 0), Q = (- 2, 10) : T = (- 0.827028, 0.562160), R = (8.380296, 2.944148)$

$P = (- 3, 0), Q = (- 3, 1) : T = (- 0.987408, 0.158192), R = (1.187435, 0.329136)$

$P = (- 10, 0), Q = (- 2, 1) : T = (- 0.959312, 0.282350), R = (0.304214, 0.321811)$

$P = (- 1024, 0), Q = (- 8, 4) : T = (- 0.970066, 0.242842), R = (7.000894, 0.244735)$

Conclusion

本实验提高精度的主要措施：

使用long double 类型；
用较多位数来表示 $\pi$ ；
将含有除法的方程交叉相乘，通过乘法代替除法，以减少在求商时的误差；
在用二分法解非线性方程时，限制条件是结果的绝对值 $\leq10^{-7}$ 。

但由于本实验的方程较为复杂，含有三角函数、根式、平方等易造成误差放大的因素，暂时还未找到其他较好的减少误差的办法，但还尝试了另一种思路，叙述如下：

设OT的延长线交PQ于R，则由角平分线定理，有 $\frac{QT}{PT}=\frac{QR}{RP}=\frac{Q_y}{R_y}-1$
设 $T(x,y),k=\frac{Qy}{Qx-Px}$ ，带入上述方程化简得：
$\frac{Q_y^2(kx-y)^2+k^2P_x^2y^2-2kQ_yPxy(kx-y)}{k^2Px^2y^2}=\frac{Q_y^2+Q_x^2+1-2yQ_y-2xQ_x}{1+P_x^2-2xP_x}$
故只需用for循环遍历或二分法解非线性方程 $Q_y^2(kx-y)^2+k^2P_x^2y^2-2kQ_yPxy(kx-y))(1+P_x^2-2xP_x)-(k^2Px^2y^2)(Q_y^2+Q_x^2+1-2yQ_y-2xQ_x)=0$
$y=\sqrt{1-x^2}$

由对称性易知
$Rx=\frac{\frac{QxTy}{Tx} + \frac{Tx(2Tx - Qx)}{Ty} + 2Ty -2Qy}{\frac{Tx}{Ty}+\frac{Ty}{Tx}}$
$Ry=Qy-\frac{Ty}{Tx(Qx - Rx)}$
但验证结果时发现，上述方法在遇到一些极端情况如 $P = (- 1.000001, 0), Q = (- 2, 2)$ 时，中间的运算过程会出现极小的浮点数导致无法继续运算，所以这种思路在细节上仍有待改进。