两参数曲面的相交曲线计算问题

lihpb00

这种有两条曲线相交而成的曲线应该不存在标准的解析式吧，可能只有对应的数值解

Jack315

下一步，看能否给双椭圆台相交曲线的参数方程？

在椭圆台参数方程中令 \(\sin{\theta}=t\)，则 \(\cos{\theta} =\pm\sqrt{1-t^2}\) ，即将椭圆台分成两半。
然后用相同的方法求解双椭圆台相交曲线的参数方程。
CPU 跑了 N 分钟后给出了结果。表达式很长很长。
俺的 CPU 是跑不动了，LZ 如果有资源可以自己试试。

Jack315

两个椭圆台的相交线位于两个平面上，因此可以改成求一个椭圆台与两个平面的相交线：
1. 用数值方法求出相交线上点的坐标。
2. 用得到的点坐标数据拟合平面。
3. 求椭圆台与两个平面相交线的参数方程。

这是取 \(h_{22}=500\) 时两个相交线的参数方程：
\(\begin{cases}
fx(t)=\frac{\cos{t}(23082.3-759.109\sin{t})-3.63798\times 10^{-12}\sin{t}}{(-15.2036+\cos{t})(-30.4071+\sin{t})}\\
fy(t)=\frac{(27705.3-911.144\sin{t})\sin{t}+\cos{t}(3.63798\times 10^{-12}-304.071\sin{t}+10\sin^2{t})}{(-15.2036+\cos{t})(-30.4071+\sin{t})}\\
fz(t)=\frac{-1083.67+\cos{t}(50678.5-1666.67\sin{t})+35.6388\sin{t}}{(-15.2036+\cos{t})(-30.4071+\sin{t})}
\end{cases}\)

\(\begin{cases}
fx(t)=\frac{-2.27374\times 10^{-13}\sin^2{t}+\cos{t}(102565+1600\sin{t}-2.66454\times 10^{-15}\sin^2{t})}{(32.0515+\cos{t}+2.22045\times 10^{-16}\sin{t})(64.103+\sin{t})}\\
fy(t)=\frac{\sin{t}[123111+1920.51\sin{t}+\cos{t}(641.03+10\sin{t})]}{(32.0515+\cos{t}+2.22045\times 10^{-16}\sin{t})(64.103+\sin{t})}\\
fz(t)=\frac{-5512.43+\cos{t}(-106838-1666.67\sin{t})-85.9934\sin{t}-2.27374\times 10^{-13}\sin^2{t}}{(32.0515+\cos{t}+2.22045\times 10^{-16}\sin{t})(64.103+\sin{t})}
\end{cases}
\)

可视化验证：


这是 Mathematica 代码文件供参考：
参数方程.rar (24.96 KB, 下载次数: 1)

2024-7-27 21:11 上传

点击文件名下载附件

文档中画的图在保存前已被删除以压缩文件大小。重新计算一遍就可以得到这些图。

数学星空

根据楼上的计算对比知道，2#的方法是因为误认为\(t_1=t_2=t\)，所以算出来的不对

重新计算如下：

设两椭圆台曲面A和B的参数方程为：

[(-30*k1 + 80)*cos(t1), (-30*k1 + 90)*sin(t1), 1000 - 1000*k1],

[-500*(-1 + k2)*sqrt(3)/2 + 5*sin(t2)*(-8 + k2), (-10*k2 + 70)*cos(t2) + 30 - 30*k2, -5*sin(t2)*(-8 + k2)*sqrt(3) - 500*(-1 + k2)/2]

隐函数方程为：

(10000*y^2)/(9*(2000 + z)^2) + 10000*x^2/(5000 + 3*z)^2 - 1,

2500*(sqrt(3)*x - 3*z)^2/(sqrt(3)*z + 7000*sqrt(3) + 3*x)^2 + (100*sqrt(3)*y - 3*sqrt(3)*z - 9*x)^2/(sqrt(3)*z + 6000*sqrt(3) + 3*x)^2 - 1

由\(x_1=x_2,y_1=y_2,z_1=z_2\)消元得到：

[-1000*(13*cos(t2) - 3*sin(t1) - 15)*(9*cos(t1) + 100*sqrt(3))*cos(t1)/((30000 + 27*(cos(t2) + 3)*cos(t1)^2 + 10000*cos(t2) - 30000*sin(t1))*sqrt(3) + 1800*(cos(t2) - (3*sin(t1))/2 + 3)*cos(t1)),

9*(((cos(t2) + 3)*cos(t1)^2 - (40000*cos(t2))/9 + 20000/3)*sqrt(3) + cos(t1)*(-(1100*cos(t2))/3 + 700))*sin(t1)/(10*((cos(t2) - 3*sin(t1) + 3)*sqrt(3) + (9*cos(t1)*(cos(t2) + 3))/100)*(sqrt(3)*cos(t1) + 100/3)),

(((-45000*cos(t2) - 135000)*cos(t1)^2 - 60000000*cos(t2) + 60000000*sin(t1))*sqrt(3) - 6900000*(cos(t2) - (18*sin(t1))/23 + 15/23)*cos(t1))/(300*((cos(t2) - 3*sin(t1) + 3)*sqrt(3) + (9*cos(t1)*(cos(t2) + 3))/100)*(sqrt(3)*cos(t1) + 100/3))]

其中约束方程为：

sin(t2) = 50*(-5600*cos(t1)*(cos(t2) - (51*sin(t1))/56 + 15/14)*sqrt(3) + 9*(-60 - 2*cos(t2) - 3*sin(t1))*cos(t1)^2 - 180000*cos(t2) + 180000*sin(t1))/((210000 + 27*(21 + cos(t2) + sin(t1))*cos(t1)^2 + 10000*cos(t2) - 150000*sin(t1))*sqrt(3) + 1800*cos(t1)*(cos(t2) - 7*sin(t1) + 21))

画图得到相交曲线图像：




将曲面A,B和相交曲线图像得到：



显然相交曲线计算的不完整，可能是因为约束方程的数值计算还要优化一下

现在剩下的问题是如何简化约束方程，使参数表达更简洁？是否有其它的更简洁的表达方案（参数方程）？

Jack315

“标准椭圆”参数方程：
\(\begin{cases}
f_x(t)=a\cos(t)\\
f_y(t)=b\sin(t)\\
f_z(t)=0
\end{cases}\)

将上述椭圆曲线绕 y 轴旋转一个角度 \(\theta\)，然后再移位
就能得到所需的椭圆台相交曲线的参数方程。
所以应该有“标准”的参数方程形式，或许更简洁。

参考：三维坐标变换--旋转矩阵与旋转向量
https://blog.csdn.net/mightbxg/article/details/79363699
\(\begin{array}{cc}
R_x=\begin{pmatrix}1&0&0\\0&\cos\theta&-\sin\theta\\0&\sin\theta&\cos\theta\end{pmatrix}&
R_y=\begin{pmatrix}\cos\theta&0&\sin\theta\\0&1&0\\-\sin\theta&0&\cos\theta\end{pmatrix}&
R_z=\begin{pmatrix}\cos\theta&-\sin\theta&0\\\sin\theta&\cos\theta&0\\0&0&1\end{pmatrix}
\end{array}\)

Jack315

用三角函数级数拟合得到的参数方程：
\(\begin{cases}
fx(t)=2.1537+50.1492\sin(1.0021t-0.0075)+1.1724\sin(1.9642t-1.3388)\\
fy(t)=-0.0549+60.0152\sin(1.0009t+1.5679)+0.9832\sin(2.0004t+0.0364)\\
fz(t)=2.3816+109.8762\sin(1.0006t-0.0021)+1.8069\sin(1.9895t-1.4882)
\end{cases}\)

\(\begin{cases}
fx(t)=-1.5424-49.972sin(1.0010t+43.9793)\\
fy(t)=0.0092+59.9685\sin(0.9998t+1.5713)+0.1729\sin(2.0062t-3.1244)\\
fz(t)=-1.0723+51.8990\sin(0.9987t+0.0052)+0.2599\sin(1.8436t-3.8483)
\end{cases}\)

拟合参数方程中二次谐波的出现，意味着曲线不是完全在一个平面上。
只是与基频比较，系数小了约一个数量级，所以看起来像是在一个平面上。

另外，用 Chebyshev 多项式拟合也是一个可能的选项，近似计算可能会更方便。

mathe

看起来像在一个平面是因为两个图看起来像两个椭圆锥相交，而且椭圆的长短轴位置对应在一起。对于这种情景，查看过四个相交在一起的长短轴顶点的平面，可以知道这个平面割两个椭圆锥的结果都是椭圆，而且这个相交出来的椭圆长短轴相等，得到两个椭圆轨迹重叠；所以两者相交在这个平面上。

Jack315

使用切比雪夫多项式拟合三角函数时，从拟合精度角度看，宜选择
第一类切比雪夫多项式 \(ChebyshevT(n,x)\) 拟合 \(\sin(x)\) 函数。
相关数据如下表所示：
\(\begin{array}{|c|c|c|}
\hline T(n,x)\text{ 的 }n\text{ 参数}&\text{多项式项数}&\text{拟合 }\sin(x)\text{ 的最大误差}\\
\hline 3&2&0.020017\\
\hdashline 5&3&0.007117\\
\hdashline 7&4&0.003619\\
\hdashline 9&5&0.002186\\
\hdashline 11&6&0.001462\\
\hdashline 13&7&0.001047\\
\hdashline 15&8&0.000786\\
\hdashline ...&...&...\\
\hdashline 43&22&0.000096\\
\hline \end{array}\)

用第一类切比雪夫多项式 \(ChebyshevT(13,x)\) 拟合 \(\sin(x)\) 的函数如下：

1. ChebyshevTSin[x_] := Module[
   
2.   {sol, r, a, t, sign},
   
3.   
   
4.   sol = NSolve[D[ChebyshevT[13, r], r] == 0, r] // Sort;
   
5.   a = r /. sol[[7]];
   
6.   
   
7.   t = Mod[Abs[x], 2 \[Pi]];
   
8.   (*确定正弦函数值的符号*)
   
9.   sign = Sign[x] If[0 < t < \[Pi], 1, -1];
   
10.   
    
11.   (*将 t 移到 [0,\[Pi]/2] 范围内对应的位置*)
    
12.   t = (2 a)/\[Pi] If[t < \[Pi]/2, t,
    
13.      If[t < \[Pi], \[Pi] - t,
    
14.       If[t < (3 \[Pi])/2, t - \[Pi], 2 \[Pi] - t]]];
    
15.   
    
16.   (*输出拟合值*)
    
17.   sign  ChebyshevT[13, t]
    
18.   ](*End Module*)

复制代码

将 \(ChebyshevTSin(x)\) 作为 \(\sin(x)\) 代入 16# 的公式即可。