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

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)\)  函数？

Jack315

【数值积分】

自己写插值、数值积分……并不真的需要完全自己动手写。
如果编程语言是 C/C++，则有一些工具就可以帮助完成这些工作，
如 Matlab 就可以生成数值积分函数：数值积分。
除了 Matlab ，或许还有其他工具，如 Mahematica，Maple, Python ……
根据自己手上的资源，逐个查验一下。

注：网页底部显示 “C/C++ 代码生成“表示可以生成 C/C++ 代码。
Matlab 是基于数组的形式进行计算的，可以修改成单个变量。

所以问题归结到当前的计算平台：

• 有没有足够的内存？
• 运行速度是否够快？
  

如果当前的计算平台资源不够，能否升级当前的计算平台？

数学星空

对于

\(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\)

做第一步代换：

m1 = -1/2*a^2*c^2 + 1/2*b^2*c^2,
n1 = a^2*b^2 + a^2*c^2/2 + b^2*c^2/2,
m2 = a^2/2 - b^2/2,
n2 = b^2/2 + c^2 + a^2/2

可以将原积分式化简为：

\(\frac{\sqrt{m1*\cos(2*t)+n1}}{m2*\cos(2*t)+n2}\)


然后积分得到：

\(a0=\sqrt{n1-m1}m2(m2-n2)=\frac{a(b^2+c^2)^{3/2}(a^2-b^2)}{2}\)

\(a1=m1n2-m2n1=-\frac{(b^2+c^2)(a^2+c^2)(a^2-b^2)}{2}\)

\(a2=m1(m2-n2)=\frac{c^2(a^2-b^2)(b^2+c^2)}{2}\)

\(x=\cos(t)\)

\(b1=\frac{2m2}{m2-n2}=\frac{(b^2-a^2)(b^2+c^2)}{2}\)

\(c1=\sqrt{\frac{2m1}{m1-n1}}=\frac{c}{a}\frac{\sqrt{a^2-b^2}}{b^2+c^2}\)

最终得到：

\(F1(t)=\frac{a1(EllipticPi(x,b1,c1)-EllipticPi(b1,c1))+a2(EllipticPi(x,c1)-EllipticK(x,c1))}{a0}\)

例如：

\(a = 3, b = 2, c = 5, m1 = -125/2, m2 = 5/2, n1 = 397/2, n2 = 63/2,a0 = (435*sqrt(29))/2, a1 = -2465, a2 = 3625/2, b1 = -5/29, c1 = (5*sqrt(145))/87, x = cos(t)\)

代入得到：

F2(t)=-2.104547165*EllipticPi(cos(t), -0.1724137931, 0.6920456655) + 1.547461151*EllipticF(cos(t), 0.6920456655) + 0.7047654707

数值计算：\(t =\frac{\Pi}{2}\frac{k}{20},k=1..20\)

[F1,F2]=[0.02697249116, 0.0269724907], [0.05414525413, 0.0541452537], [0.08171330529, 0.0817133067], [0.1098614892, 0.1098614897], [0.1387601509, 0.1387601517], [0.1685615051, 0.1685615067], [0.1993966715, 0.1993966737], [0.2313732537, 0.2313732547], [0.2645733101, 0.2645733107], [0.2990515907, 0.2990515917], [0.3348339855, 0.3348339867], [0.3719162088, 0.3719162097], [0.4102628287, 0.4102628296], [0.4498068149, 0.4498068155], [0.4904498159, 0.4904498166], [0.5320633799, 0.5320633804], [0.5744912936, 0.5744912942], [0.6175531354, 0.6175531360], [0.6610490316, 0.6610490321], [0.7047654702, 0.7047654707]

级数展开：

F2( x)=-(a1*EllipticPi(b1, c1) + EllipticK(c1)*a2)/a0 + (a1 + a2)*x/a0 + (a1*(c1^2/6 + b1/3 + 1/6) + (c1^2/6 + 1/6)*a2)*x^3/a0 + (a1*(3/40*c1^4 + 1/10*b1*c1^2 + 1/5*b1^2 + 1/20*c1^2 + 1/10*b1 + 3/40) + (3/40*c1^4 + 1/20*c1^2 + 3/40)*a2)*x^5/a0 + (a1*(5/112*c1^6 + 3/56*b1*c1^4 + 1/14*b1^2*c1^2 + 3/112*c1^4 + 1/7*b1^3 + 1/28*b1*c1^2 + 1/14*b1^2 + 3/112*c1^2 + 3/56*b1 + 5/112) + (5/112*c1^6 + 3/112*c1^4 + 3/112*c1^2 + 5/112)*a2)*x^7/a0 + (a1*(35/1152*c1^8 + 5/144*b1*c1^6 + 1/24*b1^2*c1^4 + 5/288*c1^6 + 1/18*b1^3*c1^2 + 1/48*b1*c1^4 + 1/9*b1^4 + 1/36*b1^2*c1^2 + 1/64*c1^4 + 1/18*b1^3 + 1/48*b1*c1^2 + 1/24*b1^2 + 5/288*c1^2 + 5/144*b1 + 35/1152) + (35/1152*c1^8 + 5/288*c1^6 + 1/64*c1^4 + 5/288*c1^2 + 35/1152)*a2)*x^9/a0 + (a1*(63/2816*c1^10 + 35/1408*b1*c1^8 + 5/176*b1^2*c1^6 + 35/2816*c1^8 + 3/88*b1^3*c1^4 + 5/352*b1*c1^6 + 1/22*b1^4*c1^2 + 3/176*b1^2*c1^4 + 15/1408*c1^6 + 1/11*b1^5 + 1/44*b1^3*c1^2 + 9/704*b1*c1^4 + 1/22*b1^4 + 3/176*b1^2*c1^2 + 15/1408*c1^4 + 3/88*b1^3 + 5/352*b1*c1^2 + 5/176*b1^2 + 35/2816*c1^2 + 35/1408*b1 + 63/2816) + (63/2816*c1^10 + 35/2816*c1^8 + 15/1408*c1^6 + 15/1408*c1^4 + 35/2816*c1^2 + 63/2816)*a2)*x^11/a0 + (a1*(9/832*b1^2*c1^4 + 5/416*b1^2*c1^6 + 3/208*b1^3*c1^2 + 1/26*b1^5*c1^2 + 5/208*b1^3*c1^6 + 5/416*b1^2*c1^2 + 35/3328*b1*c1^2 + 1/52*b1^4*c1^2 + 3/104*b1^4*c1^4 + 15/1664*b1*c1^6 + 15/1664*b1*c1^4 + 105/13312*c1^8 + 25/3328*c1^6 + 63/3328*b1*c1^10 + 35/3328*b1*c1^8 + 3/208*b1^3*c1^4 + 63/6656*c1^10 + 35/1664*b1^2*c1^8 + 105/13312*c1^4 + 63/6656*c1^2 + 1/26*b1^5 + 35/1664*b1^2 + 5/208*b1^3 + 1/13*b1^6 + 3/104*b1^4 + 63/3328*b1 + 231/13312*c1^12 + 231/13312) + (231/13312*c1^12 + 63/6656*c1^10 + 105/13312*c1^8 + 25/3328*c1^6 + 105/13312*c1^4 + 63/6656*c1^2 + 231/13312)*a2)*x^13/a0 + (a1*(35/6144*c1^8 + 35/6144*c1^6 + 1/128*b1^2*c1^4 + 7/1024*b1*c1^8 + 1/128*b1^2*c1^6 + 1/96*b1^3*c1^2 + 3/320*b1^3*c1^4 + 143/10240*c1^14 + 63/10240*c1^4 + 21/1280*b1^2 + 7/384*b1^3 + 1/48*b1^4 + 5/768*b1*c1^6 + 1/96*b1^3*c1^6 + 7/768*b1^2*c1^8 + 1/15*b1^7 + 1/40*b1^5 + 77/10240*c1^12 + 77/5120*b1 + 143/10240 + 1/30*b1^6 + 7/384*b1^3*c1^8 + 21/2560*b1*c1^10 + 1/30*b1^6*c1^2 + 21/2560*b1*c1^2 + 1/60*b1^5*c1^2 + 21/1280*b1^2*c1^10 + 7/768*b1^2*c1^2 + 63/10240*c1^10 + 1/80*b1^4*c1^4 + 7/1024*b1*c1^4 + 1/48*b1^4*c1^6 + 77/5120*b1*c1^12 + 1/40*b1^5*c1^4 + 1/80*b1^4*c1^2 + 77/10240*c1^2) + (143/10240*c1^14 + 77/10240*c1^12 + 63/10240*c1^10 + 35/6144*c1^8 + 35/6144*c1^6 + 63/10240*c1^4 + 77/10240*c1^2 + 143/10240)*a2)*x^15/a0

代入上面的相关参数，得到

[F1,F2]=[0.02697249116, 0.07122775680], [0.05414525413, 0.08127508887], [0.08171330529, 0.09721574295], [0.1098614892, 0.1180523822], [0.1387601509, 0.1427261684], [0.1685615051, 0.1703034082], [0.1993966715, 0.2000823349], [0.2313732537, 0.2316116508], [0.2645733101, 0.2646452166], [0.2990515907, 0.2990699796], [0.3348339855, 0.3348378550], [0.3719162088, 0.3719168525], [0.4102628287, 0.4102629092], [0.4498068149, 0.4498068222], [0.4904498159, 0.4904498166], [0.5320633799, 0.5320633806], [0.5744912936, 0.5744912943], [0.6175531354, 0.6175531358], [0.6610490316, 0.6610490323], [0.7047654702, 0.7047654708]

这个级数展开在t较小时，误差很大，对于基本的椭圆积分，应该有相关渐近展开公式，有谁可以提供一下??  
@wayne

\(EllipticPi(z,v,k)=\int_0^z\frac{1}{(1-vt^2)\sqrt{(1-t^2)(1-k^2t^2)}}\dif t\)

\(EllipticPi(1,v,k)=EllipticPi(v,k)\)

\(EllipticF(z,k)=\int_0^z\frac{1}{\sqrt{(1-t^2)(1-k^2t^2)}}\dif t\)

\(EllipticF(1,k)=EllipticK(k)\)



以下内容为微信公众号：数学的情怀提供

Jack315

数学星空 发表于 2024-3-28 09:40
对于

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

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

\(\frac{\sqrt{m1\cos(2t)+n1}}{m2\cos(2t)+n2}=\sqrt{2}\frac{\sqrt{b^2c^2+a^2(2b^2+c^2)-(a^2-b^2)c^2\cos(2t)}}{a^2+b^2+2c^2+(a^2-b^2)\cos(2t)}\)

两个式子对不上，不知道是哪里搞错了。

其它的项可能有误，但余弦项应该是没对上。

数学星空

你的分子计算都不对，你可以再算一遍

Jack315

\(\cos^2t=\frac{1}{2}(1+\cos 2t)=\frac{1}{2}(1+x)\)
\(\sin^2t=\frac{1}{2}(1-\cos 2t)=\frac{1}{2}(1-x)\)
代入下式：
\(\sqrt{\frac{a^2\cos^2t+b^2\sin^2t}{a^2\cos^2t+b^2\sin^2t+c^2}-[\frac{(a^2-b^2)\sin(t)\cos(t)}{a^2\cos^2t+b^2\sin^2t+c^2}]^2}\)
\(=\sqrt{2}\frac{\sqrt{(a^2-b^2)^2x^2+(a^2-b^2)(a^2+b^2+c^2)x+(2a^2b^2+b^2c^2+c^2a^2)}}{(a^2-b^2)x+(a^2+b^2+2c^2)}\)

Matlab 代码验证：

1. f[t_] := Sqrt[((a^2 Cos[t]^2 + b^2 Sin[t]^2)/(a^2 Cos[t]^2 + b^2 Sin[t]^2 + c^2) - (((a^2 - b^2) Sin[t] Cos[t])/(a^2 Cos[t]^2 + b^2 Sin[t]^2 + c^2))^2)]
   
2. h[t_] := Sqrt[2] Sqrt[((a^2 - b^2)^2 x^2 + (a^2 - b^2) (a^2 + b^2 + c^2) x + (2 a^2 b^2 + b^2 c^2 + c^2 a^2))/((a^2 - b^2) x + (a^2 + b^2 + 2 c^2))^2]
   
3. f[t] - h[t] /. x -> Cos[2 t] // Simplify

复制代码

结果为 0 。

数学星空

为了得到较精确的渐近展开，我们放弃直接级数展开，利用逐层渐近展开计算方案：

对于

\(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\)

做第一步代换：

m1 = -1/2*a^2*c^2 + 1/2*b^2*c^2, n1 = a^2*b^2 + a^2*c^2/2 + b^2*c^2/2, m2 = a^2/2 - b^2/2, n2 = b^2/2 + c^2 + a^2/2

可以将原积分式化简为：

\(\frac{\sqrt{m1*\cos(2*t)+n1}}{m2*\cos(2*t)+n2}=p0\frac{\sqrt{k1*\cos(2*t)+1}}{k2*\cos(2*t)+1}\)

sqrt(1+k1x)=1 + 1/2*k1*x - 1/8*k1^2*x^2 + 1/16*k1^3*x^3 - 5/128*k1^4*x^4 + 7/256*k1^5*x^5 - 21/1024*k1^6*x^6 + 33/2048*k1^7*x^7 - 429/32768*k1^8*x^8 + 715/65536*k1^9*x^9            

1/(k2*x + 1)=-k2^9*x^9 + k2^8*x^8 - k2^7*x^7 + k2^6*x^6 - k2^5*x^5 + k2^4*x^4 - k2^3*x^3 + k2^2*x^2 - k2*x + 1

其中

\(k1=\frac{m1}{n1}=\frac{-a^2c^2 + b^2c^2}{2a^2b^2 + a^2c^2 + b^2c^2}\)

\(k2=\frac{m2}{n2}=\frac{a^2 - b^2}{a^2 + b^2 +2c^2}\)

\(p0=\frac{\sqrt{n1}}{n2}=\frac{\sqrt{4a^2b^2 +2a^2c^2 +2b^2c^2}}{a^2 + b^2 +2c^2}\)


（1）若我们直接利用根式与分式展开代入原积分式可以得到（由于表达式太长，我们只取关于k2^i,i=0..6 项，因为k2一般为零点几，对于精度在10^(-6)工程应用 高次项可以忽略）

F11(t)=p0*((-(315*sin(4*t))/131072 - (63*sin(8*t))/131072 - (7*sin(12*t))/131072 - (105*t)/16384)*k1^6 + ((35*sin(6*t))/24576 + (7*sin(10*t))/40960 + (35*sin(2*t))/4096)*k1^5 + (-(5*sin(4*t))/1024 - (5*sin(8*t))/8192 - (15*t)/1024)*k1^4 + (sin(6*t)/384 + (3*sin(2*t))/128)*k1^3 + (-sin(4*t)/64 - t/16)*k1^2 + k1*sin(2*t)/4 + t + ((-(21*sin(8*t))/32768 - (105*sin(4*t))/32768 - (35*t)/4096 - (7*sin(12*t))/98304)*k1^5 + ((25*sin(2*t))/2048 + sin(10*t)/4096 + (25*sin(6*t))/12288)*k1^4 + (-sin(8*t)/1024 - sin(4*t)/128 - (3*t)/128)*k1^3 + ((3*sin(2*t))/64 + sin(6*t)/192)*k1^2 + (-sin(4*t)/16 - t/4)*k1 - sin(2*t)/2)*k2 + ((-(75*sin(4*t))/16384 - (25*t)/2048 - (15*sin(8*t))/16384 - (5*sin(12*t))/49152)*k1^4 + (sin(10*t)/2560 + (5*sin(2*t))/256 + (5*sin(6*t))/1536)*k1^3 + (-sin(4*t)/64 - (3*t)/64 - sin(8*t)/512)*k1^2 + ((3*sin(2*t))/16 + sin(6*t)/48)*k1 + sin(4*t)/8 + t/2)*k2^2 + ((-(15*sin(4*t))/2048 - (3*sin(8*t))/2048 - sin(12*t)/6144 - (5*t)/256)*k1^3 + ((5*sin(2*t))/128 + (5*sin(6*t))/768 + sin(10*t)/1280)*k1^2 + (-sin(4*t)/16 - sin(8*t)/128 - (3*t)/16)*k1 - (3*sin(2*t))/8 - sin(6*t)/24)*k2^3 + ((-(15*sin(4*t))/1024 - (3*sin(8*t))/1024 - (5*t)/128 - sin(12*t)/3072)*k1^2 + (sin(10*t)/320 + (5*sin(6*t))/192 + (5*sin(2*t))/32)*k1 + sin(4*t)/8 + sin(8*t)/64 + (3*t)/8)*k2^4 + (((35*sin(2*t))/1024 + (7*sin(6*t))/1024 + (7*sin(10*t))/5120 + sin(14*t)/7168)*k1^2 + (-(15*sin(4*t))/256 - (3*sin(8*t))/256 - sin(12*t)/768 - (5*t)/32)*k1 - (5*sin(2*t))/16 - (5*sin(6*t))/96 - sin(10*t)/160)*k2^5 + (((7*sin(10*t))/1280 + sin(14*t)/1792 + (7*sin(6*t))/256 + (35*sin(2*t))/256)*k1 + sin(12*t)/384 + (3*sin(8*t))/128 + (15*sin(4*t))/128 + (5*t)/16)*k2^6)

对于a = 3, b = 2, c = 5, k1 = -125/397, k2 = 5/63, m1 = -125/2, m2 = 5/2, n1 = 397/2, n2 = 63/2

代入得到：

F11(t1)=-0.05337798447*sin(2.*t1) + 0.0003465662833*sin(4.*t1) - 1.683699352*10^(-8)*sin(12.*t1) - 1.470725462*10^(-6)*sin(8.*t1) - 1.946628469*10^(-7)*sin(10.*t1) - 0.00004744759329*sin(6.*t1) - 1.607758395*10^(-13)*sin(14.*t1) + 0.4486677593*t1
作差F11(t1)-F2(t1)函数得到(结果很完美!)

      
更进一步，我们 设a=10,b=10*sqrt(1-(k/10)^2),c=1~10,20~100 共计171个 样本函数（椭圆离心率e=sqrt(a^2-b^2)/a=0.1~0.9，锥体母线斜率c/a=0.1~10)

作差F11(t1)-F2(t1)函数得到：（计算表明仅当e>0.9时，且t->PI/2时误差较大，其余几乎控制在10^(-4)内，要想计算更准确，需要更多的项，但不实用）





（2）为了更精准的计算，我们保留分式（不做展开），只将根式做渐近展开，代入积分可以得到

F12(t)=p0*((((21*sin(t)*cos(t)^7)/160 - (77*sin(t)*cos(t)^5)/640 + (7*sin(t)*cos(t)^3)/128 - (21*sin(t)*cos(t)^9)/320 - (21*sin(t)*cos(t))/1024)/k2 + ((63*t)/8192 + (21*sin(t)*cos(t)^7)/512 - (63*sin(t)*cos(t)^5)/1024 + (189*sin(t)*cos(t)^3)/4096 - (105*sin(t)*cos(t))/8192)/k2^2 + (-(21*sin(t)*cos(t))/1024 - (7*sin(t)*cos(t)^5)/256 + (7*sin(t)*cos(t)^3)/256)/k2^3 + ((21*t)/2048 + (21*sin(t)*cos(t)^3)/1024 - (21*sin(t)*cos(t))/2048)/k2^4 - (21*sin(t)*cos(t))/(1024*k2^5) + (-(21*k0)/1024 + (21*t)/1024)/k2^6)*k1^6 + (((21*t)/2048 + (7*sin(t)*cos(t)^7)/128 - (21*sin(t)*cos(t)^5)/256 - (35*sin(t)*cos(t))/2048 + (63*sin(t)*cos(t)^3)/1024)/k2 + (-(7*sin(t)*cos(t)^5)/192 + (7*sin(t)*cos(t)^3)/192 - (7*sin(t)*cos(t))/256)/k2^2 + ((7*t)/512 - (7*sin(t)*cos(t))/512 + (7*sin(t)*cos(t)^3)/256)/k2^3 - (7*sin(t)*cos(t))/(256*k2^4) + ((7*t)/256 - (7*k0)/256)/k2^5)*k1^5 + ((-(5*sin(t)*cos(t)^5)/96 - (5*sin(t)*cos(t))/128 + (5*sin(t)*cos(t)^3)/96)/k2 + ((5*t)/256 + (5*sin(t)*cos(t)^3)/128 - (5*sin(t)*cos(t))/256)/k2^2 - (5*sin(t)*cos(t))/(128*k2^3) + ((5*t)/128 - (5*k0)/128)/k2^4)*k1^4 + ((t/32 - sin(t)*cos(t)/32 + sin(t)*cos(t)^3/16)/k2 - sin(t)*cos(t)/(16*k2^2) + (t/16 - k0/16)/k2^3)*k1^3 + (-sin(t)*cos(t)/(8*k2) + (t/8 - k0/8)/k2^2)*k1^2 + (-k0/2 + t/2)*k1/k2 + k0)
      
其中\(t0=k0=\frac{\arctan(\frac{\tan(t)\sqrt{b^2+c^2}}{\sqrt{a^2+c^2}}(a^2+b^2+2c^2)}{2\sqrt{(b^2+c^2)(a^2+c^2)}\)

对于a = 3, b = 2, c = 5, k1 = -125/397, k2 = 5/63, m1 = -125/2, m2 = 5/2, n1 = 397/2, n2 = 63/2, p0 = sqrt(794)/63, t0 = (63*arctan(tan(t)*sqrt(29)*sqrt(34)/34)*sqrt(29)*sqrt(34))/1972

F12(t)=26.40188106*t - 25.95321331*arctan(0.9235481454*tan(t)) + 0.002604742358*sin(t)*cos(t)^7 - 0.02219743895*sin(t)*cos(t)^5 + 0.1869524419*sin(t)*cos(t)^3 - 2.256838506*sin(t)*cos(t) - 0.0003603517212*sin(t)*cos(t)^9

我们绘制F1(t1)-F12(t1)的误差函数得到(结果很完美!)：



更进一步，我们 设a=10,b=10*sqrt(1-(k/10)^2),c=1~10,20~100 共计171个 样本函数（椭圆离心率e=sqrt(a^2-b^2}/a=0.1~0.9，锥体母线斜率c/a=0.1~10)

作差F11(t)-F2(t)函数得到：（计算表明有部分曲线(直观看仅有8条曲线误差较大)超出了误差范围，需要做进一步的理论分析，要想计算更准确，需要更多的项，但不实用）

数学星空

找到 苏步青 编写的<实用微分几何引论>中有关斜圆锥体的展开角计算

数学星空

接楼上

Jack315

数学星空 发表于 2024-3-30 12:10
接楼上

推导下苏老前辈教材中的公式。

圆锥底面圆的参数方程：
\(x(\varphi)=R\cos\varphi\)
\(y(\varphi)=l+R\sin\varphi\)
\(z(\varphi)=a\)
其中 \(\varphi\) 是参变量，\(0\le\varphi\le2\pi\) 。

极径及其微分：
\(r(\varphi)=\sqrt{x^2(\varphi)+y^2(\varphi)+z^2(\varphi)}\)
\(r=\sqrt{l^2+a^2+R^2+2lR\sin\varphi}\)
\(dr=\frac{1}{r}lR\cos\varphi d\varphi\)

弧长相等：
\(Rd\varphi=\sqrt{(dr)^2+(rd\theta)^2}\)

极角微分：
\(d\theta=\frac{1}{r}\sqrt{(Rd\varphi)^2-(dr)^2}\)
\(=\frac{R}{r}\sqrt{1-(\frac{l}{r}cos\varphi)^2}d\varphi\)
\(=\frac{R}{r^2}\sqrt{a^2+R^2+2lR\sin\varphi+l^2\sin^2\varphi}d\varphi\)

极角积分公式：
\(\theta=\int_0^\varphi\frac{R\sqrt{a^2+R^2+2lR\sin\varphi+l^2\sin^2\varphi}}{l^2+a^2+R^2+2lR\sin\varphi}d\varphi\)

网友 Huxley 推导用的是求向量夹角的方法，是否有问题还真说不好。
余弦展开为级数时没有奇次幂项，而向量点积展开为级数时则有奇次幂项。
因而在理论上或存在误差。至于对结果有多大影响没深究。

用极径相等和弧长相等在理论上是严谨的。