あらきけいすけのメモ帳

あらきけいすけの雑記帳2

壁に耳なし芳一、あるいは阪大の出題ミスを高校で習う範囲の数式できちんと説明する試み

例の入試ミス関連で、いくつかのブログやツイッターやPDFを見たけれど、高校数学の範囲内の数式を使って計算をきっちりと書いたものに出会えなかったので書いた。くだんのSEGの先生のメール*1やPDF*2を見たが、いかにも「受験物理のお行儀の良い答案」みたいな記述で、これでは天下の阪大の大先生たちを説得できはしないだろうと思う。たぶん「12月により詳細な指摘」(朝日新聞,1月6日)で数式を使った議論が出て、ようやくねじ伏せられたのだろう。理系の研究者は日頃から「論文の査読」というシステムでバトルしているので、そうそう簡単に説き伏せることはできないものである。

鉄1kgと綿1kg重いのは…

阪大の出題ミスで「阪大の言い訳」を朝日新聞で読んだとき、こんなクイズを思い出した「鉄1kgと綿1kg、重いのはどっち?」。こいつに「どちらも同じ」と答えたときに、「大気から受ける浮力の分だけ綿の方が軽い(kgって質量の単位だよね)」と返されたような気分になった。*3どこまでを「前提」としてるんだ?朝日新聞の図を引用しよう。

f:id:arakik10:20180114142242j:plain

この「音叉の問題」の全体を見ると、しょっぱなの問題がこれである

f:id:arakik10:20180114143958j:plain

きっちりとAモードの振動で疎密波の方向分布を問うている。ここから始めておいて、最後の問題だけ「Bモードもあるよ」ってのはちゃぶ台返しである。*4

f:id:arakik10:20180115041545p:plain

Google検索の結果(2018年1月15日アクセス)

愚痴、あるいはなぜ阪大ははまったのか?

出題者の学部学科なんて絶対に表に出ない。でも、これだけは言えるが、実は日本の大学の「理学部・物理学科」で流体力学を研究してる先生はほとんどいないのだ。いまの「物理学科」は量子力学が大好きで、流体力学を本格的に学びたいならば「宇宙物理学」「地球物理学」「機械工学」「航空工学」「核融合」といった学科、分野に進む方がいい(それとも東大理、京大理をめざす?)。今回も流体力学に慣れていないのでドツボにはまったのでは、と邪推している(てゆーか、受験物理で音波の反射は鬼門だという意見もツイッターとかで見た)。

波ってどのように式で表すのだろう?

まずは次のグラフを見て、これが「右に進む波」なのか「左に進む波」なのか分かるだろうか?

f:id:arakik10:20180114133912j:plain

正解は「わからない」である。波の山や谷が次の瞬間にどちらに動くかという情報が無いので何とも言えない。ではどちらに動くかということを数式ではどのように表すのだろう?答えは音源の位置をx=0として次のようになる:

 右(x>0の方向)に進む波:\sin\left[2\pi f\left(t-\frac{x}{c}\right)\right],

 左(x<0の方向)に進む波:\sin\left[2\pi f\left(t+\frac{x}{c}\right)\right]

ここでfは音源の周波数、cは音速。この式の意味は、波の山がxの位置にあったとすると、音源がその山を作ったのは|x|/c秒前だということ。波を表す関数は時刻と位置の関数だから「2変数関数」という高校数学の範囲外の関数が現れる。でも、どちらも三角関数の引数の中にお行儀良くあらわれているので、慣れれば大したことないと思う。

波の伝搬と反射の設定と計算

 まず音叉から出る右(正方向)と左(負方向)に進む音波をA-1(問1)の状況をきっちりと考慮に入れ、

f:id:arakik10:20180115060712p:plain

密度の分布が音叉に対して左右対称になるように、速度を次式で与える(振幅を1と置き、減衰はしないとした)*5

 u_{右}(x,t)=\sin\left[2\pi f\left(t-\frac{x}{c}\right)\right] (ただしx>0),

 u_{左}(x,t)=-\sin\left[2\pi f\left(t+\frac{x}{c}\right)\right] (ただしx<0),

ここでfは音源の周波数、cは音速。ここで音が反射する壁での状況を考える。壁(x=-d)での「壁への入射波」u_{左}

 u_{左}(-d,t)=-\sin\left[2\pi f\left(t-\frac{d}{c}\right)\right]

となる。ここで「壁からの反射波(右に進む波)」u_{反}

 u_{反}(x,t)=\sin\left[2\pi f\left(t-\frac{x}{c}\right)+A\right]

と置く(音速、周波数は反射で変わらないことを使った)。反射波の振動は山と谷がどこに来るかわからないので、位相にAを入れておく。ここで未知数Aは、壁(x=-d)での条件である「空気は壁を突き抜けて進むことはない」ことより

 u_{左}(-d,t)+u_{反}(-d,t)=0

が常に成り立つこととなるから*6

 2\cos\left(2\pi ft+\frac{A}{2}\right)\sin\left(2\pi f\frac{d}{c}+\frac{A}{2}\right)=0

 \Longrightarrow 2\pi f\frac{d}{c}+\frac{A}{2}=n\pi,

ここでnは整数。これより反射波を表す関数は

 u_{反}(x,t)=\sin\left[2\pi f\left(t-\frac{x}{c}-\frac{2d}{c}\right)+2n\pi\right]

となる。この反射波が音叉から右方向に出る音波と共鳴する条件は次のようになる:

 u_{右}(x,t)=u_{反}(x,t) (ただしx>0),

 \Longrightarrow \sin\left[2\pi f\left(t-\frac{x}{c}\right)\right]=\sin\left[2\pi f\left(t-\frac{x}{c}-\frac{2d}{c}\right)+2n\pi\right],

 \Longrightarrow -2\pi f\frac{2d}{c}+2n\pi=2k\pi,

ここでkは整数。 ここで波長\lambda\lambda=c/fで与えられることを考慮すると、求める条件は

 2d=m\lambda (ここでm=n-kは整数)

となる。

SEGの参考書は書店で拝見したことがあるが、SEGで学ぶ国公立理系上位校ねらいの生徒たちなら、これくらいの数式くらい簡単にフォローできるだろう。

入試出題採点経験者として言うが、出題者は学習指導要領に縛られるが、解答者は縛られる必要は全くないのである。

*1:https://www.dropbox.com/s/2vlzte1toqrgsf2/大阪大.pdf?dl=0

*2:https://www.dropbox.com/s/z5x1pq7k7i8daoy/osaka_u.pdf?dl=0

*3:こうゆうのをネットスラングで「小並感」というwwwww

*4:気づくのが遅すぎて、採点のやり直しもへったくれもない状況だったので、強弁せざるをえなかったんだろう。また、問5の設定から、問4の解答を逆算することもできるので、入試問題としてはいまいちかもしれない。

*5:出題者は左右の速度の設定でこけたのか?

それともう一つ。ここでの音叉から出る音は4重極放射なので、波の振幅が変わらないとした仮定はかなり非現実的かもしれない。

*6:これは学習指導要領の高校物理の範囲内では「固定端反射」に相当する。ちなみにここで壁は音が当たっても振動しないことになっているので、壁の向こう側に音は伝わらない。だからこの実験をレオパレスでやることはお勧めできない…ような気がする。

C言語の atan2(y,x) を arccos と sgn を使って書いてみた

授業用の覚書:関数電卓もそうだが、座標(a,b)偏角を求めるときに、\theta(b,a)=\arctan\frac{b}{a}では第3象限-\pi\lt\theta\leq-\frac\pi2, 第4象限\frac\pi2\leq\theta\leq\piの値を返さないので、C言語atan2(y,x) と同じ動作の関数を標準的な数学関数arccos, sgnの組み合わせで作ってみた。定義域は(a,b)\in\mathbb{R}^2\backslash\{(0,0)\}.

\displaystyle\theta(b,a):= \underbrace{{\rm{sgn}}(b)}_{(-\pi,0)に拡張}\overbrace{\arccos\frac{a}{\sqrt{a^2+b^2}}}^{[0,\pi]に値を持つ}+\overbrace{R(b,a)}^{0,\pi を担う}

\displaystyle R(b,a):=\frac{\pi}{2}(1-{\rm{sgn}}(a))\underbrace{|{\rm{sgn}}(a)|(1-|{\rm{sgn}}(b)|)}_{a\ne0,\ b=0\ を取り出す係数}

なぜ初等関数にないのだろう?

 

2023.7.4 desmosを使って実装してみた

www.desmos.com

C言語の仕様に「グローバル変数」という名前は出てこない

授業の資料作成で調べ物関連のメモ。C99の仕様書 "ISO/IEC 9899:1999 - Programming languages C" を参照していたら「グローバル変数」という語は全く無かった;その代わり "scope" (6.2.1), "linkage" (6.2.2) という語が丁寧に定義されている。いわゆる「グローバル変数」に相当する部分を拾い読みする:

If the declarator or type specifier that declares the identifier appears outside of any block or list of parameters, the identifier has file scope, which
terminates at the end of the translation unit. (6.2.1 Scopes of identifiers #4)

関数の定義のうち関数本体(function body)も波括弧で括られた複文 "block" なので、"any block" で「そのファイルのすべての関数の外」が含意されるようだ。

An identifier declared in different scopes or in the same scope more than once can be made to refer to the same object or function by a process called linkage. (6.2.2 Linkages of identifiers #1)

"linkage" とはメモリ上で同じ場所にあるものを読み(書き)できること。

If the declaration of a file scope identifier for an object or a function contains the storage-class specifier static, the identifier has internal linkage. (6.2.2 Linkages of identifiers #3)

変数のスコープがそのソースコードファイルの全体になる場合に、その変数は "internal linkage" を持つ。

定義の合わせ技を読み取るのは難しい。

 

Raspberry Pi 3 の SPI を用いてA/Dコンバーターからデータを読む python 関数 readadc() の各行の動作

教育用の覚え書き。ラズパイ3で標準のSPIバス(+RPi.GPIOライブラリ)を用いてA/DコンバーターICからセンサ値を読み出す解説記事に頻出する python 関数 readadc() の各行の動作を確認したのでメモを残す。output()でbit1個分のピンへのON/OFF信号の送出, input()でbit1個分のデータ読み出しを実行しているので、呼び出し時の関数の動作時間は多分、一定していない。コメントはA/DコンバーターとしてMCP3208を想定して書いている。

def readadc(adcnum, clockpin, mosipin, misopin, cspin):
    # 宣言しておく import 文
    #   import RPi.GPIO as GPIO
    # 入力(すべてint型)
    #   adcnum:   データを読み出すMCP3208のチャネル番号 0, ..., 7
    #   clockpin: クロック(CLK)信号を出力するラズパイのピン番号
    #   misopin:  Master In Slave Out のラズパイのピン番号
    #   mosipin:  Master Out Slave In のラズパイのピン番号
    #   cspin:    Channel Selector のラズパイのピン番号
    # 戻り値
    #   int 型, 0, 1, 2, ..., 4095: adcnum チャネルのA/D変換値
    #   -1: adcnum が 0, ..., 7 の範囲にない

    # 読み出しチャンネル番号が 0, ..., 7 に無ければ処理しない
    if adcnum > 7 or adcnum < 0:
        return -1

    # CLK への出力を0に初期化しておく
    #   この関数では CLK のパルス生成を GPIO.output(*,1), 
    #   GPIO.output(*,0) の連続呼出で実行している
    #   CLKパルスの立ち下がりのタイミングで1bitの入出力が起こる
    #
    GPIO.output(cspin, GPIO.HIGH)   # I/Oが起きないようにCSを1にする
    GPIO.output(clockpin, GPIO.LOW) # CLK を0に設定しておく
    GPIO.output(cspin, GPIO.LOW)    # I/Oを始めるためにCSを0にする

    # この3行は commandout に MSB IN のデータ(8bitで上位5bitがデータ)を
    # 作成する部分(秋月の MCP3204/3208 のPDFの p.20, Fig. 5-1.)
    #
    # タイミングチャートの3行目の D_IN に書いてある
    #
    #  [Start][SGL/DIFF][D2][D1][D0]
    #     |       |      |   |   |
    #     |       |      読み出しチャンネルの指定
    #     |   0:差動入力
    #     |   1:シングル・エンド
    #  1:スタートビット
    #
    # を作成する
    commandout = adcnum # [ 0][ 0][ 0][ 0][ 0][D2][D1][D0] adcnum
    commandout |= 0x18  # [ 0][ 0][ 0][ 1][ 1][D2][D1][D0] 0x18と論理和
    commandout <<= 3    # [ 1][ 1][D2][D1][D0][ 0][ 0][ 0] 3bitシフト

    # このループは作成した commandout の上位5bitを1bitずつ送り込むループ
    for i in range(5):

        # commandout の最上位1bitが1ならば、MOSIに1を送り込み、
        # さもなければ、MOSIに0を送り込む
        if commandout & 0x80:
            GPIO.output(mosipin, GPIO.HIGH)
        else:
            GPIO.output(mosipin, GPIO.LOW)

        # MOSIに送り込んだら commandout の内容を1bitずらす
        commandout <<= 1

        # CLKに1パルス送り込む
        GPIO.output(clockpin, GPIO.HIGH)
        GPIO.output(clockpin, GPIO.LOW)

    # 出力データ adcout を1bitずつ読み込むループ
    adcout = 0          # データ保存変数 adcout の初期化
    for i in range(13): # 全部で 13 bit 読み込む

        # A/Dコンバーターの CLK に1パルス送り込む
        GPIO.output(clockpin, GPIO.HIGH)
        GPIO.output(clockpin, GPIO.LOW)

        # adcoutを1bit左にずらして、最下位bitに空きを作る
        adcout <<= 1

        # もし MISO に1が入っていたら、adcout の最下位bitに1を書き込む
        if i>0 and GPIO.input(misopin)==GPIO.HIGH:
            adcout |= 0x1

    # データ通信が終わったら CS をHIGHにして通信を止める
    GPIO.output(cspin, GPIO.HIGH)

    # 読み出しデータを呼び出し側に戻す
    return adcout

複数のGPIOピンからのpwm信号出力が欲しいがゆえにpigpioライブラリをインストールしたが、ライブラリ中のSPIの関数pigpio.spi_*()たちの引数の与え方が分からずに*1、このコードをpigpio.read()/write()を使って書き換えた関数を使ってセンサ値を読み出しているのは内緒。

非同次1階線形常微分方程式の初期値問題の解法

授業のための覚書。y, p(x), q(x)xの関数とするときの\displaystyle\frac{dy}{dx}+p(x)\,y=q(x)の解を求める。

 

大学教養程度の薄い教科書では不定積分を使うものしか見かけないのでメモ。

 

まず同次方程式\displaystyle\frac{dy}{dx}+p(x)\,y=0を変数分離法でx=0からx=tまで積分すると、\displaystyle\int_0^t\frac{1}{y}\frac{dy}{dx}dx=-\int_0^tp(x)dxである。ここで部分積分法より左辺の積分変数をyに変数変換すると\displaystyle\int_{y(0)}^{y(t)}\frac{1}{y}dy=-\int_0^tp(x)dxとなる。よって積分の結果は\displaystyle\ln|y(t)|-\ln|y(0)|=-P(t)となる(ここで\displaystyle P(t):=\int_0^tp(x)dxとする)。これより\displaystyle\left|\frac{y(t)}{y(0)}\right|=e^{-P(t)} \Longrightarrow \displaystyle\frac{y(t)}{y(0)}=\pm e^{-P(t)} だが、t=0の条件より負号は棄却され、\displaystyle{y(t)}={y(0)}e^{-P(t)} となる。*1

次に非同次の解だが、まずy(x)e^{P(x)}微分すると\displaystyle\left(y(x)e^{P(x)}\right)' \displaystyle=y'(x)e^{P(x)}+y(x)P'(x)e^{P(x)} \displaystyle=y'(x)e^{P(x)}+y(x)p(x)e^{P(x)} \displaystyle=q(x)e^{P(x)}であるから、この最左辺と最右辺をx=0からx=tまで積分すると、\displaystyle y(t)e^{P(t)}-y(0)e^{P(0)}=Q(t)となる(ここで\displaystyle Q(t):=\int_0^tq(x)e^{P(x)}dxとする)。P(0)=0に留意して整理すると、求める解は次のようになる:\displaystyle y(t)=y(0)e^{-P(t)}+e^{-P(t)}Q(t).

 

初期値問題の解のx=tでの値:\displaystyle y(t)=y(0)e^{-P(t)}+e^{-P(t)}\int_0^tq(x)e^{P(x)}dx, ただし \displaystyle P(t):=\int_0^tp(x)dx

*1:2019.4.19改稿; 旧版:y(t)y(0)から連続につながっているので、y(0)と同符号になるから、左辺は\displaystyle\ln|y(t)|-\ln|y(0)|=\ln\frac{|y(t)|}{|y(0)|}=\ln\frac{y(t)}{y(0)}となる。以上より両辺の\lnをはらって、y(t)=y(0)e^{-P(t)}となる。

e関連の極限の式を避けて対数関数と指数関数を微分する

 授業のための覚書。高校の数学IIIの合成関数の微分の知識を用いる。数学IIIの知識が一通りあれば、論理で押し切れそう。

 対数関数は対数法則 f(xy)=f(x)+f(y) を満たす可微分な関数なので、合成関数の導関数の公式を使い、この両辺を y微分して xf'(xy)=f'(y) となる。これに y=1 を代入して xf'(x)=f'(1) となる。ゆえに \displaystyle f'(x)=\frac{f'(1)}{x} である。ここで f'(1)=1 となる特別な対数関数を「自然対数」と呼び、 \ln x と表記することにすると、自然対数の導関数\displaystyle (\ln x)'=\frac{1}{x} となる(ところで対数の底の値は?)。

 自然対数の逆関数となる指数関数*1y=f^{-1}(x))の導関数\ln y=x の両辺を x について微分することで y'=y と分かる。この指数関数を便宜的に f^{-1}(x)=e^x と表記することにすると、(e^x)'=e^x となる。

 対数法則に x=y=1 を代入して f(1)=0 となること、x\gt0 のとき x^{-1}\gt0 を使って、対数関数の x での値は定積分 \displaystyle \ln x=\int_1^x\frac{dt}{t} で定義される数である。この定義が破綻しないためには \ln x の定義域を x\gt0 にしないといけない。

 対数の底の具体的な値は \displaystyle\int_1^e\frac{dx}{x}=1 を満たす1より大きい数 e数値計算で求めればよい(よい精度で求められるアルゴリズムがあればよいのだが、ボクは知らない)。

 自然対数以外の底の対数関数は\displaystyle f(x)=f'(1)\int_1^x\frac{dt}{t}=f'(1)\ln xで与えらえる。ここでこの対数関数の底をaとするとf(a)=f'(1)\ln a=1となるので、\displaystyle f'(1)=\frac{1}{\ln a}となる。これよりこの対数関数を\log_aと書くと、\displaystyle\log_ax=\frac{\ln x}{\ln a}となる。これは底の変換の式になっている。これより\displaystyle(\log_ax)'=\frac{1}{x\ln a}となる。

*1:f(x)=X, f(y)=Y と置くと、対数法則の式を指数法則の式 f^{-1}(X)f^{-1}(Y)=f^{-1}(X+Y) に書き換えることができる。