作业捷径！ scipy库实现曲率、挠率、Frenet标架计算！

---

本周作业算的曲率、挠率、Frenet标架实在算的让人秃头，废话不多说，我们直接上代码！

#曲率、挠率、Frenet标架计算
from sympy import *
def OuterProduct(A,B):
    #计算向量外积
    C=Matrix([A[1]*B[2]-A[2]*B[1],A[2]*B[0]-A[0]*B[2],A[0]*B[1]-A[1]*B[0]])
    return C
def MixedProduct(A,B,C):
    #计算混合积
    D=Matrix([[A[0],A[1],A[2]],[B[0],B[1],B[2]],[C[0],C[1],C[2]]])
    return D.det()
def magnitude(A):
    #计算向量模长
    return sqrt(A[0]*A[0]+A[1]*A[1]+A[2]*A[2])

def curvature(A):
    #计算曲率
    son=magnitude(OuterProduct(diff(A,t,1),diff(A,t,2)))
    mom=magnitude(diff(A,t))**3
    return son/mom
def torsion(A):
    #计算挠率
    return MixedProduct(diff(A,t,1),diff(A,t,2),diff(A,t,3))/magnitude(OuterProduct(diff(A,t,1),diff(A,t,2)))**2
def Frenet(A):
    #计算Frenet标架
    T=diff(A,t)/magnitude(diff(A,t))
    N=(1/curvature(A))*diff(T,t)/magnitude(diff(A,t))
    B=OuterProduct(T,N)
    return T,N,B

这是全部的函数，包含有计算外积、混合积、曲率、挠率、模长、Frenet标架的代码。我们下面对一个题目进行计算：


我们直接把相应函数输入：

#函数曲线
t=symbols('t')
C=Matrix([exp(-1/t**2),t,0])
factor(Frenet(C))

得到相应的Frenet标架：
$$\left(\left[\begin{array}{c}\frac{2e^{-\frac{1}{t^2}}}{t^3\sqrt{1+\frac{4e^{-\frac{2}{t^2}}}{t^6}}}\\ \frac{1}{\sqrt{1+\frac{4e^{-\frac{2}{t^2}}}{t^6}}}\\ 0\end{array}\right],\left[\begin{array}{c}\frac{\left(1+\frac{4e^{-\frac{2}{t^2}}}{t^6}\right)\left(\frac{2\left(\frac{12e^{-\frac{2}{t^2}}}{t^7}-\frac{8e^{-\frac{2}{t^2}}}{t^9}\right)e^{-\frac{1}{t^2}}}{t^3{\left(1+\frac{4e^{-\frac{2}{t^2}}}{t^6}\right)}^{\frac{3}{2}}}-\frac{6e^{-\frac{1}{t^2}}}{t^4\sqrt{1+\frac{4e^{-\frac{2}{t^2}}}{t^6}}}+\frac{4e^{-\frac{1}{t^2}}}{t^6\sqrt{1+\frac{4e^{-\frac{2}{t^2}}}{t^6}}}\right)}{2\sqrt{\frac{{\left(-3+\frac{2}{t^2}\right)}^2{e}^{-\frac{2}{t^2}}}{t^8}}}\\ \frac{\frac{12e^{-\frac{2}{t^2}}}{t^7}-\frac{8e^{-\frac{2}{t^2}}}{t^9}}{2\sqrt{\frac{{\left(-3+\frac{2}{t^2}\right)}^2{e}^{-\frac{2}{t^2}}}{t^8}}\sqrt{1+\frac{4e^{-\frac{2}{t^2}}}{t^6}}}\\ 0\end{array}\right],\left[\begin{array}{c}0\\ 0\\ -\frac{\sqrt{1+\frac{4e^{-\frac{2}{t^2}}}{t^6}}\left(\frac{2\left(\frac{12e^{-\frac{2}{t^2}}}{t^7}-\frac{8e^{-\frac{2}{t^2}}}{t^9}\right)e^{-\frac{1}{t^2}}}{t^3{\left(1+\frac{4e^{-\frac{2}{t^2}}}{t^6}\right)}^{\frac{3}{2}}}-\frac{6e^{-\frac{1}{t^2}}}{t^4\sqrt{1+\frac{4e^{-\frac{2}{t^2}}}{t^6}}}+\frac{4e^{-\frac{1}{t^2}}}{t^6\sqrt{1+\frac{4e^{-\frac{2}{t^2}}}{t^6}}}\right)}{2\sqrt{\frac{{\left(-3+\frac{2}{t^2}\right)}^2{e}^{-\frac{2}{t^2}}}{t^8}}}+\frac{\left(\frac{12e^{-\frac{2}{t^2}}}{t^7}-\frac{8e^{-\frac{2}{t^2}}}{t^9}\right)e^{-\frac{1}{t^2}}}{t^3\sqrt{\frac{{\left(-3+\frac{2}{t^2}\right)}^2{e}^{-\frac{2}{t^2}}}{t^8}}\left(1+\frac{4e^{-\frac{2}{t^2}}}{t^6}\right)}\end{array}\right]\right)$$
这确实不像一个让人算的东西，我们来看看t从正方向趋于0时的极限：

[[limit(Frenet(C)[0][0],t,0,dir='+'),limit(Frenet(C)[0][1],t,0,dir='+'),limit(Frenet(C)[0][2],t,0,dir='+')],[limit(Frenet(C)[1][0],t,0,dir='+'),limit(Frenet(C)[1][1],t,0,dir='+'),limit(Frenet(C)[1][2],t,0,dir='+')],[limit(Frenet(C)[2][0],t,0,dir='+'),limit(Frenet(C)[2][1],t,0,dir='+'),limit(Frenet(C)[2][2],t,0,dir='+')]]

$$[[0,1,0],[1,0,0],[0,0,-1]]$$
这里就友好的多了，相同的方式，我们对t小于0部分梅开二度：

D=Matrix([0,t,exp(-1/t**2)])
factor(Frenet(D))

$$\left(\left[\begin{array}{c}0\\ \frac{1}{\sqrt{1+\frac{4e^{-\frac{2}{t^2}}}{t^6}}}\\ \frac{2e^{-\frac{1}{t^2}}}{t^3\sqrt{1+\frac{4e^{-\frac{2}{t^2}}}{t^6}}}\end{array}\right],\left[\begin{array}{c}0\\ \frac{\frac{12e^{-\frac{2}{t^2}}}{t^7}-\frac{8e^{-\frac{2}{t^2}}}{t^9}}{2\sqrt{\frac{{\left(-3+\frac{2}{t^2}\right)}^2{e}^{-\frac{2}{t^2}}}{t^8}}\sqrt{1+\frac{4e^{-\frac{2}{t^2}}}{t^6}}}\\ \frac{\left(1+\frac{4e^{-\frac{2}{t^2}}}{t^6}\right)\left(\frac{2\left(\frac{12e^{-\frac{2}{t^2}}}{t^7}-\frac{8e^{-\frac{2}{t^2}}}{t^9}\right)e^{-\frac{1}{t^2}}}{t^3{\left(1+\frac{4e^{-\frac{2}{t^2}}}{t^6}\right)}^{\frac{3}{2}}}-\frac{6e^{-\frac{1}{t^2}}}{t^4\sqrt{1+\frac{4e^{-\frac{2}{t^2}}}{t^6}}}+\frac{4e^{-\frac{1}{t^2}}}{t^6\sqrt{1+\frac{4e^{-\frac{2}{t^2}}}{t^6}}}\right)}{2\sqrt{\frac{{\left(-3+\frac{2}{t^2}\right)}^2{e}^{-\frac{2}{t^2}}}{t^8}}}\end{array}\right],\left[\begin{array}{c}\frac{\sqrt{1+\frac{4e^{-\frac{2}{t^2}}}{t^6}}\left(\frac{2\left(\frac{12e^{-\frac{2}{t^2}}}{t^7}-\frac{8e^{-\frac{2}{t^2}}}{t^9}\right)e^{-\frac{1}{t^2}}}{t^3{\left(1+\frac{4e^{-\frac{2}{t^2}}}{t^6}\right)}^{\frac{3}{2}}}-\frac{6e^{-\frac{1}{t^2}}}{t^4\sqrt{1+\frac{4e^{-\frac{2}{t^2}}}{t^6}}}+\frac{4e^{-\frac{1}{t^2}}}{t^6\sqrt{1+\frac{4e^{-\frac{2}{t^2}}}{t^6}}}\right)}{2\sqrt{\frac{{\left(-3+\frac{2}{t^2}\right)}^2{e}^{-\frac{2}{t^2}}}{t^8}}}-\frac{\left(\frac{12e^{-\frac{2}{t^2}}}{t^7}-\frac{8e^{-\frac{2}{t^2}}}{t^9}\right)e^{-\frac{1}{t^2}}}{t^3\sqrt{\frac{{\left(-3+\frac{2}{t^2}\right)}^2{e}^{-\frac{2}{t^2}}}{t^8}}\left(1+\frac{4e^{-\frac{2}{t^2}}}{t^6}\right)}\\ 0\\ 0\end{array}\right]\right)$$

[[limit(Frenet(D)[0][0],t,0,dir='+'),limit(Frenet(D)[0][1],t,0,dir='+'),limit(Frenet(D)[0][2],t,0,dir='+')],[limit(Frenet(D)[1][0],t,0,dir='+'),limit(Frenet(D)[1][1],t,0,dir='+'),limit(Frenet(D)[1][2],t,0,dir='+')],[limit(Frenet(D)[2][0],t,0,dir='+'),limit(Frenet(D)[2][1],t,0,dir='+'),limit(Frenet(D)[2][2],t,0,dir='+')]]

$$[[0,1,0],[0,0,1],[1,0,0]]$$

这样我们发现：t从正向趋于0和从负向趋于0的对应Frenet标架是不连续的。

更重要的是：这东西终于不用自己算了！

---

（阅览量破1000，就算累趴，也出下集！）