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

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

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

根据曲线论基本定理一节课内容的知识我们知道了在曲率与挠率确定的条件下，空间曲线是唯一的（在刚体变换等价类的意义下）。那么，当我们知道空间曲线的曲率与挠率，我们理应可以按照以下常微分方程求出对应的曲线的参数表示：
{r′(s)=t(s)t′(s)=κ(s)∗n(s)n′(s)=−κ(s)∗t(s)+τ(s)∗b(s)b′(s)=−τ(s)∗n(s)\left\{ \begin{aligned} r'(s)=&t(s) \\ t'(s)=&\kappa(s)*n(s)\\ n'(s)=&-\kappa(s)*t(s)+\tau(s)*b(s)\\ b'(s)=&-\tau(s)*n(s)\end{aligned} \right. ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧r′(s)=t′(s)=n′(s)=b′(s)=t(s)κ(s)∗n(s)−κ(s)∗t(s)+τ(s)∗b(s)−τ(s)∗n(s)
用矩阵形式表示即为：
dds(r(s)t(s)n(s)b(s))=(010000κ(s)00κ(s)0τ(s)00−τ(s)0)(r(s)t(s)n(s)b(s))\begin{aligned} \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{aligned} dsd⎝⎜⎜⎛r(s)t(s)n(s)b(s)⎠⎟⎟⎞=⎝⎜⎜⎛000010κ(s)00κ(s)0−τ(s)00τ(s)0⎠⎟⎟⎞⎝⎜⎜⎛r(s)t(s)n(s)b(s)⎠⎟⎟⎞
对于非常系数线性微分方程，我们并没有较好的解决方法，但是我们可以通过计算机进行求解。

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

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

我们发现这里只需要求平面上的曲线，而不需要求空间上的曲线，因而我们只需要假定挠率τ=0\tau=0τ=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,br,t,n,br,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函数解出相应的方程组，并取出rrr的成分r1,r2,r3r1,r2,r3r1,r2,r3.
得到的结果为：
r1(s)=asinh(s);r2(s)=s2+1−1;r3(s)=0.r_1(s)=asinh(s);r_2(s)=\sqrt{s^2 + 1} - 1;r_3(s)=0.r1(s)=asinh(s);r2(s)=s2+1−1;r3(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))

得到了解的积分曲线：

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

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

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

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

1-2 基于MATLAB的空间曲线曲率挠率的数值计算
1-2 基于MATLAB的空间曲线曲率挠率的数值计算 1.工具向量函数:设曲线r(s)=(x(s),y(s))r(s)=(x(s),y(s))r(s)=(x(s),y(s))是一条正则曲线,其中ss ...
B样条数据点反求控制点绘制曲线（源码）
一.软件功能需求 1)所设计的软件应具有图形化用户界面(GUI): 2)用户在软件界面上可用随机数方式或手工方式输入若干曲线或曲面的数据点,例如起点.终点.列表型值点等,对于曲线,还可设置步长参数:对 ...
高数 | 【微分方程】已知常系数微分方程特解，反求原方程
方法一:待定系数法(要求事先知道方程是哪种形状) 方法二:相消法(要求只要知道是几阶就能做) 万能方法三步走 1.写通解 2.求导再求导(几阶方程求导几阶) 3.联立三个方程消掉C1,C2(从后两个 ...
计算机反求设计的一般步骤,逆向设计的概念和基本步骤
一.逆向设计的概念: 讲逆向设计前,先来看下传统产品开发的流程:构思-设计-产品原型.顾名思义,所谓逆向设计理念恰好与正向设计相反. 逆向设计,又反求设计,逆向工程,是一种基于逆向推理的设计,通过对现 ...
矢量叉乘能否反求矢量
这两天纠结一个问题,矢量叉乘能否反求矢量? a ⃗ × b ⃗ = c ⃗ \vec{a}\times\vec{b}=\vec{c} a ×b =c 如果已经知道矢量 a ⃗ \vec{a} a 和矢 ...
计算机反求设计的一般步骤,反求设计
(重定向自逆向设计) 反求设计(Inverse Design / Reversion Design) [编辑] 什么是反求设计反求设计是以先进的产品或技术为对象,通过深人分析,掌握其关键技术,在消化 ...
计算机反求设计的一般步骤,第七章反求工程概述.pptx
第七章反求工程概述与创新设计; 7.1反求工程概述;7.1.1反求设计反求工程首先是要进行反求分析,反求设计与传统的产品正向设计方法不同,他根据已存在的产品或零件来构造产品的工程设计模型或概念模型 ...
已知基础解系反求有效方程(矩阵)
已知基础解系反求有效方程(矩阵) @(数学) 这个是很有趣的推导过程,原理需要弄清楚. 即:已知Ax = 0的基础解系,由Ax = 0的系数行向量与解向量的关系可以反过来求解A. 具体推导如下: 齐次 ...
三、C++反作弊对抗实战 (实战篇 —— 3.如何获取游戏中角色人物角色的名称坐标、血量、武器信息(非CE扫描))
提示:本章节将介绍如何获取CS1.6游戏中角色人物角色的名称.坐标.血量.武器信息(非CE扫描前言在上一章节中<三.C++反作弊对抗实战 (实战篇 -- 2.认识CS1.6常见的数据结构与流 ...
Android 过反抓包总结入门篇
Android 过反抓包总结入门篇做协议分析少不了抓包,但是对于新入门我们来说,这是一大难题,网上各种各样的工具不少.但是如果遇到反抓包,那就芭比Q了.但是一些简单的反抓包机制还是很好过的. 1.要 ...

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

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

最新文章

热门文章

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

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

最新文章

热门文章

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

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