断面二次モーメント
やろうとしている事
断面二次モーメントを数値計算的に算出しようとしている。（先々月の続き。）
計算の仕方は二つ
長さ寸法と、公式に当てはめて計算をする。
公式の公式集。書籍よりDVD版の方が明らかにお得だが高い。
- メリット
- 正確な値
- ディメリット
- 形状によっては式が用意されていない。
座標値を使い、数値計算的に求める
- メリット
- 断面形状にかかわらず、座標値さえあれば計算可能。
- ディメリット
- 点の間隔、カーブ、数値積分の手法などで誤差が生じる。
- 座標値を準備しなければならない。
1の方法が普通のやり方。できるなら1をやるべき。
断面二次モーメントの公式
$$\int _{\rm A}y^2\rm\ dA$$
上は、断面二次モーメントの式。
下は、8×8の正方形の断面二次モーメントを計算する際の積分のイメージ。
数値計算で断面二次モーメントを求めるための、計算手順
断面二次モーメントまでの道のり
座標値を手に入れる
- 断面一次モーメントを計算する。
- 断面一次モーメントから図心を計算する。
- 図心が原点となるように座標値を整える。(省略)
断面二次モーメントを計算する。
1. 座標値を手に入れる 四角形のx,y座標を得る関数
def coorsquare(r=4, center=[0,0], n=10, closed=True):
rx1=-r/2.0+center[0]
rx2=r/2.0+center[0]
ry1=-r/2.0+center[1]
ry2=r/2.0+center[1]
x = [rx1 + i*(rx2-rx1)/float(n) for i in range(n)]
y = [ry2 for i in range(n)]
#1
x = x + [rx2 for i in range(n)]
y = y + [ry2 - i*(ry2-ry1)/float(n) for i in range(n)]
#2
x = x + [rx2 - i*(rx2-rx1)/float(n) for i in range(n)]
y = y + [ry1 for i in range(n)]
#3
x = x + [rx1 for i in range(n)]
y = y + [ry1 + i*(ry2-ry1)/float(n) for i in range(n)]
#4
if closed:
x = x + [rx1]
y = y + [ry2]
return x, y
x,y = coorsquare(r=5, center=[0,0], n=8, closed=True)上の関数を使って作ったx,yを作図。
import matplotlib.pyplot as plt
fig,ax1 = plt.subplots(1)
ax1.plot(x,y)
ax1.scatter(x,y)
for i in range(len(x)):
ax1.text(x[i],y[i],"%s"%(i))
ax1.set_xlim(-4,4)
ax1.set_ylim(-4,4)
ax1.set_aspect("equal")
plt.show()
2. 面積を計算→断面一次モーメント→図心を計算。
def area(x, y):
A = sum([(y[i]+y[i+1])/2.0*(x[i+1]-x[i]) for i in range(len(x)-1)])
return A
x,y = coorsquare(r=5, center=[0,0], n=20, closed=True)
print(area(x,y))
def FirstInrt(x, y):
FIy = sum([(y[i]+y[i+1])**2/8.0*(x[i+1]-x[i]) for i in range(len(x)-1)])
FIx = sum([(x[i]+x[i+1])**2/8.0*(y[i]-y[i+1]) for i in range(len(x)-1)])
return FIx, FIy
print(FirstInrt(x, y))
def CntCS(x, y):
A = area(x, y)
FIx, FIy = FirstInrt(x, y)
Cntx, Cnty = FIx/A, FIy/A
return Cntx, Cnty
print("center(%s,%s)"%(CntCS(x, y)))3. 断面二次モーメントの計算（シンプソン則）
二次の式なので、台形則でなくシンプソン則で計算した。この関数が、この文章の一番重要な所
def simpson(x,y):#sekibun_1/3simpson(3point)
Cntx, Cnty = CntCS(x, y)
x = np.array(x)-Cntx
y = np.array(y)-Cnty
SecondInrty,SecondInrtx = 0,0
x1 = np.array(x)**3/3.0
y1 = np.array(y)**3/3.0
for i in range(len(x))[::2][:-1]:
a1 = float(x[i+1]-x[i])
# if は、0割対策。
if a1 !=0.0:
a1 = y1[i] * (x[i+2]-x[i])*(-x[i+2]-2.0*x[i]+3.0*x[i+1])/a1/6.0
a2 = float((x[i+1]-x[i])*(x[i+2]-x[i+1]))
if a2 !=0.0:
a2 = y1[i+1] * (x[i+2]-x[i])**3/a2/6.0
a3 = float(x[i+2]-x[i+1])
if a3 !=0.0:
a3 = y1[i+2] * (x[i+2]-x[i])*(2.0*x[i+2]+x[i]-3.0*x[i+1])/a3/6.0
SecondInrty = SecondInrty + a1 + a2 + a3
b1 = y[i+1]-y[i]
if b1 !=0:
b1 = x1[i]*(y[i+2]-y[i])*(-y[i+2]-2*y[i]+3*y[i+1])/b1/6.0
b2 = (y[i+1]-y[i])*(y[i+2]-y[i+1])
if b2 !=0:
b2 = x1[i+1]*(y[i+2]-y[i])**3/6.0/b2
b3 = y[i+2]-y[i+1]
if b3 !=0:
b3 = x1[i+2]*(y[i+2]-y[i])*(2*y[i+2]+y[i]-3*y[i+1])/b3/6.0
SecondInrtx = SecondInrtx - (b1 + b2 + b3)
return SecondInrtx, SecondInrty
# 断面が中心にない場合でも、図心を原点に持って行った時の値になるようにした。
x,y = coorsquare(r=5, center=[7,7], n=8, closed=True)
print("1/3simpson(3point)%s,%s"(simpson(x,y)))
# nは、偶数にしないと、正しく計算されない。
#うまくいかない。
x,y = coorsquare(r=5, center=[0,0], n=7, closed=True)
print("1/3simpson(3point)%s,%s"%(simpson(x,y)))
# 全体で偶数になっていても、一辺が奇数ではだめ。
#うまくいかない
x,y = coorsquare(r=5, center=[0,0], n=8, closed=True)
x.insert(-9,(-2.5-1.875)/2.0)
y.insert(-9,-2.5)
y.insert(-1,(2.5+1.875)/2.0)
x.insert(-1,-2.5)
print("1/3simpson(3point)%s,%s"%(simpson(x,y)))
# nが偶数なら、不等間隔でも計算される。
x,y = coorsquare(r=5, center=[0,0], n=8, closed=True)
x.insert(-9,(-2.5-1.875)/2.0)
y.insert(-9,-2.5)
x.insert(-9,-2.4)
y.insert(-9,-2.5)
print("1/3simpson(3point)%s,%s"%(simpson(x,y)))for ループは、2個飛ばしで行う。
3点使って三個飛ばしだと、間が空く。
3点使って二個飛ばしだと、間が埋まる。
ごめん、うまく書けない。こういう事。
ってなっちゃうから、3個飛ばしにすればいいわけじゃない。
3個使うけど、最後の一個は次のループでも使うから、二個飛ばしが正解。
(ちなみに、pythonで二個飛ばしの書き方は、[::2]。)
こういうindexingの所で、いつも試行錯誤してる。頭だけで出来ないとだめだろうか？
前のループでi+2だったところを、次のループではiとして扱っている。
積分は不等間隔のシンプソン則で模擬。不等間隔のシンプソン則は、ココを参考にして作った。
一辺のサンプルが奇数ならば正しく計算できる。
ひし形(45度回転。一部不等間隔。)
x,y = coorsquare(r=5, center=[0,0], n=8, closed=True)
x.insert(-9,(-2.5-1.875)/2.0)
y.insert(-9,-2.5)x.insert(-9,-2.4)
y.insert(-9,-2.5)
xy= np.dot(np.array([x,y]).T, [[np.cos(1/4.0*np.pi), np.sin(1/4.0*np.pi)],[-np.sin(1/4.0*np.pi), np.cos(1/4.0*np.pi)]]).T
x=xy[0]
y=xy[1]
print("1/3simpson(3point)%s,%s"%(simpson(x,y)))
前回は、これが正しく計算できなかったが、今回は不等間隔シンプソン則のおかげで計算できた。
（サンプル数の偶奇に気を付けなければならないが）
穴あき（正方形）
x1,y1 = coorsquare(r=5, center=[0,0], n=8, closed=True)
x2,y2 = coorsquare(r=3, center=[0,0], n=4, closed=True)
x2=x2[::-1]
y2=y2[::-1]
Ix1,Iy1 = simpson(x1,y1)
Ix2,Iy2 = simpson(x2,y2)
a = Ix1+Ix2
b = Iy1+Iy2
print("1/3simpson(3point)%s,%s"%(a,b))import matplotlib.pyplot as plt
fig,ax1 = plt.subplots(1)
ax1.plot(x1,y1)
ax1.plot(x2,y2)
ax1.scatter(x1,y1)
ax1.scatter(x2,y2)
for i in range(len(x1)):
ax1.text(x1[i],y1[i],"%s"%(i))
if i<len(x2):
ax1.text(x2[i],y2[i],"%s"%(i))
ax1.set_xlim(-3,3)
ax1.set_ylim(-3,3)
ax1.set_aspect("equal")
plt.show()
緑色の部分が反時計回りになっている所が重要。反時計回りはdx（差分）がマイナスになるので、負の面積というか負の断面二次モーメントが出てくる。
円形
N=21
x = np.cos(np.linspace(0, 2.0*np.pi, N))
y = np.sin(np.linspace(0, 2.0*np.pi, N))
print("1/3simpson(3point)%s,%s"%(simpson(x,y)))
円は、理論値からの誤差がまだある。どのくらいサンプル間隔を細かくすればいいのか考えるために、下で誤差などを確認した。
円の理論値
$$\frac{d^{4}\pi}{64} \\ 半径r=1なので直径d=2\\ \frac{16\pi}{64}=0.25\pi \\ =0.78539816339744830961566084581988 \\ (電卓でπ＝3.1415926535897932384626433832795を使い計算)$$
円周率にはnp.pi 3.141592653589793（16桁）を使うのでそれが原因の誤差が出るが、確認しない。
va = 0.25 * np.pi
Ix = []
Iy = []
Acx = []
Acy = []
Erx = []
Ery = []
for N in np.array(np.logspace(1,4,10),"int"):
#for N in np.array(np.logspace(1,4,10),"int")*4+1:
x = np.cos(np.linspace(0, 2.0*np.pi, N))
y = -np.sin(np.linspace(0, 2.0*np.pi, N))#Clockwise
I = simpson(x,y)
Ix.append(I[0])
Iy.append(I[1])
Acx.append(Ix[-1]/va)
Acy.append(Iy[-1]/va)
Erx.append(1 - Acx[-1])
Ery.append(1 - Acy[-1])
print("%s-points,\tIx:%s,\tIy:%s,\tAcx:%s,\tAcy:%s,\tErx:%s,\tEry:%s,\t"%(N,Ix[-1],Iy[-1],Acx[-1],Acy[-1],Erx[-1],Ery[-1]))import matplotlib.pyplot as plt
x = np.array(np.logspace(1,4,10),"int")#*4+1
fig,(ax1,ax2) = plt.subplots(2,1)
ax1.plot(x,[va for i in range(len(x))])
ax1.scatter(x, Acy)
ax2.plot(x,[1 for i in range(len(x))])
ax2.scatter(x, Ery)
ax1.set_ylim(0.99,1)
ax2.set_ylim(0.00000000000001, 0.1)
ax1.set_xscale("log")
ax2.set_xscale("log")
ax1.set_yscale("log")
ax2.set_yscale("log")
plt.show()
上が精度で１に近づけばいい値。理論値で割って計算した。y軸はなんか知らんがすっ飛んだ。0.99~1.0間。
下が誤差で、0に近づけばいい値。1-精度で計算。両方とも対数グラフ。
サンプル数によっては、うまく計算できないのがある。実行して、printの結果を確認すればわかるけど、10の時に断面二次モーメントに異常な値が出る。dxの±が入れ替わる所か、yの±が入れ替わる所で起こっているのではないかと思っているが、確認は次回以降にする。とにかく、サンプル数に気を遣わずにシンプソン則を実行できなければならない。あと、欲を言えば、excelで出来れば言うことない。有料ソフトなのに、どこに行っても必ずあるし。このままでは、ベジェ曲線で自由に描いた線の積分ができない。円形で起こった現象は、ひし形でも出てきそう。正方形で出てこなかったのは、dxがぴったし0だから。
実用上は、40個位の分割でも円の断面二次モーメントを10-4の誤差で計算できる事が知れてよかった。
ちなみに、データ点数を4n+1（nは整数）にすれば、異常な値は出てこない。
現在の課題
- ひし形、丸形などで、誤差が大きくなる。
- 任意形状（ひらがなの「あ」とか）でやってない。
次回以降の目標
- 【断面二次】 シンプソン則のサンプル数を気にせず実行。
【断面二次】Excelでやる。
- 一覧表を作成し、誤差等をチェック
【three.js】 pdbを表示する
絶対やるというわけではない。
