椭圆台曲面被斜截面截交的曲线求解

数学星空

关于8#网友Huxley 的推导过程中省略的部分做一点计算和说明：

\( \cos(d\phi)\)=\(\frac{r\*r+r\*r'd\theta}{\sqrt{r\*r*(r\*r+2r\*r'd\theta+(r'\*r')(d\theta)^2}}=\frac{1+xd(\theta)}{\sqrt{(1+xd(\theta))^2+(y-x^2)d(\theta)^2}}\)

其中\(x=\frac{r\*r'}{r\*r}\),\(y=\frac{r'\*r'}{r\*r}\),\(A=\sqrt{y-x^2}\)

\(d\phi=\arccos( \frac{1}{\sqrt{1+(\frac{Ad(\theta)}{1+xd(\theta)})^2}})d(\theta) \)

对右边关于\(\theta\)渐近展开得到：

\(d\phi=Ad\theta+Axd(\theta)^2+\frac{1}{3}A(-3x^2+A^2)(d(\theta))^3-xA(-x^2+A^2)(d(\theta))^4+....\)

取第一项得到：

\(d\phi=Ad\theta=\sqrt{y-x^2}d(\theta)=\sqrt{\frac{r'\*r'}{r\*r}-(\frac{r\*r'}{r\*r})^2}d(\theta)\)

对于网友Jack315 15#的计算结果：

\(\theta=\int_0^t(\sqrt{\frac{a^2\sin(t)^2+b^2\cos(t)^2}{a^2\cos(t)^2+b^2\sin(t)^2+c^2}-(\frac{(b^2-a^2)\sin(t)*\cos(t)}{a^2\cos(t)^2+b^2\sin(t)^2+c^2})^2}\dif t\)

可以利用级数展开得到：

\(\theta=\frac{b}{\sqrt{a^2+c^2}}\int_0^t(1+A_1t^2+A_2t^4+A_3t^6+...)\dif t\)

其中

\(A_1=\frac{(2b^2 + c^2)(a - b)(a + b)t^2}{2*(a^2 + c^2)*b^2}\)

\(A_2=\frac{(a + b)(a - b)(16a^2b^4 + 8a^2b^2c^2 - 3a^2c^4 - 24b^6 - 20b^4c^2 - b^2c^4)}{24(a^2 + c^2)^2b^4}\)

\(A_3=\frac{(a + b)(a - b)(272a^4b^6 + 136a^4b^4c^2 - 30a^4b^2c^4 + 45a^4c^6 - 960a^2b^8 - 896a^2b^6c^2 - 88a^2b^4c^4 - 30a^2b^2c^6 + 720b^10 + 840b^8c^2 + 182b^6c^4 + b^4c^6)}{720(a^2 + c^2)^3b^6}\)

对于｛a=3,b=2,c=5}

\(F_1(t)=\frac{b}{\sqrt{a^2+c^2}}\int_0^t(1+A_1t^2+A_2t^4+A_3t^6+...)\dif t=\sqrt(\frac{2}{17})(t+\frac{165}{3*272}t^3-\frac{32345}{5*147968}t^5+\frac{3380957}{7*40247296}t^7+...)\)

\(F_2(t)=\int_0^t(\sqrt{\frac{9\cos(t)^2+4\sin(t)^2}{9\cos(t)^2+4\sin(t)^2+25}-(\frac{25\sin(t)\cos(t)}{9\cos(t)^2+4\sin(t)^2+25})^2}\dif t\)

取20个样本点\(t=\frac{\pi}{2}*\frac{k}{20},k=1..20\)

得到：

\([F_1(t),F_2(t)]=[0.02697249116, 0.02697249116], [0.05414525431, 0.05414525413], [0.08171331124, 0.08171330529], [0.1098615660, 0.1098614892], [0.1387607016, 0.1387601509], [0.1685642254, 0.1685615051], [0.1994070420, 0.1993966715], [0.2314059407, 0.2313732537], [0.2646623758, 0.2645733101], [0.2992679249, 0.2990515907], [0.3353128063, 0.3348339855], [0.3728978377, 0.3719162088], [0.4121502213, 0.4102628287], [0.4532435323, 0.4498068149], [0.4964223004, 0.4904498159], [0.5420315594, 0.5320633799], [0.5905517542, 0.5744912936], [0.6426393807, 0.6175531354], [0.6991737474, 0.6610490316], [0.7613102389, 0.7047654702]\)

显然，当\(t<1\)时结果还不错，但大于1时，偏差有点大，需要进一步改进

数学星空

关于用三角级数展开上面的积分式，很奇怪得到答案偏差很大：

@wayne  帮忙重新算算，是否是同样的结论？

@mathe  帮忙分析一下，如何渐近展开下面积分式，使其计算更精确？

\(F_1(t)=\int_0^t(\sqrt{\frac{a^2\sin(t)^2+b^2\cos(t)^2}{a^2\cos(t)^2+b^2\sin(t)^2+c^2}-(\frac{(b^2-a^2)\sin(t)*\cos(t)}{a^2\cos(t)^2+b^2\sin(t)^2+c^2})^2}\dif t\)

首先使用三角代换\(\sin(t)=x,\cos(t)=\sqrt{1-x^2}\),并关于\(x\)渐近展开后再代入\(x=\sin(t)\)并积分得到：

F_2(t)=m + ((2*b^2 + c^2)*m*(a - b)*(a + b)*(t/2 - sin(2*t)/4))/(2*b^2*(a^2 + c^2)) + m*(a + b)^2*(8*b^4 + 4*b^2*c^2 - c^4)*(a - b)^2*((3*t)/8 + sin(4*t)/32 - sin(2*t)/4)/(8*b^4*(a^2 + c^2)^2) + m*(a + b)^3*(16*b^6 + 8*b^4*c^2 - 2*b^2*c^4 + c^6)*(a - b)^3*((5*t)/16 + (3*sin(4*t))/64 - (15*sin(2*t))/64 - sin(6*t)/192)/(16*(a^2 + c^2)^3*b^6) + ((128*b^8 + 64*b^6*c^2 - 16*b^4*c^4 + 8*b^2*c^6 - 5*c^8)*m*(a + b)^4*(a - b)^4*((35*t)/128 + (7*sin(4*t))/128 - (7*sin(2*t))/32 + sin(8*t)/1024 - sin(6*t)/96))/(128*(a^2 + c^2)^4*b^8) + ((256*b^10 + 128*b^8*c^2 - 32*b^6*c^4 + 16*b^4*c^6 - 10*b^2*c^8 + 7*c^10)*m*(a + b)^5*(a - b)^5*((63*t)/256 + (15*sin(4*t))/256 - (105*sin(2*t))/512 + (5*sin(8*t))/2048 - (15*sin(6*t))/1024 - sin(10*t)/5120))/(256*(a^2 + c^2)^5*b^10)

其中\(m=\frac{b}{\sqrt{a^2+c^2}}\)

对于｛a=3,b=2,c=5}

\(F_2(t)=0.3429971702 + 0.1061429486t - 0.05377977603\sin(2t) + 0.0004682430182\sin(4t) - 0.00008154587352\sin(6t) + 5.355467693*10^{-6}\sin(8t) - 9.937085131*10^{-7}\sin(10t)\)

\(F_1(t)=\int_0^t(\sqrt{\frac{9\cos(t)^2+4\sin(t)^2}{9\cos(t)^2+4\sin(t)^2+25}-(\frac{5\sin(t)*\cos(t)}{9\cos(t)^2+4\sin(t)^2+25})^2}\dif t\)

取20个样本点\(t=\frac{\pi}{2}*\frac{k}{20},k=1..20\)

得到：

\([F_1(t),F_2(t)]=[0.3430307266, 0.02697249116], [0.3432645548, 0.05414525413], [0.3438936713, 0.08171330529], [0.3451029204, 0.1098614892], [0.3470626483, 0.1387601509], [0.3499250760, 0.1685615051], [0.3538213552, 0.1993966715], [0.3588592096, 0.2313732537], [0.3651210441, 0.2645733101], [0.3726624232, 0.2990515907], [0.3815108450, 0.3348339855], [0.3916647344, 0.3719162088], [0.4030925949, 0.4102628287], [0.4157323021, 0.4498068149], [0.4294906392, 0.4904498159], [0.4442433457, 0.5320633799], [0.4598361191, 0.5744912936], [0.4760870785, 0.6175531354], [0.4927910990, 0.6610490316], [0.5097261239, 0.7047654702]\)

Jack315

【建模】
公式：
\(\rho(t)=\lambda\sqrt{a^2\cos^2t+b^2\sin^2t+c^2}\)
\(\theta(t)=\int_0^t{\sqrt{\frac{a^2\cos^2t+b^2\sin^2t}{a^2\cos^2t+b^2\sin^2t+c^2}-[\frac{(b^2-a^2)\sin t\cos t}{a^2\cos^2t+b^2\sin^2t+c^2}]^2}}\dif t\)
令：
\(f(t)=\sqrt{\frac{a^2\cos^2t+b^2\sin^2t}{a^2\cos^2t+b^2\sin^2t+c^2}-[\frac{(b^2-a^2)\sin t\cos t}{a^2\cos^2t+b^2\sin^2t+c^2}]^2}\)
其中：\(0\le\lambda\le1, 0\le t\le 2\pi\) 。
给定参数：
\(\lambda=1, a=3, b=2, c=5\)
\(\rho(t)\) 和 \(f(t)\) 的函数曲线如下图所示：

注：周期为 \(\pi\) 。
这两个函数都可以用正弦函数来拟合，
而且只用一次和两次谐波就能得达到一定的精度 (1E-4)。
但如要展开成三角级数却并不容易，
因为三角级数系数的计算涉及到了函数在一个周期内的积分，
而这正是现在所面临的难题。

但办法仍然是有的：实验设计建模。
假设 \(\rho(t)\) 和 \(f(t)\) 的函数的形式确定，只是参数 \(\lambda,a,b,c\) 会发生变化。

• 确定参数  \(\lambda,a,b,c\)  的变化范围。
• 在参数  \(\lambda,a,b,c\)  的变化范围内，用正交试验的方法确定一组参数 \(Xs\)。
• 对于每一组参数 \(Xs\)进行”实验“，得到一组拟合参数 \(Ys\)。
• 根据实验结果，用 \(Xs\) 和 \(Ys\) 建立（多项式）模型。
  

得到的模型就是用参数  \(\lambda,a,b,c\)  计算拟合参数的表达式。
注：所谓的”实验“，就是用数值积分计算三角级数的系数。

Jack315

参考 23#
建模中的一次“实验”示例：

\(f_\rho(t)=\lambda\sqrt{a^2\cos^2t+b^2\sin^2t+c^2}\)
\(f_\theta(t)=\sqrt{\frac{a^2\cos^2t+b^2\sin^2t}{a^2\cos^2t+b^2\sin^2t+c^2}-[\frac{(b^2-a^2)\sin t\cos t}{a^2\cos^2t+b^2\sin^2t+c^2}]^2}\)
不难看出这两个函数都是周期偶函数，且周期为 \(T=\pi\) 。

傅立叶级数的公式为：
\(f(t)=\frac{a_0}{2}+\sum_{i=0}^{\infty}[a_n\cos(n\omega t)+b_n\sin(n\omega t)]\)
其中：
\(\omega=\frac{2\pi}{T}\)
\(a_0=\frac{2}{T}\int_0^Tf(t)\dif t\)
\(a_n=\frac{2}{T}\int_0^Tf(t)\cos(n\omega t)\dif t\)
\(b_n=\frac{2}{T}\int_0^Tf(t)\sin(n\omega t)\dif t\)
\(n=1,2,3,...\)
两个函数的傅立叶级数系数 \(b_n=0, n=1,2,3,...\)

给定参数 \(\lambda=1,a=3,b=2,c=5\)，
计算得到 \(a_n, n=0,1,2,3,...\) 的系数为……
\(\rho\)：
11.220546616558773

0.222849522518958
-0.00221386367935080
4.39909085014873e-05
-1.09270306744829e-06
3.03996077794157e-08

\(\theta\)：
0.890194089661798

0.0715502824818127
-0.00216353849418902
3.76183728794146e-06
9.52612614705229e-06
-1.48434981671988e-06

若取 \(n=0,1,2\) ，则
\(\rho(t)=f_\rho(t)\approx 5.61+0.22\sin(2t)-0.0022\sin(4t)\)

\(f_\theta(t)\approx 0.45+0.072\sin(2t)-0.0022\sin(4t)\)
再求积分 \(\int_0^t f_\theta(t)\dif t\) 得：
\(\theta(t)\approx0.45t-0.3545-0.036\cos(2t)+0.00055\cos(4t)\)
说明：

• 上述三角级数系数未经验证。
• 也可以在得到表达式模型后再求积分。
  

至此，由椭圆锥体参数 \(Xs\): (1,3,2,5) 得到拟合参数 \(Ys\): (5.61, 0.22, -0.0022), (0.45,-0.3545,-0.036, 0.00055) 。
如果有 N 组这样的实验数据，则有可能得到模型 \(Ys=f(Xs)\) ，
即 7 个用椭圆锥体参数构成的（多项式）表达式，用于计算 7 个拟合参数。

Jack315

参考 24#

【建模示例】
由椭圆锥体参数\((\lambda,a,b,c)\)，计算三角级数系数\((a_0,a_1,a_2)\)的多项式模型。
模型质量未验证，仅为展示建模方法而写的代码。

Matlab 代码：
FitModel.rar (124.89 KB, 下载次数: 4)

2024-3-25 22:43 上传

点击文件名下载附件

数学星空

没想到，maple 求出了积分表达式，然而表达式包含了许多根号和椭圆函数，让人望而生畏~对下一步简化工作，完全没有信心啦~



表达式见下：

Jack315

表达式在 Matlab 文档里：

其中：\(\lambda, a, b, c\) 是椭圆锥体的参数。

这些表达式不一定正确或者最佳，
但结果的形式就是这样的多项式表达式。
不太清楚这个是不是所要的形式？

Jack315

理一下思路，确保大家是在同一个方向使劲……

给定曲面参数方程：
\(r(t)=\{\lambda a\cos t,\lambda b\sin t, \lambda c\}, 0\le\lambda\le1,0\le t\le2\pi\) 。
得到展开曲线方程：
\(\rho(t)=\lambda\sqrt{a^2\cos^2t+b^2\sin^2t+c^2}=\lambda f_1(t)\)
\(\theta(t)=\int_0^t{\sqrt{\frac{a^2\cos^2t+b^2\sin^2t}{a^2\cos^2t+b^2\sin^2t+c^2}-[\frac{(b^2-a^2)\sin t\cos t}{a^2\cos^2t+b^2\sin^2t+c^2}]^2}}\dif t=\int_0^tf_2(t)\dif t\)
由于计算能力的限制，\(f_1(t), f_2(t)\) 都无法进行实时计算。
因此，再将 \(f_1(t), f_2(t)\) 展开为傅立叶级数：
\(f_1(t)\approx\frac{a_{10}}{2}+a_{11}\cos(n\omega t)+a_{12}\cos(n\omega t)\)
\(f_2(t)\approx\frac{a_{20}}{2}+a_{21}\cos(n\omega t)+a_{22}\cos(n\omega t)\)
由于存在快速的 Cordic 算法用于计算三角函数，因此假设到这一步计算能力已能满足要求。
现在的问题是要找出 \((a_{10},a_{11},a_{12})=F_1(\lambda,a,b,c),(a_{20},a_{21},a_{22})=F_2(\lambda,a,b,c)\) 的表达式。
我们通过建模的方法找出在计算能力范围内的 \(F_1,F_2\) ，一般为多项式表达式。
这个思路方向正确吗？

如果这个思路确实能解决计算能力的问题，建议 LZ 给出一个具体的、接近实际情况的设计规范，包括：

• 设计输入——曲面参数方程及各参数的变化范围（上、下限）。
• 设计输出——展开曲线各参数的误差范围（上、下限）。
  

建模（优化）过程比较繁琐费时，模型也有多种选择。
如果这个思路方向正确的话，我再花时间做一个能满足精度要求的模型，
并为 LZ 提供一个相对完整的文档供参考。LZ 意下如何？

Jack315

【误差估计】

傅立叶级数截尾误差估计在 1E-5 数量级：


模型响应曲面：
椭圆锥体参数 a =10, e = 0.1~0.9, c = 1~100 。

从图上可以看出响应曲面相对“光滑”，但模型函数的形式可能会比较复杂。
如果椭圆锥体参数变化范围比较小，则用多项式模型就能得到比较理想的结果。

接下来先在几个不同位置但椭圆锥体参数范围较小的区域建模，
然后比较函数 \(f_1(t), f_2(t)\) 的圴方根误差 (RMSE) ，对总体误差做个初步估计。
如有需要，再寻找椭圆锥体参数范围较大时的模型形式。

Jack315

CORDIC：https://baike.so.com/doc/7592544-7866639.html
CORDIC 算法可用于三角函数和双曲函数的计算，
因此假设所有基本初等函数的计算都不受计算能力的限制。
基本初等函数：https://baike.so.com/doc/3863106-4055726.html

• 常数
• 幂函数
• 指数函数
• 对数函数
• 三角函数
• 反三角函数
  

如果假设成立，并且在代码中能实现数值积分，则可以直接对 \(\rho(t),\theta(t)\) 做“精确”计算。
反之如果假设不成立，则建模所能使用的函数也会有所限制。
目前估计可使用的基本初等函数包括：

• 常数
• 幂函数（幂次为整数）
• 三角函数
• 反三角函数
• 双曲函数
• 反双曲函数
  

在椭圆锥体全体参数范围内，建模仍有可能找到最佳的拟合函数形式。

建议 LZ 先确定以下问题，再决定下一步的方向：

• 基本初等函数的计算是否有受限部分，如有则不受限的函数有哪些？
• 如果能计算所有的基本初等函数，数值积分计算是否受限？
  

从计算精度角度来看，能直接“精确”计算 \(\rho(t),\theta(t)\)  函数为佳；
从计算速度角度来看，以放弃部分精度，计算三角级数或使用插值的方案为佳。

另外开根号可以用迭代求根的算法，是否有更高效的算法？
能不能先尝试数值积分直接“精确”计算 \(\rho(t),\theta(t)\)  函数？