工学野郎の思う微分_2 [ddt³さんの部屋]
工学野郎の思う微分_2
起動しなくなった旧PCより外付け装置に移植したHDですが、「フォーマットされていません」が出やがったので、データ復旧店に持ち込みました。全てではなかったけれど、概ねデータは救出できました。
・・・という訳で、ほぼ復旧記念です(^^)。
2.概念としての微分
工学野郎の思う微分_1の1.で、「概念上の話としては、微分可能な関数の一点には、理想化すればその接線が住んでると認めるのと同じです」と書きました。では、微分可能な関数の一点にその接線が住んでるというような、証拠はあるんでしょうか?。あるんですよ、それが。もちろん状況証拠ですけどね(^^;)。以後の喋りは、たぶんにイメージ優先という事になります。
関数の一点に接線がいる事を見るには、いったいどうしたら良いでしょう?。単純に考えれば、一点の内部を覗けば良いわけです。だけど点は大きさを持ちません(無限に小さい)。大きさを持たなければ内部構造もないはずだから、そもそも一点の内部なんてのは概念矛盾だと、プロゲラさんに怒られそうですが、でも[大きさを持たない]=[無限に小さい]が正しいなら(← 正しいという保証はない(^^;))、無限に拡大すりゃぁ~内部を覗けれんじゃねっ?、って発想になります。
そうすると拡大対象は関数なので、関数の拡大と縮小をまず復習しておきましょう。関数、
があった時、kを任意の実数として関数、
は、原点を中心に(1)を1/k倍に拡大した関数でした。どうしてかというと、
という変換を考えてみると、(3)を(2)に代入し、
が得られて、(4)は(1)と同じ曲線を描きます。しかし(X,Y)の中身は(kx',ky')なので、(4)が(1)と同じ曲線を描くためには、(x,y)=(kx',ky')でなければならない事になり、(x',y')の数値は、
と(x,y)の1/k倍になるためです。
高校では(x',y')と(x,y)をごっちゃにして扱う傾向があるために、高校生は(2)を見て非常に悩むんですよね(^^;)。
変換、
は等方的でもあります。それを確認するには、円で試せば十分です。円の方程式、
に変換(6)を施せば、
なので、(6)が原点から(x,y)までの距離を等方的に1/k倍にしてるのは明らかですし、原点は移動しないので、原点を中心とした拡大/縮小である事もわかります。等方性は、今の場合の関数の一点の拡大に望まれる性質です。
不便なのは原点中心でしか拡大/縮小できない事ですが、それなら拡大したい関数の一点を原点に移動させればOKです。このとき関数のグラフの形を変えてはいけないので、平行移動という話になります。
が、(1)をx方向にa,y方向にbだけ平行移動した関数を表すのでした。ここでも変換、
を考えれば、(x,y)→(x+a,y+b)である事は明らかです。
以後、高校方式で、(x,y)と(x',y')を明確に区別するのはやめます。
やりたい事は、関数の一点(x0,f(x0))が原点に来るように関数を平行移動し、そこで1/h倍(0<h<1)に拡大する事です。
最初に、(x0,f(x0))が原点に来るように関数を平行移動します。これは(a,b)=(-x0,-f(x0))の平行移動です。
次に(9)を原点中心で1/h倍に拡大します。
無限に拡大するには、h→0の極限をとれば良いわけですが、そのために(10)を少々変形します。
h≠0とすると、
(11)の最後の形でh→0の極限をとります。
最後の結果を出すために、(普通の(^^))微分の定義を使ってます。すなわち、
です。(12)は、(x0,f(x0))が原点に来るように平行移動した結果です。もとの位置に戻しましょう。それには、(a,b)=(x0,f(x0))なる平行移動をすれば良いので、
です。・・・うふふ、(x0,f(x0))における接線の方程式じゃ~ないですかぁ~(^^)。接線はいたんですよ。
ここで、微分を使わずに接線を定義するにはどうするかという(不要に?)頭の痛い問題が生じますが、それは無視します(^^;)。
次に、不要ではない哲学的批判に対処します。恐らく、ネムネコファミリーの中でこの批判を出すとしたら、プロゲラさんです(^^)。批判とは、
「点とは大きさを持たないこと。それが点の数学的定義だ」
「したがって、点は内部構造を持たない」
「よって内部構造を持たない点内部に、接線は存在できない」
「(12)で行った極限計算はその経過から明らかなように、点外部の拡大を行ったに過ぎず、点内部を覗いたとは、まったく言えない」
「これは点の数学的定義から最初から明らかで、点内部にその接線がいるという考えは、概念としてまさしく概念矛盾を起こしている。(12)~(14)は、概念的に無意味な数学的トリックに過ぎない」
はっきり言って論理的に完璧な批判です。反論は恐らく無理です(^^;)。ただ数学は、(12)~(14)を無意味な数学的トリックとは考えません。概念矛盾であると知りつつ、というか知ってるからこそ(本当によぉ~く考えた人たちは)、次のように身をかわします。
「現実の分析に便利なように、VR(バーチャルリィアリティー)なツールを作っただけなんですよね」
・・・と(^^;)。
(14)は(12)の洗礼の後に出てきたものなので、(14)の(x-x0)はいくら大きくとも無限倍拡大の縮尺で考えると、実スケールでは0です。だからそれは、VRな点内部の距離と考えてもいいじゃないかと(大きさを持たないことが点の数学的定義)。つまり微分可能な曲線は、各点ごとの接線の集合だと考えられる。
曲線の各点ごとの接線を近似するものはなんでしょう?。十分細かい区間分割で考えられた、曲線の折れ線近似ですよ。
よって、関数の一点の内部に接線がいるとVRとして認める事は、曲線は区分的な一次関数の集まりだと認めるのと同じです。これはさっき言ったように、哲学的な概念分析としては非常に危ういのかも知れませんが、認めてみると非常に役に立ちました。
したがって微分とは、一点における(局所的な)関数の一次関数化(線形化)、すなわち関数の局所線形化という人間の意志の産物なのです。
(執筆:ddt³さん)
[無限遠までの積分] by ddt³さん [ddt³さんの部屋]
[無限遠までの積分]
>それはそれとして、ネムネコが考えるに、数値計算屋のddt³さんが、きっと、これに関係する記事を投稿してくれると思う・・・
・・・と言われましても、いったい何を書けばいいんだろう?(^^;)。ネコ先生よりかっちょ良い変数変換でも行えばいいのだろうか?。いや、そうじゃない。「数値計算屋」のご指定があるからには(否定はしないけど(^^;))、ネコ先生のお望みはきっと台形公式のみを使った強行突破だ!。変数変換なんてもっての外の、とんでもはっぷん!・・・、なのかな?(^^)。
問題の被積分関数は、
でしたね(図-1)。
いやぁ~、これはまた性質の良さそうな関数です。
1から滑らかに0を収束する優美な姿(^^)。
でも問題は、(1)を0~∞で積分せよ、というものでした。
鬱陶しいのは、積分区間の上限が∞である事です。積分区間の上限が∞だと、どこで数値的足し算を打ち切って良いのかを判断しなきゃなりません。自分は有限幅の積分区間の場合、区間幅の100分割を標準にしてます。理由は、100等分して折れ線で結んでやると、大概の関数は十分に滑らかに描けるから。ところで折れ線近似の面積とは、台形公式による結果ですよね?。だから100分割で十分と・・・(^^)。
図-1をみると0~∞の積分の大勢は、0~10くらいで決まりそうです。なので区間分割幅δは、δ=10/100=0.1とします。図-1の白丸は、0.1刻みです。十分ですよね?。
次に積分の打ち切り範囲を考えます。図-1から明らかなように、x=10から先は、ほぼ傾き0の直線で近似できそうです。そこで10<xでのyの高さから、その点における接線を引き、x軸と交わる点x+Lを計算して、図-2に示す三角形の面積を計算します。
なので
から、
三角形の面積をδSとすれば、
このδSで、xから先の残りの面積を近似してやっても、そんなに間違いではないでしょう。δSが十分小さくなった時、計算を打ち切れば良い訳です。小ささの基準は次のように決めます。
xの分割幅をδとすると、個々の台形の面積Ajには、Aj×δ2程度の誤差が伴われると考えられます(←ネコ先生がそう言っていた(^^)←それは、(有界な区間の)大域的な誤差で、局所的な誤差はδ³程度ですよ(^^))。
全体では∑j Aj×δ2~F(x)・δ2です。ここでF(x)はxまで積分した時の正確な面積。これをF(x)=∑j Ajとして台形公式による面積で代用してやると、δ=0.1だったので、0.01・F(x)。ここにF(x)は、台形公式でxまで積分した結果。
要するにF(x)の1%くらいの誤差は許容できる、という事になりますので、残りの面積の近似δSについては、F(x)の0.1%くらいになれば十分じゃないでしょうか?。
つまり、F(x)=∑j AjとδS(x)をモニターしながら台形公式による積算を行い、δS(x)<F(x)/1000の時点で計算を打ち切る。
実際に0~10を台形公式で計算すると、F=1.471なので、この結果をもとに打ち切り工数を概算すると、
から、169.938<xで、169.938/0.1=1699~1700ステップくらいになります。
エッ?、1700ステップは多いって?。そもそも数値計算をやろうって奴が、1000や2000や10000でガタガタ言うんじゃありません。エッ?、計算はExcelでやりたいって?。だからぁ~、現在のExcelは1,048,576行まで使えるんだって・・・。
エッ?、10,000,000ステップだったら、どうするって?。100万行の列を10列つくりゃぁ~いいだろう!。現在のExcelのフル性能は、1,048,576行×16,384列=17,179,869,184ステップ分もあるんだ、どうだ参ったか?(← Microsoftの狗(^^;))。
なに~っ!、シートを全部ステップで埋めたら、結果を書く場所がないって?。隣のシートを利用せんかい!。もぉ~うるさいなぁ~、100万行も数式をマウスでコピーできないって?。数値計算をやるなら、Excelでマウスは使うな。[Cnt]+[Shift]+[上下左右カーソル]のショートカットキーをおぼえろぉ~!。
・・・もっとも50万行の列が1個あっただけで、重くて嫌になりますけどね(^^;)。で、実際にやってみると1598ステップ,x=159.8で、F=1.564539になります。|F-π/2|/(π/2)=0.004。Fに残りの予想値δS=0.001564を加えるとF=1.566103、|F+δS-π/2|/(π/2)=0.003です(^^)。収束状況を、図-3に示します。
さっき言ったように、現在の表計算ソフトは1000や2000のステップ数は問題にしませんから、このような素直な方法は、とてもよろしいと思うんですよね。ddt^3は時々ネコ先生に促されて、区間の端で∞になるような可積な関数、例えば(図-4)、
のような関数の0~1での数値積分を扱いますが、その時は、台形公式の増分精度をあらかじめ決めておいて、それに合わせてxが1に近づくほど、分割幅を短くするような方法をとりました(※)。
いやぁ~、ネコ先生が数値計算屋って指定するもんだから、ついつい強行突破をはかっちゃって・・・(^^;)。でも分割幅を可変にするのは、どちらかというと禁じ手です。分割幅が変わったために起きる種々の不便さの弊害なんかが考えられますし、分割幅を変える手続きが入るので素直さに欠けます。
それならいっその事(面積を求めるだけなら)、(5)を積分するんでなく、その逆関数の上側の面積を今回の方法で0~∞で計算するのもいいのかもしれない、と思っています。
(※) 参考記事
【面倒なことはしたくない】
https://nekodamashi-math.blog.ss-blog.jp/2018-09-10-1
この手の積分の近似計算は、逆に、積分区間を(半)無限区間に適当に変換し、たとえば、
という形にし、この半無限積分の近似値を求めるなんてことが行われたりするようです。
たとえば、
二重指数関数型積分公式
工学野郎の思う微分_1 【ddt³さん、半回復記念、特別寄稿記事】 [ddt³さんの部屋]
【ddt³さん、半回復記念特別寄稿】 工学野郎の思う微分_1
先日「自宅のPCが起動しなくなったぁ~!」と悲鳴をあげましたが、新しいPCを購入し(マザーボードがお釈迦という診断でした)、旧PCのHDは抜き出して外付け装置に移植しデータ救出。目出度しめでたし・・・のはずだったのですが、昨日外付けHDのデータを内蔵HDに移そうと思ってアクセスしたら、「フォーマットされていません」が出やがったんですよ(^^;)。「これは外付け装置の外周りの故障に違いない。中に移植したHDは無事に違いない!(そう思いたい!)」という状況になっちゃいました(^^;)。
・・・という訳で、半復活記念です(^^)。
1.二つの微分可能の定義
とりあえず普通の微分可能の定義をあげます。
[定義-1 微分可能]
関数f(x)に対して、
という計算が可能なとき、関数f(x)はxで微分可能という。関数f(x)がxで微分可能なら、(1)の結果をf'(x)と書く。すなわち、
ところが現代の(20世紀以降の)数学は、次の微分可能の定義を正式採用します。
[定義-2 微分可能]
hに依存しないa(x)が存在し、任意のhで、
が成り立つとき、関数f(x)はxで微分可能という。ただしO(h2)は、
を満たす(※)。
(※)
(3)、(4)ではなく、
関数fが
を満たすとき、関数fは点xで微分可能である、と定義するほうがメジャー。たとえば、高木貞治の「解析概論」。
また、ε−δ論法との兼ね合いで、
をもって、微分可能を定義する流儀もある。
英単語”where”は、日本語の「ここで」くらいの意味。
Ο(h²)、ο(h)に登場する、記号Οとοは、ともにLandau(らんだう)の記号と呼ばれるもの。
Ο(h²)、ο(h)の括弧()の中のh²、hにはそれぞれ意味があるのだけれど、これをやりだしたら、際限がなくなるので、ここでは触れない(^^ゞ。
これには、ちょっと、特殊な演算規則、色々なお約束事があるので。
「これでは曖昧すぎる」というヒトには、ε−δ論法による定義を。
任意の正数ε>0に対し、ある正数δ>0が存在し、0<|h|<δを満たす全ての実数hに関して、
が成り立つ、hに無関係な関数(?)aが存在するとき、関数fは点xで微分可能であるといい、
で表す、
といった感じになる。
(※終)
[定義-2]は、さすがに数学科で正式採用されただけあって不要に一般化されていて、一回読んだだけでは、いったい何を言ってんだかわかんないと思います(^^;)。しかし[定義-1]と[定義-2]は同等なのです。ここで同等とは、[定義-1]が成り立つなら必ず[定義-2]が成り立つし、[定義-2]が成り立てば必ず[定義-1]が成り立つ、という事です。[定義-1]と[定義-2]は、互いに互いの必要十分条件だ、という事です。
[定義-2]を採用するメリットこそ、ここで話したい事なのですが、まずは[定義-1]と[定義-2]は同等である事を証明します。それがないと話が進まないので。
[定理-1]
[定義-1]と[定義-2]は、同等な微分の定義である。
[証明]
1)[定義-2]⇒[定義-1]
[定義-2]が成り立つとすれば、(3)に次の変形が可能である。
h≠0のとき、
hは任意だったから、h→0でもかまわない。
上記に(4)を使えば、
となり、[定義-2]が成り立てばa(x)は存在するから、上記左辺の計算は可能である。これは[定義-1]。
2)[定義-1]⇒[定義-2]
[定義-1]が成り立つとすれば、(2)で定義されるf'(x)が存在する。そこで、
と定義する。上記に(2)を考慮すれば、
を導ける。さらにf'(x)はhに依存せず存在するので、それをa(x)と書いても良い。f'(x)をa(x)でおきかえ、(5)を(3)の形に移項すれば(6)((4))も成り立つ事から、これは[定義-2]である。
[証明終]
という訳で微分の定義としては、どっちでも良いんですよ。にも関わらず20世紀の数学科が、[定義-2]の表現を選んだのは、それがまず概念的,理論的に非常に有用だったからです。ここからは、その有用性について語るつもりです。
ちなみに「[定義-2]が不要に一般化されてる」というのは、[定理-1]の証明過程から明らかなように、[定義-2]は本来、h→0という状況で考えるのが目的なのは明瞭だからです。にも関わらず、任意のhでとなってしまうところが数学科(^^;)。
以後、hが0に近い事をh~0と書きます(hはほとんど0と読む)。[定理-1]の証明過程から(3)でa(x)をf'(x)と書いてもOKなのも明らかでしょう。
(10)は実質的には、h→0という状況を念頭においた式でした。h→0ならh~0であり、hの大きさはカスです。しかしここで皆さんに問いたいのは、大きい/小さいはどうやって判断してますか?、という話です。
h=0.1なら十分小さいと感じますけれど、それは普通h=1を基準にするからです。もしh=0.01を基準にしたら、h=0.1は基準より10倍も大きいです。何を言いたいかというと、大きい/小さいを言うためには、基準の大きさが必要ですよね?、という事です。基準の大きさでhを割って、その比率でhの大きさを判断してませんか?、という事です。
[定義-2]ではO(h2)をhで割ってました。という事は、O(h2)の大きさをhを基準に測ると言ってるのと同じです。ところが(4)が成り立つのでh~0でカスなhに対し、O(h2)は、h~0でさらにカスだという事になります。カス過ぎるものは無視したって、とりあえずOKなはずです。問題は人間が、どこまでの精度を求めるかでしょう。しかし最初は最低限の近似ですよね?。そして理屈の上では最低限の近似も、無限に精度を上げる事ができます。最低限の近似を考慮するだけで一般的には理論を組み立て可能です。よって理論屋な微分ユーザーの欲求からすれば、(10)のO(h2)は無視して良いのです!(^^)。
となります。(11)の意味はなんでしょう?。こうなります(^^)。
(D) 微小入力では、出力は入力に比例する.
・・・です。
hはh~0なので、hが微小入力です。その出力はf'(x)・hです。だってh=0ならf(x+h)=f(x)なんだから、という訳です。そしてf'(x)・hは、確かにhに比例します。比例定数は微分すればf'(x)と出るよ、というのが(11)の意味です。
(D)を認めるという事は概念上の話としては、微分可能な関数の一点には、理想化すればその接線が住んでると認めるのと同じです。これは、瞬間速度f'(t)を定義に従って計算する時に薄々ついてまわる気持ち悪さ、「瞬間では動いてないのに速度あり?」という疑問に対する概念上の答でもあります。
「動いてないけど接線はいるんだ。だから傾きはある。文句あっかぁ~!」
・・・という訳です(^^;)。
h~0なので(11)は、x近傍の局所の話です。接線は一次関数ですが、一次関数は線形関数とも呼ばれます。よって微分とは、局所線形化をやりたいという人間の意志だったのです。その意志を許容する条件が、微分可能という条件だった事になります。
しかし、いくら(11)の形が理論的に使い出があり、概念的にわかりやすくても(← 本当か?)、計算手段として不便だったら、数学屋さんたちはきっと別の表現に走ったはずです。じつは(10),(11)の形は、計算手段としても非常に使い出があります。それを見るためにここでは、四則演算に対する微分公式を証明してみます。
和の微分公式は、
でした。これを通常は、まだ見ぬ結論(12)の姿を想像しながら、
h≠0のとき、
とこねくり回し、(12)を導くわけですが、もし(10)を使えば、
と目標物であるf'(x)とg'(x)が、初っ端から眼前に現れるので、後は微分の定義式を思い出すだけで、
という変形はたやすいはずです。上記に条件(4)を使えば(12)です。なれてくると(4)は当然として、O(h2)の項は、書く事すらしなくなります(※)。
(※)
としてもよい。
そして、Ο(h²)の特殊な演算器則
を使うと、
h≠0のとき
としてよい。
上記の演算規則を使わずに、
としてもいい。
ランダウの記号を用いた算法では、記号の胡散臭さ、疚(やま)しさ――記事を書く数学屋が感じる疚しさ、後ろめたさ!!――を少しでも軽減するために、極限を→で表すことが多い(^_^;)。
(※終)
定数倍の微分公式も同様です。
積の微分公式は、(13)のO(h2)の項を失くした形で((11)の形で)辺々の積をとり、
最後は(14)と(11)の形を直接比べ、hで割る操作もlimh→0の操作も省略して結論を出しちゃいました(^^)。
面倒くさいので省略しますが、商の微分公式も合成関数の微分公式も、同じくらい簡単に導けます。最後に、ネコ先生の大嫌いなロピタルの定理を証明しましょう。
[定理-2]
関数f(x)とg(x)が、あるxでf(x)=g(x =0、かつg'(x)≠0とする。このとき、
[証明(?)]
(11)を使い、f(x)=g(x)=0かつg'(x)≠0であるから、
より明らか。
[証明終]
なれてくると(16)を見た瞬間に、それは明らか!と思えるようになります(^^)。しかしロピタルの定理の運用には、十分な注意が必要です。それはネコ先生の仰る通りなのです(^^;)。
(おまけ)
h≠0のとき、
オイラー法とシンプレクティック法 その2 [ddt³さんの部屋]
オイラー法とシンプレクティック法 その2
3.しかし逆変換はある
しかし(6)の逆変換はあります。それが(7)です。(7)を繰り返し用いれば、(p(1),q(1))=(1.12253,1.47121)から初期条件(p(0),q(0))=(0,1)に戻れるのは明らかです。時間反転とは、時間に関して解の進行を逆向きに進める事なので、これは時間反転ではないのでしょうか?。この定義からは、時間反転ではあるでしょう。
ただこれは、オイラー法による時間反転ではないのです。それは(6)と(7)でτを-τにかえただけの同じ計算手続きになっていないからです。ここから、オイラー法は時間反転を満たさない、と言われます。
4.シンプレクティック法
オイラー法の時間更新手続き(6)を、以下のようにほんのわずか変更します。
(11)
明らかに(11)の局所精度は、(6)とほとんど変わらないはずです。しかし(11)は等速運動させてから等加速度運動させる計算手順になっています。このような運動を許容する運動方程式はあるはずです。じっさいシンプレクティック変換の標準的構成法によって、(11)を与える事はたやすいです。たやすいのでその詳細は省略します(^^;)。そういう訳で(11)は、シンプレクティック変換です(^^)。
シンプレクティック変換ならば、何らかの運動方程式に従った運動を表します。その運動方程式は、(11)の局所精度が(6)と同等程度という事を考えれば、もとの系に対する近似運動方程式になるはずです。しかも近似エネルギー曲線も持つはずです。もとの系のエネルギー曲線の表式と近似エネルギーの表式を比較すれば、大域的精度までいっきょに明らかになるかも知れません。こうしてなんとなく性質の良さそうな近似解法への展望が開けます(^^)。
要するに問題の理論的背景は、(11)の計算手順が厳密解を与える運動方程式をみつける問題に帰着されます。シンプレクティック法は、少なくとも(2)のような自励系に対しては既に確立された方法です。オーソライズされた方法は非常にスタイリッシュで、少なくとも初見では、とても読みたくなるようなものではありません(^^;)。こんな具合です・・・。
まず(11)の1段目を、運動エネルギー1/2×p2に適当な偏微分作用素Dp()が作用したものとみなし、同様に2段目をポテンシャルエネルギー-1/2×q2にDq()が作用したものとみなします。Dp()とDq()は非可換作用素です。(11)を繰り返す計算手続き全体は、Dp()とDq()から構成される指数演算子を用いて、
と表されます。右辺のH(p,q)が求める近似運動方程式のエネルギーです。エネルギーがわかればその形から、運動方程式も導けます。普通の指数法則から、
と思いたくなりますが、Dp()とDq()は非可換なのでそうはならず、非可換作用素代数(リー代数)のBCH公式(ベイカー・キャンベル・ハウスドルフ公式)を使ってH(p,q)を求めます。BCH公式は、クソ力任せな計算で証明されますが、H(p,q)を求める計算もクソ力任せです(^^;)。
オーソライズされた方法は、(11)が厳密解を与える系を力任せに構成するという方法です。もう一つのやり方としては、(11)が厳密解を与える系を探すという方向もあり得ます。こっちの方が多少初等的(?)と思えるのですが、非常に地道な大量の計算をしなければならないので、やっぱり詳細は省略します(^^;)。こんな具合です・・・。
まず(3)のポテンシャルエネルギー:-1/2×q2を、一般にU(q)で表します。要はもとの系として、
の形のエネルギーを持つ自励系を想定する、という意味です。そうするとシンプレクティック変換(11)は、
と書けます。シンプレクティック変換(13)が厳密解を与える系を探すために、近似系のエネルギーに、
の形を仮定し、これに変換に対して運動方程式が不変になるという条件(S-2)を要求します。要するに(13)で(14)を(P,Q)に書き換えた時に、同じ形になるという要求です。
で、・・・初等的なムチャ長ぁ~い計算を実行すると・・・、(14)の未知関数f0,f1,f2,・・・に対して、次の偏微分方程式の系列が得られます。D1,D2を偏微分作用素、
として、
ここにn=0,1,2,・・・で、n=0のとき右辺は0とします。また0≦k,mに対し、
は、D1,D2に含まれる偏微分演算子の係数を、偏微分演算子に対して定数と考えてとった積です。(15)もなかなか力任せの計算になりそうですが、D1作用素が線形演算子なので、なんとか処理できます。
n=0の時は、
なので、f0として、
の形の解を選べるのは、けっこうすぐわかります。すなわちf0として想定した系のエネルギーの表式をとれます。n=1では、
ですが、右辺に(17)を代入し少し頑張ると、
を得ます。n=2で、
(17)と(18)を代入し、かなり頑張ると、
が出ます。以下同様に(15)を解き続ける事になりますが、3≦nでは右辺の項数も増えて、非常に非常に頑張らないといけない事になり、一般項を与えるためには、やはりクソ力任せ計算となります(^^;)。しかし(2)のような線形系においては、ここまでで十分です。
(17),(18),(19)でU(q)=-1/2×q2とします。そうすると、
(20)と(22)を比べると、(22)は(20)の定数倍になっているのがわかりますが、(16),(17)から明らかなようにこれは、斉次方程式D1g=0の特解で、1≦nで非斉次になる(15)の解の不定性を表すものだと解釈できます。従って(22)を採用するか否かは、人間が判断して良いものです。目的は、(12)に出来るだけ近い(14)を定める事でした。よって(22)の不定性を無視し、f2=0を採用します。
3≦nで(15)右辺のf0,f1に関する和は、pまたはqの3階微分以上の項の和になるので、それらは0です。n=3で右辺に残るのはf2に関する和ですが、f2=0と決めました。従ってD1f3=0となり、f3=0を解に選べます。以下帰納的に3≦nでD1fn=0が得られ、f2=f3=f4=・・・=0を選べます。
以上まとめれば、
が求める近似系のエネルギーです。近似運動方程式は正準方程式の手順によって、
と得られますが、(24)の厳密解が(11)で与えられるとわかってる以上、(24)を使うことはほとんどありません。それよりも、重要なのは(12)と(23)の比較です。(23)と(12)の差は、|pq|のτの1次オーダーです。最も荒い見当で|p2|~|q2|~|pq|程度とすれば、(11)では相対誤差をτの1次オーダーに保ったまま、どこまでも計算できる事になります。こうして、数値解の大域的性質がいっきょに手に入ります。
相対誤差なので今回のような単調増加する解では、絶対誤差は増え続け、局所精度ではルンゲクッタ法やオイラー法にさえ及ばない事態も起きますが、厳密解が大きくなりすぎ他の解法では計算不能になっても、シンプレクティック法は相対誤差を保ちながら計算可能な場合も良くあります。
現実には無限大に逝っちゃう解にはあまり興味がありません(少なくとも工学では)。重要なのはどこまでいっても有界な解です。このときポアンカレの再帰定理より、解は1度通過した点の近傍に必ず帰ってきます。という事は相対誤差一定なので、いったん増えた絶対誤差が減少もします。このような性質は、他の解法にはないものです。ここから特に周期解(振動系)において、抜群の性能を発揮します。こういう事もあり、シンプレクティック法に数値誤差の蓄積はないといわれます。
5.シンプレクティック法の宣伝
ちょっとシンプレクティック法の宣伝をしちゃいますね(^^)。1次の陽的シンプレクティック法の一般形は以下です。系のエネルギーを、
として、
ちなみにネコ先生のあげたリープ・フロッグ法は、2次の陽的シンプレクティック法の一種です。系のエネルギーが(25)のように、運動量だけの関数と変位だけの関数に分離されるとき、8次の陽的シンプレクティック法まで知られています(本当はどこまでもいける)。大域的精度はn次解法ならばτのn次精度です。
(26)は本当は陰解法です。時間更新手続きの両辺に(t+τ)が出てくるからです。にもかかわらず陽的と言われるのは、あたかもオイラー法のごとく順番に代入計算してけばOKだからです。これはV(p)やU(q)が非線形であっても同じで、非線形陰解法に憑いてまわる、時間更新ステップごとの繰り返し計算はいっさい不要になります。つまり非線形現象に対しても計算時間は、陽解法とほとんど同じです(高速)。かつプログラム簡単(Excelでもできちゃう!)。
しかも無条件安定の精度保証付きです!。陽解法の宿命である、解の発散や縮退は起こりません!。使ってみたくなりませんか?(^^)。
6.シンプレクティック法における時間反転
1次の陽的シンプレクティック法(11)を、初期条件(p(0),q(0))=(0,1),τ=0.1で計算したものが図-3です。シンプレクティック法の方がオイラー法より若干はずれた解を与えますが、近似エネルギー曲線にはしっかり載ってます。
(p(1),q(1))=(1.17309,1.48394)です。さらにτ=-0.1で時間反転を行うと・・・。
あれ~(^^;)もとに戻らないじゃん。しかもオイラー法の時より、はずれが大きい。・・・そりゃそうなんですよ。
τを-τにしたという事は、近似エネルギー曲線、
を、
に変えたって事なんですから。C=-0.32592は初期条件(p(1),q(1))=(1.17309,1.48394)から計算しました。オイラー法は異なる初期条件の解曲線に数値解の点列が飛んでく解法でしたが、シンプレクティック法では逆に、人間の方で異なる初期条件の解曲線へ数値解を乗せ換えてやったようなものです。実際にやってみるとわかりますが、アクア色の点列は、後の方のエネルギー曲線にしっかり載ってるのです(^^;)。
つまりシンプレクティック法でも解を逆進させるためには、逆変換を使うしかないという事です。逆変換は(27)になります。これは(11)でq(t),p(t)が左辺に来るように、移項しただけのものです。(27)を見ると、なんだか(11)と似ています。違いは、pの更新を先にやってからqの更新を行う点です。qとpの立場が入れ替わっています。じつはシンプレクティック変換では、変位qと運動量pの区別はないのです。qとpの事を人間が勝手にそう名付けただけ、という事になります。qとpを入れ替えるシンプレクティック変換が存在するからです。だから(27)もシンプレクティック変換なのです(^^)。従って(同じではないけれど)シンプレクティック法による、正確な解の逆進が可能です。この事態をさして、シンプレクティック法には時間反転性があるといわれます。
数値解法で時間反転を考えると、このようにちょっとわかりにくい事態になりますね(^^;)。
7.q'=qの場合
最後に、懸案の(p(0),q(0))=(1,1)のケースです(q'=qに相当)。結果は図-4になります。
q'=qの場合はエネルギー曲線が直線で、オイラー法の接線飛び出しが直線の傾きに一致して点列は正確に厳密解q=etのエネルギー曲線に載りますが、シンプレクテッィク変換ではないため、誤差は蓄積して行きます。シンプレクテッィク法では最初からエネルギー曲線がずれてますが、t=1ではオイラー法より厳密解に近いです。この場合の評価は微妙ですなぁ~(^^;)。
オイラー法とシンプレクティック法 その1 [ddt³さんの部屋]
オイラー法とシンプレクティック法 その1
以前このプログで、時間に関する1階偏微分方程式を時間に関する2階偏微分方程式に直して、強引にシンプレクティック法を適用してさんざんな目に合ったおぼえがあるので、とりあえず今回はどうなんだろぉ~とは思いました。また数値解法上で時間反転を考えると、ちょっとわかりにくい事態になるんですよね(^^;)。
ちなみにシンプレクティック法は基本的に、2階微分方程式に関する専用解法です。ここで2階微分方程式を一般に運動方程式と呼ぶとすると、運動方程式に従う運動はシンプレクティック変換を満たさねばならない事がわかります。具体的には位置と運動量の組を(q(t),p(t))として、時間がτ進んだ状態(q(t+τ),p(t+τ))への変換、
(q(t),p(t)) → (q(t+τ),p(t+τ))
はシンプレクティック変換になってるという訳です。
数値解法は一般に時間ステップ幅τによって、点列、
(q(t),p(t)) → (q(t+τ),p(t+τ)) → (q(t+2τ),p(t+2τ)) → ・・・
を与えます。この(q(t+nτ),p(t+nτ))から(q(t+(n+1)τ),p(t+(n+1)τ))を数値的に計算する過程がシンプレクティック変換になっているものを、シンプレクティック法と言います。
シンプレクティック法がリープ・フロッグ法(蛙飛び法)とも言われるのは、シンプレクティック変換の標準的構成法を使うと、変位場と速度(運動量)場を時間に従って交互に更新する計算になるからです。電磁場方程式の場合は、電場と磁場を交互に更新する蛙飛びです(^^)。
1.やってみますか
ネコ先生から出された宿題は、
でした。qを変位と呼びます。シンプレクティック法は2階微分方程式の専用解法なので、(1)を強引に2階微分方程式に持ち込みます。(1)を時間tで微分し、
すなわち、
になります。
(2)はエネルギー曲線を持ちます。両辺にq'をかけ、合成関数の微分公式の逆を使い、時間tで積分すれば、
が得られます。Cは初期条件で決まる積分定数。
ここでpは、
で定義される運動量(速度)です。エネルギー曲線(3)とオイラー法との関係を調べるために、(2)を連立1階微分方程式へ引き直します。
(5)の1段目は(4)そのものです。2段目は1段目の意味を考慮すれば(2)そのもので、じっさいに1段目を微分して2段目に代入すれば、(2)が得られます。(5)に対するオイラー法は時間ステップ幅をτとして、
でしょう。これはネコ先生が与えたオイラー法と本質的に同等な計算手順です。
ネコ先生の宿題は1階微分方程式(1)に対するものだったので、初期条件はt=0でq(0)=1だけでした。(5)は2階微分方程式(2)と同等なので初期条件は2個必要ですが、宿題に合わせれば(1)より、t=0でq(0)=p(0)=1です。しかしこの解ではオイラー法の特徴が見えにくいので、初期条件をq(0)=1,p(0)=0とします。
(6)でτ=0.1とした点列を(p,q)空間にプロットします。(p,q)空間は位相空間と言われ、微分方程式の解の定性的性質を見定める便利手段として、けっこう多用されます。図-1の計算はExcelで行いました。残念ながら、タイガー計算機(← 炊飯器じゃないよ(^^))などは用いておりません(^^;)。Cは-0.5となります。
時間範囲は0≦t≦1だったので、τ=0.1から(6)を10回繰り返します(赤点列)。実際に計算すると、(p(1),q(1))=(1.12253,1.47121)になります。さらにこれを初期条件とし、τ=-0.1として(6)を10回繰り返した時間反転結果が緑の点列です。最終結果は、(p(0),q(0))=(0.00000,0.90438)になります。という訳で、時間反転しても最初の初期条件には戻りません。
2.オイラー法の定性的(大域的)性質を調べる
どうしてこうなったかを考えるために、定量的な局所精度を調べるのではなく、運動方程式の不変性に注目して定性的(大域的)にオイラー法の性質を調べてみます。そこで基準になるのは、厳密解を与える解法の条件とはどんなものかという都合の良い考えです。あるとすれば理想的な解法の条件です。
(S-1) 変換:(q(t),p(t)) → (q(t+τ),p(t+τ))は、シンプレクティック変換である。
(S-2) 変換:(q(t),p(t)) → (q(t+τ),p(t+τ))は、運動方程式を不変に保つ。
(S-1)は冒頭で述べた、運動方程式に従う運動の必要条件です。(S-2)は、数値解の点列が同じ運動方程式を満たし続けるという条件です。(q(t+τ),p(t+τ))=(Q,P)と書くことにすれば、厳密解を与える限り、(q,p)は厳密な運動方程式、
を満たすはずです。同じ理由から、変換後も(Q,P)に対して厳密な運動方程式が成り立たねばなりません。
条件(S-2)さえあれば数値解法は、厳密解を与えそうな気がします。しかし同じ運動方程式を満たし続けたとしても、同じ初期条件から出発した点列であるという保証はありません。そこを制約するのが条件(S-1)です。じつはオイラー法は、条件(S-2)を満たすのです。
(6)は線形変換、
ですから逆行列をとる事により、
が得られ、(7)を時間tで微分すれば、
になるのは明らかです。
(7),(8)を運動方程式(5)、すなわち、
に代入すれば、
と(q(t+τ),p(t+τ))=(Q,P)に対して(9)と同じ形(10)が得られます。
よってオイラー法(6)がシンプレクティック変換になっていれば、オイラー法は厳密解を与える事になります!。
安心してください。ちゃんとシンプレクティック変換ではありませんよ(^^)。(S-2)を満たすオイラー法が、シンプレクティック変換なら(S-1)も満たすので、それは厳密解を与えるはずですがオイラー法の点列は図-1から、厳密解のエネルギー曲線から外れた点列になっています。よってオイラー法は自明にシンプレクティック変換ではありません。
つまりオイラー法は、運動方程式は満たし続けるけれど、異なる初期条件から出発した解にジャンプし続けているのです。実際t=0,0.1,0.2,・・・に対するオイラー法の点列(q(t),p(t))について、エネルギーの表式である(3)を計算してみると図-2となり、オイラー法においてはエネルギーが保存しません。そしてその理由は明白です。何故ならオイラー法は計算手順(6)より、エネルギー曲線の接線方向へ飛び出し続けるというくせを持つからです。そしてその理由も明白です。
(6)の時間更新手続きはqもpも同時に更新するので、等速運動させながら同時に等加速度運動もさせるというものです。そのような運動は原理的に不可能であり、当然それを許容するような運動方程式もあり得ないので、オイラー法がエネルギー保存則を満たす訳はないのです。
だったら図-1の時間反転部分も同じじゃないですか!。オイラー法は定性的(大域的)に考えると、数値解がどこにジャンプして飛んで行ってしまうかわからない解法なのです。時間反転しても、最初の初期条件にもどる保証はありません。
新潟山形県沖地震の揺れの解析 その2 [ddt³さんの部屋]
4.基線補正の検証と周波数特性
図-13に、基線補正後の加速度波形のEW,NS,UDを示します。これと原波形の図-1を比較しても、どこが変化したのか、ほとんどわからないと思います。
もしぱっと見で違いがわかったら、明らかに補正の(データ加工の)やり過ぎです。違いが見えないと思いますので、周波数領域で比較します(図-14)。
図-14では、補正前の加速度の振幅スペクトルを黄色いラインで示していますが、補正後の青ラインと重なってほとんど見えません。基線補正が地震波の特性に大きな影響を及ぼしていない事がわ
海溝型の地震の主要周波数は1~10 Hzと言われています。この地震の周波数特性ですが、まず第1ピークが5~10 Hzにあり、高周波の地震です。また第1ピークとあまり違わない大きさの第2ピークが10~11 Hzの範囲にあり、直下型の地震にかなり似ています。この前の北海道胆振東部地震(直下型)でも、同じようなスペクトル分布が得られました。今回の新潟山形県沖地震と、海溝型の典型である東北地方太平洋沖地震(3.11)を、パワースペクトル(振幅2乗値)で比較したのが、図-15です。
3.11はM9.0という超巨大地震ですので10 Hz以上のパワーもすごいですが、無視できないパワーを持つ周波帯が、1 Hz以下まで満遍なく広がっているのが特徴です。そういう意味では3.11は、低周波な地震です。総出力は概ねパワースペクトルの面積に比例すると考えられるので、3.11のパワーは圧倒的です。ちなみに牡鹿は3.11の震央から最も近い地点で121 kmです(^^;)。
新潟山形県沖地震は、直下型に近いなぁ~という印象です。ただし3.11の深さはこの規模の地震では非常に浅いと考えられる24 kmですから、直下型の性質も持ってるようです。
5.変位軌道
図-12の変位結果から、変位軌道を描いたのが図-16です。こんな風に揺れました?・・・ってわかるわけないか(^^)。
5.距離減衰
震源(震央距離)と加速度振幅の2乗値(1乗値)の両対数は、良好な線形相関を持つと言われています。その相関を表す回帰直線が、距離減衰式です。ここでは震源距離dと振幅2乗値pを用います。振幅2乗値は意味として、波動のパワーに対応するので。図-17です。図中の黒点線が、距離減衰式を表します。少々データ処理は行いましたが、どちらの直線も相関係数のR2値は0.6以上あります。回帰直線は、
log10 p=-1.6546・log10 d+1.2586,d ≦ 100 km.
log10 p=-5.1687・log10 d+8.2348,100 k m≦ d.
回帰直線が2つに折れるのは、直下型の特徴だと自分は思っています。一方、海溝型の典型である東北地方太平洋沖地震では、
という風に、一つの回帰直線で表せるケースが多いです。
距離減衰からみても、今回の地震は直下型に近いのではなかろうか?というのが、自分の意見です。実際に冒頭のgoogle earthからも明らかなように、陸地から非常に近く非常に浅い所で起こった、海底地震です
6.まとめ
なんか四川省の直下型地震と規模も深さも変わらないくせに、どうして日本の方が被害が小さいんだ?などと、中国がうらやんでるみたいな記事も見たのですが、自分の意見では偶然ですよ。本当の直下型でなかった、というだけの偶然。
加速度振幅の2乗値の対数と計測震度の間には、じつは非常に良好な相関があるのです。どちらも地震のパワーを表す代表値だからです。
加速度の振幅2乗値と計測震度には、R2=0.9294の相関があります。計測震度の計算法は、上手くできてるなぁ~と自分は思います。
図-19の回帰直線は、sを計測震度として、
s=0.8569・log10 p+5.5601
(5)
です。
もし新潟市が震央だったらと仮定してみます。要するに震源の直上です。その時の震源距離は、深さ10 kmです。これは(4)で、log10 d=log10 10=1という事です。よってlog10 p=-0.396。(5)よりs=5.22。震源距離14 kmの温海の実測震度は5.2ですから、この値は妥当そうですが、回帰直線はあくまで統計量です。
図-17を見ればわかるように、実測データは距離減衰式のまわりでばらついています。なので青ラインで示したように、実測データを包絡させた方が安全なのです。青ラインは(3)で表される回帰直線を上へ1だけシフトさせたものです。でも対数目盛ですからね。上へ1という事は10倍の差です。振幅2乗値で10倍ならば、振幅値では3.2倍の差になります。温海の最大加速度振幅は653 galですから、653×3.2=2065 gal~2.1 gにもなります。この値は日本に残っている強震観測記録の最大値の約2倍です。
このように距離減衰式は、対数目盛マジックにおんぶにだっこの方法なのです(^^;)。さらに図-19にも±1程度のばらつきがあります。これも青ライン側へシフトさせて考えるのが安全です。そうすると最悪の事態では、log10 p=-0.396+1=0.604,s=0.8586×0.604+5.5601+1=7.079~7.1。でもまだ甘いかも知れない。
2014年10月23日に起きた新潟県中越地震は、M6.8,深さ13 kmで今回の地震に似ています。ただしこれは完全な直下型。最大計測震度は7を記録しました。震源の直上にいなかった幸運に感謝すべきだと、自分は思います。もしいたら、こうなっていてもおかしくなかったと思う。
図-20 北海道胆振東部地震の斜面崩壊,厚真町
https://ja.wikipedia.org/wiki/%E5%8C%97%E6%B5%B7%E9%81%93%E8%83%86%E6%8C%AF%E6%9D%B1%E9%83%A8%E5%9C%B0%E9%9C%87#/media/ファイル:Ground_surface_in_northern_Atsuma_Town_after_2018_earthquake.jpg
こういう斜面崩壊が何10 kmの長さで、延々と起こった。
1995年1月の兵庫県南部地震(阪神淡路大震災)の少し前、1989年10月にサンフランシスコ地震(ロマ・プリータ地震、M7.1,深さ18.5 kmの直下型)が起き、倒壊する高速道路から落下する車両という衝撃の動画が報道されたものでした。
この時、日本の耐震業界の大御所は言いました。
「日本の耐震技術はアメリカより進んでるので、わが国でこんな事は起こらない」
・・・と。
図-21 兵庫県南部地震 名神高速道路(大阪灘区)
一般国民を安心させるための方便だったのかも知れませんが、特殊民である我々は、
「あ~言っちゃったよ、この先生」
と思ったものでした。ほどなくして兵庫県南部地震(M7.3,深さ16 kmの直下型)が起こり、高速道路の丸ピアが次々と根こそぎ倒れたのは記憶に新しいです(私には)。
その後、ものすごい速さで復旧仕様が出され、耐震基準はほぼ年単位で改訂されて行きました。その時の中心メンバーの一人は、先の大先生でした。
直下型の現場にはいたくないですね。
新潟山形県沖地震の揺れの解析 その1 [ddt³さんの部屋]
新潟山形県沖地震の揺れの解析 その1
先日の中国四川省の直下型地震で、M6.0,深さ「10 km」はやばい!なんてコメントしたら、18日に来ちゃいましたね。それもM6.8,深さ「10 km」が・・・。それでお見舞いというわけではないのですが、地震波を解析してみました。
1.震源情報
気象庁発表によれば、発生時刻は2019年6月18日22時22分、震央は北緯38.6°,東経139.5°、M6.8、深さ10 kmです。ここでは防災科学研究所の強震観測網K-Netのデータで、震源に最も近い温海の観測ステーションの記録を使います。震央距離は8 kmです。ところで新潟は・・・震央距離約87 kmの計測震度4となっておりました。
2.地震加速度波形
温海観測ステーションのEW(東西),NS(南北),UD(上下)の地震加速度波形を図-1にあげます。サンプリング周波数は100 Hz(0.01秒)、観測時間は148秒です。
東西方向の最大加速度は約600 gal、南北方向は600 gal超え、上下方向でも200 galです。1 gal=1 cm/s2ですから、水平方向には概ね0.6 g/cos45°、上下方向には0.2 gとなります。つまり水平方向に自分の体重の約8割がたの力で引っ張られ、同時に子供一人を肩に載せたくらいの衝撃力を、0.01秒間に食らった事になります。インタビューで誰かが言ってましたが、とても立ってはいられませんね。
3.基線補正
加速度記録(cm/s2)があれば、1回積分して速度(cm/s)、2回積分して変位(cm)です。数値積分にはふつう台形公式が使われます。台形公式は振動と周期関数のたぐい対しては精度が良好な事、加速度記録はデジタルデータなので折れ線近似以上の情報は得られないから、というのがその理由です。しかしそのまま積分を実行すると、積分結果は見事に発散するのが普通です。そこで基線補正といわれるデータ補正を行います。ただし補正のやり方には諸説あって、ここでは一番単純な方法を採用します。
まずNS加速度をそのまま積分した結果が、図-2です。
積分結果は、速度の最大が375 cm/s,変位の最大が27551 cm。最終速度は時速13.5 km/h,移動距離は275.5 mになります。明らかに嘘ですよね?(^^;)。
グラフを見直すと、速度はほとんど直線で等加速度運動状態です。対応して変位はほとんど放物線です。地震で等加速度運動はあり得ないでしょう。だって図-1の地震の主要動は概ね16~30秒の範囲で、それ以降は明らかに終息に向かってます。にも関わらず速度と変位がずっと増え続けるというのは、納得できません。宇宙空間じゃあるまいし・・・。
加速度測定値に一定値の誤差があると考えられます。これを零点シフトと言って、現実の計測では必ず起こります。零点シフトの値を推定します。
零点シフトは地震波がない時の測定値の状態を調べれば良いはずです。図-1のNS加速度を見ると、計測開始の最初の16秒間ほどは地震波が来てないので、そこの平均値をとり時系列全体から引けば良いと考えられます。しかし地震波の正確な始まりの時刻はどこでしょう?。もちろん厳密にはわかりません。また16秒くらいでやっても問題なさそうにも思えます。しかし数値積分はちょっとしたゴミのために図-2のような結果を招きます。出来る限りの事をするために、ここでは以下のようにします。
K-Netの加速度計仕様によれば、分解能は0.001 galです。0.001 gal以内の変動は、+0.001 galとも0 galとも-0.001 galとも記録される可能性があります。よって有意な変動は0.002 galよりも大きなものと考えられます。図-3は0.01秒ごとの加速度変動を計算し、0.002 galよりも大きければフラグ1を立て、以下であればフラグ0を立てて時系列順に並べたものです。
有意変動は計測時間中に万遍なく分布していますが、非有意変動がほぼ連続的に密集分布している時間帯が、計測開始からしばらく続きます。その範囲と図-1のNS波形を見比べると、この部分は地震波がまだ来てない時間帯と判断できます。そこでその時間帯の最後の時刻t0=15.94秒を波形開始時刻と決めます。0~t0にあるデータの平均値は、ajを加速度データとすると台形公式で、
(1)
で求められます。a1=a(0),aN0=a(t0)です。A0の値を加速度時系列全体から引き、積分をやり直した結果が図-4です。これを[基線補正-1]とします。
地震波が来ていないと考えたt=t0までは確かに妥当です。速度,変位ともほとんど0をたもってます。しかし地震波が始まると、また等加速度状態に見えます。という事は波形開始後は、それ以前と別の零点シフトになったのでしょうか?。零点シフトが零点ドリフト(漂流)に移行する事も、じつはごく普通に起こります。
加速度計は電気回路ですが、入力が小さい時と大きい時で回路の特性に微妙な変化が生じるのが現実です。それでも最終速度は15 cm/s,移動距離は9 m程度と大幅に改善されてます。人の歩行速度は1 m/sほどなので、[基線補正-1]は確かに効果があったと考えられます。ちなみにA0=-2.435 gal~-0.0025 gのゴミですが、速度と変位には大きな影響があります。この程度の加速度が本当に効いてるとすれば、無感地震程度で日本列島は移動しまくりです。
図-4の速度波形を参考に波形開始後の零点シフトも開始前と同様に考えると、測定後半に速度の傾きとしてそれが明確に現れてるように見えます。図-4からは100秒以降で速度は十分直線的かな?と見えるのですが、NS波形の時系列全体を示すと、図-5です。
100秒以降にも小さな有意振動が見えるんですよね。そこで非有意変動が時系列の終端でほぼ連続的に分布するようになるまで、有意の基準を上げてやります。
図-6で時系列終端にある非有意データの密集時間帯の最初の時刻t1を、いちおう地震波形の終わりとみなし、それ以降の平均A1(計算法はA0と同じ)を、加速度時系列のt0以降のデータから引きます。A0ではなくA1を引きます。これを[基線補正-2]とします。
再計算結果には、まだ少し等加速度状態が見えます。つまり0≦t0とt0≦t≦t1とt1≦tでは、本当に零点シフトが違うという事です。本当に零点ドリフトが起きてます。この中で、いちおう確実と思えるのは0≦t0の零点シフトだけです。それでも終端速度は2 cm/s,移動距離は2 mになりました。A1=-2.558 galで、A0=-2.435 galとほとんど変わりません。それでもこれだけの違いがあります。
こういう事態になると、もっと強い仮定を持ち込むべきだという意見も出てきます。今までやってきた事は、全てデータ処理に先立つ事前情報(仮定)に基づいています。地震が終息に向かえば速度と変位が増加し続けるわけないとか、計測器は時間的にほぼ一定な零点シフト誤差を必ず持つとかは、みなそうです。根本的な話として、加速度のデジタルデータの測定値を折れ線でつなぐのもそうです。これは、自然は少なくとも連続であろうという仮定です。より強い仮定の代表は、次の二つです。
1) 地震の終端で速度は0になる。つまり地震が終われば地盤は動かない。
2) 地震の終端で変位は0になる。つまり地震が終われば地盤は原位置に戻る。
一見問題なさそうな(というか当たり前の)仮定に思えますが、地震の終端を計測データの終わりにとるしかないという問題があります。計測データの終わりで地震が止んだ、なんていう保証はないんですよ。図-5をみると、なんか150秒以降も地震は続いてるようにも見えますが、1)を適用してみると図-7が得られます。この方法は、時刻t0以降のデータ平均をt0以降のデータから引いた[基線補正-2]と同じです
図-7の結果については「嘘っくせぇ~!」です(^^;)。何が気に食わないっていって、まず速度が(わずかですが)等加速度状態で0に収束するところです。自然現象がそんなに綺麗に挙動するわけねぇ~だろ~!。対応して変位は、時系列の終端でちょうど頂点に達する放物線になります。
「嘘だこれは!」「人工的過ぎる!」。ただし[基線補正-1)]で納得できる結果が得られるケースも、けっこうあるのは事実です。
さらに仮定2)まで適用すると、もっと人工的な結果になる事が多いです。仮定2)には明確な反証があります。例えば東北地方太平洋沖地震(いわゆる東日本大震災3.11)では、地盤が最大5.85 m動いた(動いたまま)と国土地理院がGPS観測に基づいて、はっきり認めています。
そろそろデータ処理の結果としてどの辺りをめざすのか、決めるべき時です。まず計測データの終端で地震が止んだとは限りません。従って時系列終端で速度と変位が0になるとは限らない。しかし収束傾向は見えて欲しいし、終端速度と変位はあまり大きな値にもなって欲しくない。こんなところでしょう。というか、その辺りが精一杯では?(^^;)。
いずれにしろ、図-6の等加速度状態を除く算段を考える必要があります。この等加速度状態は、t0≦t≦t1とt1≦tでの零点シフトのわずかな差から、t0≦tで平均的に起きた等加速度状態だと考えられます。それを言ったら、t0≦t≦t1とt1≦tの零点シフトというものだってそれぞれの区間の平均値でしょう。正確な零点シフトの時間依存性(厳密な零点ドリフト)なんてわかりっこありません。
でもいま実用的に問題になってるのは、t0≦t全体の零点ドリフトの平均値としての零点シフトです。そこで気づくのが、図-4のt0以降の変位が、ある2次関数に載りそうだという事です。変位の近似2次関数の2次の係数の2倍は、その等加速度状態の加速度です。図-4の変位のt0以降を取り出し、近似2次関数を計算してみると、以下になります。
相関係数はR2=0.9999972。ほぼ確実な近似です。この異常な相関性の良さが、零点ドリフトなどの系統誤差の存在を示唆します。
(2)のtの2乗の項の係数の2倍が、零点ドリフトも考慮したt0≦t全体における平均的な零点シフトA1'=-2×0.0562446=-0.1124892 galだと思われます。この値を先の[基線補正-1]の加速度時系列のt0以降のデータから引き、[基線補正-2']とします。
予想通り、ほぼ等速運動状態になりました。ちなみにA1=-2.558 gal、A0+A1'=-2.548 gal。両者の差は、わずかに0.01 gal=0.1 mm/s2。たったこれだけのゴミで、図-6と図-9の劇的な違いが生まれます。
さて図-9の等速運動状態の物理的解釈です。これはもちろん、(2)のtの1次の項に対応した等速度です。
もしt0≦tで今度は直線近似を行えば、再び異常な相関性を示すでしょう。これを系統誤差と考えると速度の誤差なので、加速度波形開始からの加速度誤差の累積です。という事は時系列の後半はほとんど等速度で間違いないと思われるので、誤差の累積は波形開始の最初の部分で起きているはずです。特に疑われれるのは加速度波形の主要動部分です。主要動から後の時間では、等速度誤差に影響する誤差の累積はほとんどないと考えられます。
時系列の終端t2から逆向きに変位データをtまでたどり、t~t2間の回帰直線を計算し変位データとの相関をとり、その変化をtのグラフとして表したのが図-10です。予想通り、異常な相関係数の高さです。最低値は0.9986になりました。終端t2付近の相関係数の低下部は、変位変動がデータ数の少なさを上回った部分で、あまり意味がありません。十分データが増えると直線性が回復し0.9998~0.9999の間で増減します。データの始まり方の増加は、波形の主要動の影響が大きく単なる偶然でしょう。ここで注目するのは、直線性が回復したとみなせる、図-10のt=t3の時点です。さっきの話が妥当であれば、この時点の傾き(速度)v0はかなり信頼できる等速度誤差の値とみなせます。t3=119.58秒,v=0.780254 cm/s。
カタツムリよりも遅いかもしれません(^^)。でもt3=119.58秒のタイミングは、図-6の結果であるt1=120.11秒とほとんど同じです。
図-11は、さっきとは逆にデータの始端からt0~t間の回帰直線を計算し、その傾きとv0との差をtのグラフとして表したものです。ここで注目するのは、データ始端側で相関係数が最大になる時と、回帰直線の傾きがv0に最も近くなるタイミングt4です。さっきの話が妥当であれば、その時点までの加速度誤差の累積で等速度誤差が生じたと考える事は可能です。図-11では理想的に両者がほぼ一致してますが、相関最大時の方を優先します。t4=30.33秒となります。この結果と図-1のNS波形を確証バイアス付きの目で比較してやると、t0=15.94秒~t4=30.33秒の範囲、概ねt=16~30秒の範囲が地震の主要動に見えてきたりします(^^;)。
t4~t2の回帰直線の傾きはv=0.780254 cm/sになります。[基線補正-3]として、
終端速度と終端変位は、-0.111 cm/sと2.404 cm。この結果の評価ですが、震央距離にして8 km、深さにして僅か10 kmという至近距離でM6.0なんていう大規模地震が起こったんですから、本当に新潟県と山形県の沿岸部は2 cmくらい「ずれて止まった」のかも知れません。しかしデータ終端で地震が止まったという保証はないので、まだ変位中かも知れません。そういう思いで変位グラフの末尾を見ると、少しですが減少傾向(収束傾向)が見てとります。また終端速度と終端変位の値は、決して大きなものではありません。この状況は求めていたものです。自分には、かなり自然な結果に見えます。
同様な方法で基線補正を行ったEWとUDの結果も、図-12に示します。
梁の方程式に至る長い過程 その3 [ddt³さんの部屋]
梁の方程式に至る長い過程 その3
5.梁の微分方程式はどこに?。
梁の曲げの微分方程式の話はどこいったんでしょう?。じつはじつは純粋な梁理論は、梁の曲げの微分方程式なんか扱わないのです。梁の曲げの微分方程式の話は、取って付けた話になります(^^;)。というのは応力(力)ベースではなく、変位(変形量)ベースで考えた方が便利なこともけっこうあるのです。設計業務で変形後の梁の変位量を出したい!(それは変形量から計算できる)という場面は、確実にあります。
じつは変形後に断面がどうなるか?という話も、歪み分布と同じなんですよ。そして発想も超単純。
梁は上面圧縮,下面引張なんだから、断面の上端は左へ変位して上側は縮み,下端は右へ変位して下側は伸びたんじゃない?と単純に考えてみます。微小変形理論でした。「微小入力に対しては、出力は微小入力に比例する」から、図-14のように上下の変位を直線でつないで良いでしょう(D/Lが十分小さいなら)。そうすると変形後の断面と変形前の断面の交点は、変位0の位置です。伸びも縮みもしない位置。これ中立軸ですよね?。だって中立軸はそうやって決めたんだから。
よって上下端がともに右へ、ともに左への可能性は排除されます。もしそうなら直線で上下端の変位をつないだとき、中立軸を通れないからです。超単純に考えたら変形後の断面の形までわかっちゃった(^^)。
次にやっぱり微小変形理論なのです。図-14の断面の傾きもオーバーです。図では変形後の断面は中立軸と目に見える角度を持ってますが、微小変形理論なので変形後も断面は、中立軸とほとんど直交してるはずです。微小変形なので、目に見えちゃいけないのです。
・変形後も断面は直線形状を保ち、中立軸と直交する。
これをキルヒホッフ・ラブの平面保持の仮定といいます。いかにもな「取って付けた仮定」と思えません?(^^;)。キルヒホッフさんは電気回路の電流保存則で有名なキルヒホッフさんです。ラブさんはそのすじでは有名な、ラブ波(表面波)のラブさんです。
変形後の中立軸の形状をy=y(x)とします。それを梁の変形曲線と呼びます。変形後の断面も中立軸と直交するのでした。そうすると図-15に示した微小距離dxで梁を切り出してやると、図-16になります(少々オーバー描写です)。
変形後の断面は中立軸と直交するので、微小距離dxで梁は、厚さDの円弧の一部になります。位置xで中立軸に接する円の半径をr(x)とします。中立軸の長さdxにおける(弧長dxの)円弧の角度は、
従って中立軸から縁端距離d1,d2における弧長は、
(17)と、伸び縮みしない中立軸の長さdxとの差が、距離dxにおける上下面の変形量です。
歪みは変形勾配でしたので、(18)をdxで割ったものが上下面の歪みです。
上下面の歪み差は、
従って図-16における、断面上の歪み勾配は、
になります。これは(12),(13)で算定された歪み勾配aです。
6.曲率による曲線の表現
曲率と曲率半径の正式な定式化はネコ先生の、
https://nekodamashi-math.blog.so-net.ne.jp/archive/c2305520799-1
などをお読みください。ここではddt^3風の曲率と曲率半径の定式化をご紹介し、最後に梁の曲げの微分方程式が、やっと出てきます(^^;)。
まず普通に半径rの円を考えます。
図-17のように半径rの円では、基準ラインからθ回った時の円弧の弧長sは、s=rθなので、円の曲率1/rは、
とすぐに得られます。曲率1/rは、θとsが比例するので円では当然一定で、それが外法線方向の角度を表すθと弧長sの比例定数になってます。
次のように考えたいんです。
・比例関係があれば、対応する微分による一般化がある。
例えば1次関数y=axでは、yはxに比例し比例定数aは一定です。こういう考えは駄目ですか?。
直線ではyの変化率aは場所xによらず一定なので、変化率dy/dxが場所xごとにf(x)と変化するのが、一般の曲線。同様に。
円ではθの変化率1/rは弧長sによらず一定なので、変化率dθ/dsが弧長sごとに1/r(s)と変化するのが、一般の曲線。
ではないだろうか?・・・と(^^;)。この対応は妥当だと思うんですよ。だって(x,y)座標系による曲線の表現は人間の勝手であって、本来は曲線をどのような座標系で表してやっても良く、それによって曲線の形が変わるわけないからです。
そうすると(x,y)系における場所xと、(s,θ)系における場所(?)sとの対応が気になります。曲線のの形がどんな座標系においても変わらないなら、xとsは何らかな意味で、曲線の「同じ場所」を指定する必要があります。こうですよね?。
微分方程式、
は、場所xごとに曲線の進むべき方向dy/dxをf(x)でコントロールするので、ある曲線y=y(x)を定義します。これは1階微分方程式の基本的な考えです。
微分方程式、
は、(s,θ)空間で場所sごとに外法線方向θの進むべき方向を1/r(s)でコントロ-ルし、θ=θ(s)を定義します。図-18からsとxは曲線の同じ場所を指します。外法線方向θは常にその同じ場所で、接線方向dy/dxと直交するので、これは場所x(もしくはs)でdy/dxをコントロ-ルするのと同じです。
(24)は(23)を座標変換しただけものと思えてきませんか?(微妙に違うけど(^^;))。なお曲線の内と外は最初に定義しておく必要がありますが、円にあわせて曲線の凸側を外と定義します。
(24)による曲線の表現は、曲線の自然方程式と呼ばれます。この表現の利点は、座標系に依存しない事です。弧長sは曲線の形だけで決まり、任意の座標系で同じです。外法線方向θの値は基準ラインの方向に依存しますが、(24)で必要なのはその増分dθだけで、これは基準ラインの方向に依存せず座標系と無関係です。そして曲率半径r(s)も座標系に無関係に決まってます。
(24)は確かに理論的には美しいのかも知れませんが、実用的な計算ではやっぱり(x,y)座標がわかりやすいんですよね。それで(24)を(x,y)系へ再翻訳するのです。
です。
dθ/dxを考えると、θ方向は常にy'(x)方向と左回りに直交するので、dθ/dx=dy'/dx=y"(x)。従って、
再び微小変形理論を使うと、(y'(x))2~0かつs~x。よって(24)の微小変形理論での翻訳は、
ちなみに曲率による表現(22)も土木では広く知られており、実際に曲率の概念で「お仕事」する時もあります。
梁の微分方程式に至る長い過程 その2 [ddt³さんの部屋]
梁の微分方程式に至る長い過程 その2
4.梁理論
断面力から断面応力を算定する段階です。それには梁理論が必要です。
図-6はあんまり細長くないですが、細長い棒と思って下さい。棒という用語も正確ではないのですが、とりあえず棒と言っときます(^^;)。上から鉛直荷重を受ける細長い棒は、図示したように上面側が縮み(圧縮状態で)、下面側が伸びて(引張状態)で「曲がる」ことが経験的に知られています。もっとも微小変形理論ですので、目に見えない程度の曲がりです。でも目に見えないと絵にも描けないので、図-6はオーバーに描いてます。
もう一つ経験的にわかってるのは、図のy方向の変形は無視できるです。(3)からわかるように物体のバネ定数はバネの長さが小さければ小さいほど大きくなります。図-6のy方向を支えるバネは、棒の長さLに比較して棒の深さDが小さければ小さいほど、相対的にx方向のバネより格段に強い訳です。ここからy方向の変形は無視します。もっとも大きい/小さいを言うためには、基準スケールを決めなければなりません。それで「相対的に」と言ったのですが、鉛直荷重で潰れた棒なんて、見た事ありませんよ(^^)。
あと図-6のx-y軸に直交する紙面直交方向のz軸があります。ありますが、紙面奥行き方向の状態は一様と仮定します。何故なら奥行き方向の荷重は考えてないからです。
以上まとめれば図-6の棒は、一次元のxだけで各部の状態が決まる一次元部材です。これを線材と言ったりします。
図-6の断面状態を考えます。図-6は長手方向(x方向)の挙動だけを考えれば良いので、水平方向のバネの集まりである図-3でモデル化できます。単純な棒の伸び縮みと違うのは、断面の高さ方向yによって伸び縮みの具合が変化する点です。上は圧縮,下は引張ですから、図-3に準じてモデル化すれば、
みたいになるでしょう。図-7に対応して、応力-歪み関係(6)を考慮すれば、断面上の歪み分布は、
みたいになるはずです。ここで引張歪みの値を正と決めます。圧縮歪みの値は負です。理由は、図-2でバネを伸ばす力Fを正にするのが普通だからです。バネが伸びるとは、バネは引っ張られてるって事です。
そして物理現象に関する暗黙の仮定を使います。全ての物理現象は連続でした。断面上の歪み分布は、上は圧縮で下は引張です。すなわちε(y)は断面上で負から正に変化する連続関数です。中間値の定理より、圧縮でも引張でもない(縮みも伸びもしない)、ε(y)=0の位置があります。その位置での材料の変形は0です。その位置を中立軸と呼びます。
中立軸がなぜ必要かというと、断面応力を算定するのに必要な断面2次モーメントの計算のために、欠かせないものだからです。
そして物理現象に関する暗黙の仮定を、もう一回使います。「微小入力に対しては、出力は微小入力に比例する」を。さっき言ったように、大きい/小さいを言うためには基準スケールを決める必要があります。では設計計算の目的はなんでしょう?。
それは構造系全体の任意の位置に発生する応力が、材料強度を越えない事を確認するのが目的です。目的は構造系全体なんですよ。つまり基準スケールは構造系全体です。図-4や図-6の橋長Lです。それに対して図-6のDが十分小さければ、それは設計計算の目的にとって微小と言えます。なのでD/Lが十分小さければ、図-9の近似で十分だろうとなります。具体的なイメージは図-5です。橋全体を一目で見渡せるくらい遠目で見てる状態です。図-5くらい細長ければ(D/Lが十分小さければ)、図-9の近似でもOKという気はしませんか?(しないかも知れない(^^;))。
歪みは単位長さあたりのバネの伸び(縮み)でした。バネが伸びたり縮んだりするからには、断面に力が作用してるはずです。図-7,8は図-5の左部分系に対応するものなので、断面に作用する力は右部分系に由来する断面力です。左部分系は、断面力M,S,Nと支点反力R1に支えられて荷重P1に対抗するのでした。
歪みに弾性係数Eをかければ応力です。応力と断面力の関係をつかむために、図-3のように断面上で一様に力Fで引っ張られてる状態を考えます。応力歪み関係(6)に断面積Aをかければ、
でなけりゃいけませんよね?。だってσ=F/Aで定義したんだから。そしてFは軸力Nですよね?。この関係は図-7でも同じですよね?。違いは、図-9のように断面上の応力分布が直線変化するだけです(^^)。
σ(y)=Eε(y),ε(y)=a(y-e)は、断面上での歪みの(応力の)直線変化を表します。a,eは未定定数です。単位面積当たりの力の密度が応力なんだから、それを面積で積算すれば軸力Nになるという式です。ここで断面形状を表す領域Rは、
・・・となったりするんですよ。だから2重積分なんです(^^;)。未定定数eは中立軸の位置を表します。y=eでε(e)=a(e-e)=0だからです。図-10のT型断面では、上側のz方向の幅が広いので、中立軸位置は上側に引っ張られます。
断面力の算定結果より、今はN=0でした。(10)でN=0とし整理すれば、
(11)
中立軸位置の算定結果である(11)の結果は、y方向の重心位置を計算する式と同じです。従って中立軸は、断面重心を通ります。ちなみに(11)の分子を、断面1次モーメントと言います。yの1乗の積分だからです。
同様に、曲げモーメントMに対応する断面応力の積算をとれば、未定定数aを決定できます。Mは回転力でした。従って対応する積算は、σ(y)による力のモーメントです。
ここで注意すべきは、σ(y)による力のモーメントを計算する際の回転中心の取り方です。さっき静止物体の回転中心はどこにとってもOKよ、と書きました。そして図-9も静止物体には違いないのですが、図-9のように断面が傾くという事は、現在の位置まで断面が回転したという事です。その間は静止物体ではありません。図-9は、断面が動き終わった変形後の状態という訳です。
※
図-9は歪み分布(x方向の伸び縮みの変形勾配の分布)なので、図-9から断面が傾く(回転する)というイメージをつかめない人もいると思います。これについては次節で再論するので、とりあえず図-9のように断面が回転すると思って下さい。
でも知りたいのは、断面を動かした力のモーメントです。よって本当の回転中心を選ぶ必要があります。それは明らかに中立軸です。とすれば、y軸の原点を中立軸位置にとるのが、最も便利です。
図-11の座標系で、σ(y)=Eε(y)=E・ay。よってσ(y)による力のモーメントは、
(12)
(12)左辺のy×の符号は、図-11の応力分布σ(y)が左回りの回転力を定義するので、左辺全体の符号が正になるように決めてます。これによってMの正方向と一致します。
歪み勾配aの表式の分母は、断面2次モーメントと言われ記号Iで表します。yの2乗の積分だからです。
IはInertia:イナーシャのIです。よって断面上の応力分布は、中立軸位置を原点にとると、
となります。
(14)より断面上の最大引張応力と最大圧縮応力は、断面の上下縁です。これらの値を材料強度と比較することになる訳です。そこで中立軸から上下縁までの距離を縁端距離と呼び、図-12のように表す事も多いです。縁端距離d±を使えば(14)は、
となるので、計算に便利なように断面係数を、
で定義しておきます。断面力計算から出てきた曲げモーメントに断面係数をかければ、速攻で最大/最小応力という仕掛けです(^^)。
以上は、構造系に水平力が働かず、軸力N=0が前提でした。この状態を純曲げ状態と言います。図-6に示すように、鉛直荷重でただ曲がってるだけだからです。ただ曲がってるだけの状態で断面に作用する回転力だから、Mを「曲げモーメント」と呼びます。土木の命名規則は、モロかマンマです(^^;)。
水平荷重だけが作用すれば、純引張(圧縮)状態です。その取り扱いは図-3で、たんなるバネの伸び縮みです。これも図-5のスケールで、D/Lが十分小さければそれで良いよね?、という発想です。
そういう訳で(14)で表される曲げモーメントMに関連する応力は、曲げ応力と言われます。一方σ=N/Aで与えられる軸力に関連する応力は、軸応力と言われます。曲げ応力と軸応力が両方働いたらどうするか?。じつは足すだけなんですよ(^^)。いちおう根拠はあります。
・2つ以上の微小入力がある時、その出力はそれぞれの微小入力に比例した結果を足すだけ。
上記は多変数関数の全微分に対応する意味になります(私見では)。以上が純粋な梁理論の全てです。厳密には純曲げ状態にある線材を梁(Beam)といい、純引張(圧縮)状態にある線材を棒(Bar)といいます。
・・・ええとここで、せんだん力Sはどうするの?って思いません?。じつはせんだん力は、梁理論では扱えないのです。せんだん力は図-5に示すように、図-6のy方向の力です。ところが梁理論ではまっさきに、y方向の変形は無視すると決めたのでした(^^;)。しかし変形はなくとも、せん断応力はあります。それがないと物体は釣り合えないからです。それで、
で・・・良いんじゃねっ?・・・となりました。これを平均せん断応力といいます。この発想は軸応力と同じです。軸応力もふと気づけば、「平均」軸応力ですよねぇ~(^^;)。曲げ応力についてだけ少々詳しく評価するのは、水平荷重の効果は鉛直荷重の効果に対してかなり小さく(構造系が横に長いから)、軸応力は曲げ応力よりだいぶ小さいのが現実だからです。またx方向のバネは、y方向のバネより圧倒的に弱いからです。このような前提が成り立つのは、D/L=1/5~1/10以下と言われています。
梁の微分方程式にいたる長い道程 その1 [ddt³さんの部屋]
梁の微分方程式にいたる長い道程 その1
d2y/dx2=M/(EI)のEIを移項した、
を、梁の曲げの微分方程式と言います。ネコ先生の定義、
https://nekodamashi-math.blog.so-net.ne.jp/2019-06-11-9
とは符号が違いますが、符号は座標系yの取り方や、Mの正方向の決め方に依存しますので、ここではこれで通します。でも梁の曲げの微分方程式とは、けっこう長い名前ですよね?。「の」が2回も出てくるし。
ところで梁の曲げの微分方程式にいたる道程は、けっこう長いんです。いくつかの別系統の概念や考えを総合しないといけないからです。
1.微小変形理論と弾性学
構造はある材料で造られます。経験的にいうと構造が目に見えるくらい変形したら、構造を構成する材料はぶっ壊れてるのが普通です。設計計算の目的は、材料がぶっ壊れない事の確認ですから、基本は目に見えない程度の微小な変形範囲の材料挙動という事になります。
材料にある変形が生じれば、変形に応じて材料内部には材料にストレスをかける力が生じます。これを応力といいますが、変形量から応力への関数を漠然と想像して下さい。変形量が入力で応力が関数の出力です。全ての物理現象には、暗黙の前提が潜んでいます。
・微小な入力に対しては、出力は微小入力に比例する。
これが微分可能という事の本当の意味です(個人的見解では(^^;))。全ての物理現象は微分可能なんです。
微分可能であれば当然連続です。全ての物理現象は連続でもあると思って、間違いないと思います。そうでない時には、必ず人為的な人間の手が加わっています。
物理現象が微分可能であれば、どんな材料についてもフックの法則のタイプが、微小変形には成立するはずです。
kは材料のバネ定数、uは材料の伸び(縮み)、Fはその時に材料に作用する力です。(2)の事を一般には材料構成則といい、(2)のタイプの材料構成則を扱う分野を、弾性学と言います。要するに弾性学とは、物体をバネの塊とみなす分野です。微小変形理論という用語は、弾性とは別の概念を表すために出来たものなのですが、実際上、微小変形理論と弾性学の成立範囲はほとんど重なります。
2.弾性係数,応力,歪み
材料構成則があれば材料挙動がわかるので、ある材料でつくられた任意の物体の状態は、原理的には計算可能なはずです。ところがフックの法則(2)には大きな問題があります。じつは(2)は、材料構成則ではないのです。
フックの法則はフックさんが、最も単純な構造部材として棒の伸び(縮み)をとりあげ定式化したものです。
図-3は、図-2のバネを長さがLになるようにm本連結し、断面積がAになるようにn本束ねたものだとみなせます。直列バネ,並列バネの関係より、このときバネ定数はkはk×n/mに変化します。まずバネ定数は材料定数ではなかったのです。材料定数は材質だけで決まるべきものですが、バネ定数は長さと断面積という物体寸法の影響を受けてます。これでは駄目です。これでは同じ材料、同じuとFであったとしても、物体形状ごとに材料状態が違うので、物体ごとに違った材料構成則が必要になります。そんなの、材料構成則ではありません。
フックに続いてヤングさんが現れます。ヤングさんは、自分では意識せずに人を小ばかにした態度をとってしまう、非常に腹立つタイプの天才だったらしいのですが、天才らしくコロンブスの卵に気づきます。
m本直列に連結すればバネ定数は1/m倍になり、n本並列に束ねればn倍になるなら、「単位長さ,単位断面積のバネ」でバネ定数を定義すれば良いと。それを弾性係数といい、記号Eで表します。
「単位長さ,単位断面積のバネ」の事をここでは「標準試験体」と呼んでおきます。長さLのバネは標準試験体をL本連結したものであり、断面積Aのバネは標準試験体をA本束ねたものですから、図-2が標準試験体でk=Eなら、図-3の棒のバネ定数は、
とすぐにわかります。図-3の棒に関するフックの法則は(2)より、
です。(4)をみると、長さLと断面積Aが明示されるので、物体寸法(物体形状)を後付けで考慮できる形になっており、Eは材質だけで決まる定数です。すなわちEは材料定数です。
しかし欲しいのは材料構成則です。図-3の棒のフックの法則を(4)のように表しても、内容はなんら(2)と変わりません。AもLもuもFも、物体寸法の影響をモロ受けだからです。要するに欲しいのは「標準試験体」に関するフックの法則なんですよ。手掛かりは、長さや断面積がどうあろうとバネはバネ。それと、直列バネ/並列バネの結果です。
バネはバネなので、標準試験体には(2)のkをEとした形が成り立つはずです。直列バネの結果より、標準試験体の伸びは長さLのバネの伸びの1/Lです。並列バネの結果より、標準試験体が分担する力は断面積Aのバネに作用する力の1/Aです。よって、(4)から簡単に誘導できる次の形が重要なのだと気づきます。
(5)は(2)のkをEに変えた形をしています。u/Lは長さLのバネの伸びを長さに応じて修正した、標準試験体の伸びです。F/Aは断面積Aのバネへの作用力を断面積に応じて修正した、標準試験体の分担力です。これって図-3に示した、棒の中にある標準試験体の状態を表してるのは明らかですよね?。これらは材質だけで決まります。
(5)が材質だけで決まる材料構成則である事を明示するために、u/Lをεで表し歪みと呼びます。F/Aもσで表し応力と呼びます。(5)から物体寸法に関する表示を追放できました。
(6)は、材料構成則である事を強調して特に応力-歪み関係と呼ばれます。歪みεの単位は、それが単位長さ当たりの伸びなので、[ε]=m/m=1、すなわち無単位です。応力σは単位断面積当たりの力なので、単位は[σ]=N/m2となり圧力といっしょです。よって弾性係数Eの単位も、[E]=N/m2になります。
(6)は図-3のような理想化された一様引張(圧縮)のモデルで得られたものです。現実の材料では一様な状態ではなく、場所ごとにεやσは違うでしょう。しかしEが材料定数であるために、歪みは変形勾配として、応力は力の密度として扱えば良いことになります。これらは物体形状とは無関係に定義可能な量です。
ちなみにコロンブスの卵に気づいたヤングさんを記念して、弾性係数をヤング率とも呼びます。弾性係数は、単位長さ,単位断面積当たりのバネ定数なので、まさに「率」なんですよ(^^)。応力と歪みを最初に実質的に導入したのもヤングさんです。
3.構造力学
構造力学は設計計算の実用理論です。主に設計計算を手計算で行う際に広く用いられます。構造力学の手順は以下です。
1) 荷重情報を用い、全体系の釣り合い条件から支点反力を算定する。
2) 荷重情報と算出した支点反力を用い、部分系の釣り合い条件から断面力を算定する。
3) 算出した断面力から断面に作用する断面応力を算定する。
4) 算出した断面応力を材料強度と比較する。
材料強度を測定するのが、いわゆる標準試験体に対する材料試験です。強度はもちろん応力で表されます(材質のみで決まるから)。ただ材料強度は設計基準といった人為的な「決め」の側面も大きいので、ここでは省略します。もちろん構造力学はこれだけではありませんが、最も多用されるのは、1)2)3)4)です。
全体系の例は、図-4です。
図-4は、両端に支点を持つ橋長Lの橋の模式図です。荷重P1~P3は、例えば橋を走行する車両重量で設計条件として与えられます。本当はこれに橋の自重が加わるのですが、面倒くさいので省略します(^^;)。橋の自重も、橋梁形式やコンクリート橋だとか鋼橋だとかは事前に決まっているのが普通なので、設計条件の一部です。R1とR2は支点反力と呼ばれ、荷重を支えるために支点で発生した力ですが、これらは未知数です。支点反力を算定するために、全体系の釣り合い条件を使います。
なぜ支点反力を算定するかというと、物体(橋)に作用する力を全て決定できなければ、内部の応力状態なんてわかりっこないからです。なぜ釣り合いかというと、構造物は静止してるからです。うろうろ歩く橋なんて見た事ないですよね?。間違ってそうなってたら、危険で近寄る事すらできません(^^;)。
作用力の合力が0と、回転力の合力が0となる条件は、明らかに以下です。
ただし鉛直上向きの力を正とし、回転力は左回りを正として左端の支点を回転中心としました。静止物体で回転力の釣り合いを考える際には、回転中心はどこにとっても良いので計算に便利なように選べます。また(7)下段の[力]×[回転中心からの距離]という量が、回転力だというのはご存じと思います。機械系ではトルクと言われますが、旧態然とした土木では「力のモーメント」と呼ばれます。L1,L2,L3は荷重P1~P3の載荷位置で、左端から測った距離であり、これらも設計条件の一部です。(7)を連立方程式として解けば、
が得られ、支点反力が算定されます。以後、支点反力は算定されたものとして、R1とR2を直接使用します。
断面力を算定します。そのために左端から距離xで構造を2つに分割して考え、部分系の釣り合いをとります。全体系が静止してる以上、部分系も当然静止してなけりゃいけません。
回転力M(x)も同様です。さっきと同じ理由でM(x)がなければP1の作用により、左部分系は左端の支点を回転中心として右回りに回転するはずですが、静止する必要があります。回転しないように回転力M(x)が働きます。断面に回転力が働くというのには、違和感を持つ人もいると思いますが、物体の釣り合い条件からの論理的な帰結です。なぜなら他に、回転力を働かせる場所がないからです。
N(x)については全体系に水平荷重は働いていないので、N(x)=0です。ただし一般的にはN(x)も無視できません。
M(x),S(x),N(x)を、曲げモーメント,せんだん力,軸力と言います。力のモーメントであるM(x)に、なぜ「曲げ」という冠がつくかは後でわかります。これらの力の由来は、もちろん図-5の右部分系です。左部分系は右部分系に支えられている訳です。一方右部分系には、左部分系と逆向きの大きさの等しい断面力が作用します。作用・反作用の法則です。つまり左部分系と右部分系は互いに支えあっています。これらの事実は、断面力が物体内部の内力であり、作用・反作用の法則は内力の伝達を表している事を示します。
やる事はいっしょです。左部分系の釣り合い条件は、
∴
と、位置xでの断面力が得られます。いまxはP1とP2の中間にとっていますが、xをP2とP3の中間にとったとすると、(8)にP2の効果も含めて同じことをするだけです。支点反力さえ算定できれば、任意の位置xでの断面力を算出できるのは明らかと思います。以後、断面力は算定されたものとし、M,S,Nを直接使用します。
(執筆:ddt³さん)