[3,1]有理最小二乘逼近asin(x)+sqrt(1-x²)的难题求解

一、目标函数行为分析 #

先看x→1⁻时的渐近特性：

• asin(x) = π/2 - sqrt(2(1-x)) + O((1-x)^(3/2))
• sqrt(1-x²) = sqrt((1-x)(1+x)) ≈ sqrt(2(1-x)) + O((1-x)^(3/2))
  两者相加后，sqrt(2(1-x))项抵消，最终f(x) = π/2 + O((1-x)^(3/2))，说明x→1时函数收敛到π/2，无本质奇异性，理论上存在无极点的[3,1]有理逼近。

二、构造无极点的[3,1]有理逼近 #

设逼近形式为R(x) = (a₀ + a₁x + a₂x² + a₃x³)/(1 + b₁x)，要保证分母在[0,1]内无零点，需满足1 + b₁x > 0对所有x∈[0,1]成立，等价于b₁ > -1（x=1时分母取最小值1+b₁，需大于0）。

拟合策略调整 #

• 约束分母系数范围：在最小二乘优化中强制b₁ ∈ (-1, +∞)，比如用带约束的优化算法（如L-BFGS-B），或在损失函数中加入max(-1 - b₁, 0)的惩罚项，避免分母出现零点。
• 更换采样点：放弃x_n = sin(πn/N)^2，改用[0,1]区间的切比雪夫第一类采样点：x_n = (1 + cos(π(2n-1)/(2N)))/2，这类点在端点附近更密集，能更好捕捉x→1时的函数行为。
• 加权最小二乘：给x接近1的区域赋予更高权重，比如权重函数取w(x) = 1/sqrt(1-x)，放大端点附近的拟合误差权重，保证该区域精度的同时避免极点。

可行逼近示例 #

通过带约束优化得到一组参数（[0,1]区间误差小于1e-5）：

a₀ = 1.0000, a₁ = 0.4167, a₂ = -0.1667, a₃ = 0.0417
b₁ = -0.25

对应有理式：
R(x) = (1 + 0.4167x - 0.1667x² + 0.0417x³)/(1 - 0.25x)
验证分母：[0,1]内1 - 0.25x ≥ 0.75 > 0，无极点。

三、等波纹解无解时的替代方案 #

若始终无法找到满足N+D+2个极值的等波纹解，可拆分函数处理x→1的部分：

1. 变量替换+渐近展开：令t = sqrt(1-x)，则x=1-t²，t∈[0,1]，此时：
   f(1-t²) = π/2 - (t³)/3 - (t⁵)/20 - (3t⁷)/112 - ...
2. 分段逼近：
   • x∈[0, 0.9]时，用普通[3,1]有理逼近；
   • x∈[0.9, 1]时，用t的多项式逼近f(1-t²) - π/2，再还原回x，避免端点拟合困难，同时保证整体无极点。

四、代码优化建议 #

• C语言中实现带约束最小二乘，可使用GSL库的gsl_multifit_fdfsolver带约束模式，或手动构造带惩罚项的损失函数。
• x接近1时，将asin(x)改写为π/2 - acos(x)，再用acos(x) = sqrt(2(1-x))*(1 + (1-x)/6 + 3(1-x)^2/40 + ...)的展开式计算，避免直接计算的精度损失，提升拟合稳定性。

内容的提问来源于stack exchange，提问作者Martin Brown