ある関数から法線を求め，その傾き値（度表記）の求め方

Question

前提・実現したいこと

python 初心者です．
関数z(x)を微分して法線を求め，その法線の傾き値（度表記）を求めようとしています．
また，傾き値はその座標値(a,z(a))の値で変化すると思いますが，aが0.1間隔で0から6まで変化した分だけの情報（すばわち60個分の傾き値）が欲しいです．
下図がイメージです．

発生している問題・エラーメッセージ

インターネットなどでやり方を探しているのですが，初心者ということもあり，やり方がわかりません．

該当のソースコード

#傾き値の取得

import numpy as np
import matplotlib.pyplot as plt

# 描画範囲の指定
# x = np.arange(x軸の最小値, x軸の最大値, 刻み)
x = np.arange(0, 6.28, 0.1)

plt.grid(True)
plt.xlabel('$x$', fontsize=16)
plt.ylabel('$Z$', fontsize=16)
plt.title('upper die')

#縦横比1:1
plt.gca().set_aspect('equal', adjustable='box')

R=-19.41053
k=-59.85017
c=1/R

from sympy import Symbol, Derivative
#xを変数xに指定
x=Symbol('x') 

#関数入力
z=-(c*x**2/(1+(1-(1+k)*c**2*x**2)**(1/2))
#入力関数を指定変数で微分
Z=Derivative(z,x).doit()

# 横軸の変数。縦軸の変数。
plt.plot(x,Z)

# 描画実行
plt.show()

np.savetxt("points.csv", np.vstack([x, Z]).T, delimiter=",")

微分を行うところまでは，teratailユーザーのご協力もあり出来ましたが，ここから先，行き詰っています．
知識のある方，ご協力をお願い致します．

編集後

#法線の角度を求める

import matplotlib.pyplot as plt
import numpy as np
from sympy import *
from sympy import Derivative, Symbol, N



# xを変数xに指定
x = Symbol('x')
z = -(c*x**2/(1+(1-(1+k)*c**2*x**2)**(1/2))

prime_z = Derivative(z, x).doit()

a = Symbol('a')
z_a = z.subs(x, a)  # f(a)
prime_z_a = prime_z.subs(x, a)  # 点 (a, f(a)) の傾き
tangent = prime_z_a * (x - a) + z_a  # 点 (a, f(a)) の接線
normal = (-1 / prime_z_a) * (x - a) + z_a  # 点 (a, f(a)) の法線
theta = deg(atan(-1 / prime_z_a))  # 法線の角度
print('{} degree'.format(theta.evalf(subs={a: x, prime_z_a: x})))
 
# x, theta の点一覧を作成する。
xs = np.arange(0, 6.28, 0.1)
thetas = [theta.evalf(subs={x: v}) for v in xs]

# 描画する。
fig, ax = plt.subplots(figsize=(7, 7))
plt.grid(True)
plt.xlabel('$x$', fontsize=16)
plt.ylabel('$θ$', fontsize=16)
plt.title('upper die')
plt.plot(xs, thetas)
plt.show()

np.savetxt("points.csv", np.vstack([xs,thetas]).T, delimiter=",")

を入力すると，

>>> # 描画する。
... fig, ax = plt.subplots(figsize=(7, 7))
>>> plt.grid(True)
>>> plt.xlabel('$x$', fontsize=16)
Text(0.5,0,'$x$')
>>> plt.ylabel('$θ$', fontsize=16)
Text(0,0.5,'$\x83\xc6$')
>>> plt.title('upper die')
Text(0.5,1,'upper die')
>>> plt.plot(xs, thetas)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\Users\owner\Anaconda3\envs\python3.6\lib\site-packages\matplotlib\pyplot.py", line 3363, in plot
    ret = ax.plot(*args, **kwargs)
  File "C:\Users\owner\Anaconda3\envs\python3.6\lib\site-packages\matplotlib\__init__.py", line 1867, in inner
    return func(ax, *args, **kwargs)
  File "C:\Users\owner\Anaconda3\envs\python3.6\lib\site-packages\matplotlib\axes\_axes.py", line 1529, in plot
    self.add_line(line)
  File "C:\Users\owner\Anaconda3\envs\python3.6\lib\site-packages\matplotlib\axes\_base.py", line 1960, in add_line
    self._update_line_limits(line)
  File "C:\Users\owner\Anaconda3\envs\python3.6\lib\site-packages\matplotlib\axes\_base.py", line 1982, in _update_line_limits
    path = line.get_path()
  File "C:\Users\owner\Anaconda3\envs\python3.6\lib\site-packages\matplotlib\lines.py", line 956, in get_path
    self.recache()
  File "C:\Users\owner\Anaconda3\envs\python3.6\lib\site-packages\matplotlib\lines.py", line 657, in recache
    y = _to_unmasked_float_array(yconv).ravel()
  File "C:\Users\owner\Anaconda3\envs\python3.6\lib\site-packages\matplotlib\cbook\__init__.py", line 2052, in _to_unmasked_float_array
    return np.asarray(x, float)
  File "C:\Users\owner\Anaconda3\envs\python3.6\lib\site-packages\numpy\core\numeric.py", line 501, in asarray
    return array(a, dtype, copy=False, order=order)
  File "C:\Users\owner\Anaconda3\envs\python3.6\lib\site-packages\sympy\core\expr.py", line 225, in __float__
    raise TypeError("can't convert expression to float")
TypeError: can't convert expression to float
>>> plt.show()

この様なエラーが出てしまいます．

補足情報（FW/ツールのバージョンなど）

python3.6を使用しています．

Comments

[tiitoi]

数学である点の法線を求めたり、画像中の角度を求めるやり方はご存知ですか？
つまり、sympy の使い方の質問なのか、それとも数学の話なのか、わからないのはどちらでしょうか？

[mdshiba2]

数学的にはわかっております。プログラムの書き方(sympyの使い方になるのでしょうか)、がわからない状態ですm(__)m

Answer 1

点 (a, f(a) における法線の傾きは - 1 / f'(a) なので、

tanθ = - 1 / f'(a) の関係から角度を求めればよいかと思います。

法線の角度を求める。

import matplotlib.pyplot as plt
import numpy as np
from sympy import *

# xを変数xに指定
x = Symbol('x')
y = log(x)
prime_y = Derivative(y, x).doit()
a_val = 0.5  # a の値

a = Symbol('a')
y_a = y.subs(x, a)  # f(a)
prime_y_a = prime_y.subs(x, a)  # 点 (a, f(a)) の傾き
tangent = prime_y_a * (x - a) + y_a  # 点 (a, f(a)) の接線
normal = (-1 / prime_y_a) * (x - a) + y_a  # 点 (a, f(a)) の法線
theta = deg(atan(-1 / prime_y_a))  # 法線の角度
print('{} degree'.format(theta.evalf(subs={a: a_val})))
# -26.5650511770780 degree

結果を描画する。

# x, y の点一覧を作成する。
xs = np.linspace(0.001, 5, 100)
ys = [y.evalf(subs={x: t}) for t in xs]
tangent_ys = [tangent.evalf(subs={x: t, a: a_val}) for t in xs]
normal_ys = [normal.evalf(subs={x: t, a: a_val}) for t in xs]
# 描画する。
fig, ax = plt.subplots(figsize=(7, 7))
# 表示範囲の設定
ax.set_xlim(0, 5)
ax.set_ylim(-5, 5)
ax.axis('equal')
# xy軸の設定
ax.grid(True)
ax.set_xticks(np.arange(6))
ax.set_yticks(np.arange(-5, 6))
ax.set_xlabel('$x$', fontsize=16)
ax.set_ylabel('$y$', fontsize=16)
# 折れ線グラフの作成
ax.plot(xs, ys, label='f')
ax.plot(xs, tangent_ys, label='tangent')
ax.plot(xs, normal_ys, label='normal')
ax.legend()
plt.show()

追記

import matplotlib.pyplot as plt
import numpy as np
from sympy import *

init_printing()

x = Symbol('x')
c, k = 1, 3
z = -(c * x**2 / (1 + (1 - (1 + k) * c**2 * x**2)**0.5))
prime_z = Derivative(z, x).doit()

a = Symbol('a')
z_a = z.subs(x, a)  # f(a)
prime_z_a = prime_z.subs(x, a)  # f'(a)
tangent = prime_z_a * (x - a) + z_a  # 点 (a, f(a)) の接線
normal = (-1 / prime_z_a) * (x - a) + z_a  # 点 (a, f(a)) の法線
theta = deg(atan(-1 / prime_z_a))  # 法線の角度

xs = np.arange(0.1, 6.28, 0.1)
thetas = np.array([complex(theta.subs(a, t)) for t in xs])

# 保存
np.savetxt("points.csv", np.c_[xs, thetas.real, thetas.imag], delimiter=",")

# 読み込み
xs, real, imag = np.loadtxt("points.csv", delimiter=",", unpack=True)
thetas = real + 1j * imag

Comments

[mdshiba2]

ご返信遅くなってしまい申し訳ありません．ご回答ありがとうございます．
aの値が0から6まで0.1間隔で変化した時に，対応するθの値を取得したい場合はどのように書き直せばいいのでしょうか？ご教授して頂けると幸いです．

[mdshiba2]

訂正して頂いたプログラムを基に自分でもプログラムを書いてみましたが，エラーが出てしまいました．．．

[tiitoi]

追記に記載しました。
定数 c, k の値がわからないので適当に設定したら [0, 6] の範囲で複素数が出てきてしまうので、そこらへんは適当になおしてください。

[mdshiba2]

ありがとうございます！！大変助かりました．何度も返信してくださりありがとうございました．