用曲率，挠率反求曲线方程！（作业捷径篇 续集）

---

时隔两天，前一篇文章阅览量竟然还未破1000（555），但是为了能够最快速度完成本周微分几何作业，我加工1h整出来了本次文章的计算代码。

上次文章我们解决的问题是给定曲线下曲率与挠率的计算，这次我们反其道而行之。

根据曲线论基本定理一节课内容的知识我们知道了在曲率与挠率确定的条件下，空间曲线是唯一的（在刚体变换等价类的意义下）。那么，当我们知道空间曲线的曲率与挠率，我们理应可以按照以下常微分方程求出对应的曲线的参数表示：
$$\{\begin{array}{rl}{r}^{\prime }(s)=& t(s)\\ t^{\prime }(s)=& \kappa (s)\ast n(s)\\ n^{\prime }(s)=& -\kappa (s)\ast t(s)+\tau (s)\ast b(s)\\ b^{\prime }(s)=& -\tau (s)\ast n(s)\end{array}$$
用矩阵形式表示即为：
$$\begin{array}{r}\frac{d}{ds}\left(\begin{array}{c}r(s)\\ t(s)\\ n(s)\\ b(s)\end{array}\right)=\left(\begin{array}{cccc}0& 1& 0& 0\\ 0& 0& \kappa (s)& 0\\ 0& \kappa (s)& 0& \tau (s)\\ 0& 0& -\tau (s)& 0\end{array}\right)\left(\begin{array}{c}r(s)\\ t(s)\\ n(s)\\ b(s)\end{array}\right)\end{array}$$
对于非常系数线性微分方程，我们并没有较好的解决方法，但是我们可以通过计算机进行求解。

我们举下面一个习题作为例子：

---

13. (1)求曲率$\kappa (s)=a/(a^2+s^2)$($s$为弧长参数)的平面曲线;
    (2)求曲率$\kappa (s)=1/\sqrt{a^2-s^2}$($s$为弧长参数)的平面曲线.

---

我们发现这里只需要求平面上的曲线，而不需要求空间上的曲线，因而我们只需要假定挠率$\tau =0$(挠率为0时，曲线为平面曲线)。

我们直接上代码：

from sympy import *
from sympy.plotting import plot3d_parametric_line,plot_parametric
import re

我们这次要用的库有上次使用过的符号计算库sympy,以及其对应的绘制三维和二维图像的两个函数，以及进行数据清洗的re库。

#自变量
s = symbols('s', real=True)
# r的三个分量
r1 = Function('r1')
r2 = Function('r2')
r3 = Function('r3')
#r=Matrix([r1,r2,r3])
# t的三个分量
t1 = Function('t1')
t2 = Function('t2')
t3 = Function('t3')
#t=Matrix([t1,t2,t3])
# n的三个分量
n1 = Function('n1')
n2 = Function('n2')
n3 = Function('n3')
#n=Matrix([n1,n2,n3])
# b的三个分量
b1 = Function('b1')
b2 = Function('b2')
b3 = Function('b3')
#b=Matrix([b1,b2,b3])

我们用sympy定义出相应的自变量和函数。需要注意的是，我们上面列出的微分方程中$r,t,n,b$都是三维的列向量，因而总方程数要有更多，在这里我们也是用分量的形式来进行定义的。

# 曲率与挠率
a=1
k = a/(a**2+s**2)# 曲率
tao = 0  # 挠率

我们输入题目中给定的曲率与挠率。(因为不会带参数的微分方程求解，因而暂时取参数a为1)（会的朋友欢迎加我好友进行交流）

eq = [r1(s).diff(s) - t1(s),
      r2(s).diff(s) - t2(s),
      r3(s).diff(s) - t3(s),
      t1(s).diff(s) - k * n1(s),
      t2(s).diff(s) - k * n2(s),
      t3(s).diff(s) - k * n3(s),
      n1(s).diff(s) + k * t1(s) - tao * b1(s),
      n2(s).diff(s) + k * t2(s) - tao * b2(s),
      n3(s).diff(s) + k * t3(s) - tao * b3(s),
      b1(s).diff(s) + tao * n1(s),
      b2(s).diff(s) + tao * n2(s),
      b3(s).diff(s) + tao * n3(s)]
con = {r1(0): 0, r2(0): 0, r3(0): 0, 
	   t1(0): 1, t2(0): 0, t3(0): 0, 
	   n1(0): 0, n2(0): 1, n3(0): 0, 
	   b1(0): 0, b2(0): 0, b3(0): 1}

我们将微分方程组和初值条件分别用列表和字典属性进行储存，

a = dsolve(eq, ics=con)
b=a[0:3]
print(str(b[0]),str(b[1]),str(b[2]),sep='\n')

用dsolve函数解出相应的方程组，并取出$r$的成分$r1,r2,r3$.
得到的结果为：
$$r_1(s)=asinh(s);r_2(s)=\sqrt{s^2+1}-1;r_3(s)=0.$$

我们还能对其曲线进行绘制：

R1=eval(re.findall(',.*\)',str(b[0]))[0][1:-1])
R2=eval(re.findall(',.*\)',str(b[1]))[0][1:-1])
R3=eval(re.findall(',.*\)',str(b[2]))[0][1:-1])
if R1==0:
    plot_parametric(R2,R3,(s,-100,100),title='R2-R3')
elif R2==0:
    plot_parametric(R1, R3, (s, -100, 100),title='R1-R3')
elif R3==0:
    plot_parametric(R1, R2, (s, -100, 100),title='R1-R2')
else:
    plot3d_parametric_line(R1,R2,R3,(s,-100,100))

得到了解的积分曲线：

同样的方法我们可以解答第二问的答案：
$$r_1(s)=s\ast \sqrt{1-s^2}/2+asin(s)/2;r_2(s)=s^2/2;r_3(s)=0.$$
其积分曲线的图像为：

---

原创不易，欢迎大家一键三连！

---

文案：锦帆远航
代码：锦帆远航
图片：锦帆远航 百度百科