python实现隐函数曲率求解

一个大作业的实现，功能为对隐函数进行单位切向量以及曲率求解，在此分享出来，源码如下，测试可运行，希望对您有帮助：

from sympy import *

#声明函数与变量
c = symbols(‘c x y z’,cls=Function)
t = symbols(‘t’)

#定义隐函数集
xt = 2 * t
yt = t2 + 1
zt = t3 + 2 * t + 1

t0 = 1

#原函数
print(“xt原函数：”, xt)
print(“yt原函数：”, yt)
print(“zt原函数：”, zt)
#代入t = 1
xt_ = xt.evalf(subs={‘t’:t0})
yt_ = yt.evalf(subs={‘t’:t0})
zt_ = zt.evalf(subs={‘t’:t0})
print(‘代入t=’, t0)
print(“xt原函数求解：”, xt_)
print(“yt原函数求解：”, yt_)
print(“zt原函数求解：”, zt_, ‘\n’)

#求一阶偏导
#dxt = diff(xt, t, 2)
dxt = diff(xt, t)
dyt = diff(yt, t)
dzt = diff(zt, t)
print(“x对t的一阶偏导：”, dxt)
print(“y对t的一阶偏导：”, dyt)
print(“z对t的一阶偏导：”, dzt)

#代入t = 1
xt1 = dxt.evalf(subs={‘t’:t0})
yt1 = dyt.evalf(subs={‘t’:t0})
zt1 = dzt.evalf(subs={‘t’:t0})
print(‘代入t=’, t0)
print(“xt一阶偏导函数求解：”, xt1)
print(“yt一阶偏导函数求解：”, yt1)
print(“zt一阶偏导函数求解：”, zt1, ‘\n’)

#求二阶偏导
dxt2 = diff(xt, t, 2)
dyt2 = diff(yt, t, 2)
dzt2 = diff(zt, t, 2)
print(“x对t的二阶偏导：”, dxt2)
print(“y对t的二阶偏导：”, dyt2)
print(“z对t的二阶偏导：”, dzt2)

#代入t = 1
xt2 = dxt2.evalf(subs={‘t’:t0})
yt2 = dyt2.evalf(subs={‘t’:t0})
zt2 = dzt2.evalf(subs={‘t’:t0})
print(“xt二阶偏导函数求解：”, xt2)
print(“yt二阶偏导函数求解：”, yt2)
print(“zt二阶偏导函数求解：”, zt2)

#求解单位切向量
#C’(t) = [xt1,yt1,zt1]
normal = sqrt(xt1xt1 + yt1yt1 + zt1*zt1)
if(normal == 0.0):
print(‘单位切向量不存在。’)
else:
print(‘单位切向量为:’, xt1/normal, yt1/normal, zt1/normal)

#求解C(t)在t0出处的曲率
#C’‘(t) = [xt2,yt2,zt2]
#计算||C’(t)||2^3
normal_3 = normalnormalnormal
#计算C’(t)叉乘C’'(t)
aXb1 = yt1zt2 - yt2zt1
aXb2 = -(xt1zt2 - xt2zt1)
aXb3 = xt1yt2 - xt2yt1
normal_up = sqrt(aXb1aXb1 + aXb2aXb2 + aXb3*aXb3)
if(normal_3 == 0.0):
print(‘normal_3 == 0.0，此曲率公式无法用于本场景计算。’)
#其它方法：函数极限逼近

else:
kt = normal_up / normal_3
print(‘所求曲率为:’, kt)

运行结果如下：