【数值积分】
自己写插值、数值积分……并不真的需要完全自己动手写。
如果编程语言是 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\)
=, , , , , , , , , , , , , , , , , , ,
级数展开：
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
代入上面的相关参数，得到
=, , , , , , , , , , , , , , , , , , ,
这个级数展开在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-29 07:21 编辑
数学星空 发表于 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)}\)
两个式子对不上，不知道是哪里搞错了。
其它的项可能有误，但余弦项应该是没对上。
你的分子计算都不对，你可以再算一遍
\(\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 代码验证：
f := Sqrt[((a^2 Cos^2 + b^2 Sin^2)/(a^2 Cos^2 + b^2 Sin^2 + c^2) - (((a^2 - b^2) Sin Cos)/(a^2 Cos^2 + b^2 Sin^2 + c^2))^2)]
h := Sqrt 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]
f - h /. x -> Cos // 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条曲线误差较大)超出了误差范围，需要做进一步的理论分析，要想计算更准确，需要更多的项，但不实用）
找到 苏步青 编写的<实用微分几何引论>中有关斜圆锥体的展开角计算
接楼上
数学星空 发表于 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 推导用的是求向量夹角的方法，是否有问题还真说不好。
余弦展开为级数时没有奇次幂项，而向量点积展开为级数时则有奇次幂项。
因而在理论上或存在误差。至于对结果有多大影响没深究。
用极径相等和弧长相等在理论上是严谨的。
重新推导椭圆锥体的公式。
椭圆锥底面椭圆的参数方程：
\(x(\varphi)=a\cos\varphi\)
\(y(\varphi)=b\sin\varphi\)
\(z(\varphi)=c\)
其中 \(\varphi\) 是参变量，\(0\le\varphi\le2\pi\) 。
极径及其微分：
\(r(\varphi)=\sqrt{x^2(\varphi)+y^2(\varphi)+z^2(\varphi)}\)
\(r=\sqrt{a^2\cos^2\varphi+b^2\sin^2\varphi+c^2}\)
\(dr=\frac{1}{r}(-a^2+b^2)\sin\varphi\cos\varphi d\varphi\)
弧长相等：
\(ds=\sqrt{a^2\cos^2\varphi+b^2\sin^2\varphi}d\varphi=\sqrt{(dr)^2+(rd\theta)^2}\)
极角微分：
\(d\theta=\frac{1}{r}\sqrt{(ds)^2-(dr)^2}\)
\(=\frac{1}{r}\sqrt{(a^2\cos^2\varphi+b^2\sin^2\varphi)-(\frac{-a^2+b^2}{r}\sin\varphi\cos\varphi)^2}d\varphi\)
\(=\sqrt{\frac{a^2\cos^2\varphi+b^2\sin^2\varphi}{a^2\cos^2\varphi+b^2\sin^2\varphi+c^2}-[\frac{(a^2-b^2)\sin\varphi\cos\varphi}{a^2\cos^2\varphi+b^2\sin^2\varphi+c^2}]^2}d\varphi\)
极角积分公式：
\(\theta=\int_0^\varphi\sqrt{\frac{a^2\cos^2\varphi+b^2\sin^2\varphi}{a^2\cos^2\varphi+b^2\sin^2\varphi+c^2}-[\frac{(a^2-b^2)\sin\varphi\cos\varphi}{a^2\cos^2\varphi+b^2\sin^2\varphi+c^2}]^2}d\varphi\)
与之前的结果一致，表明网友 Huxley 推导出来的公式至少在这个题目上是对的。