2007年12月30日

三角関数の微分

この解法はうまくて間違いがないように見えるけれども,どうしたらそれを思いつくことができるだろうか.    この実験はうまくて事実を示すように思われるが,どうしたらそれが発見できたであろうか.    どうしたら私は自分でそれを思いついたり発見したりできるであろうか.

さあ一人ぼっちの年末なので三角関数の微分(といってもsineとcosineだけなんだけど)でもしよう。

で、そのために加法定理

20071230eq1

と次の極限の式

20071230eq2

の2つを途中で使うのだけれど、証明がよく分からないめんどくさいので天下り的に使う。

加法定理はお絵かきしてもらうかオイラーの公式使うかすればすぐ証明できるのでいいとして(いいんだろうか)、 2つ目が若干分かりにくい。単位円の接線は角度0付近でy軸と平行なので、角度0付近ではsin xの変化と角度xの変化量が等しくなる(=両者の比が1になる)、と言えば想像しやすいだろうか。普通は「はさみうちの原理」 を使って証明する、はず、多分、ググれ。

じゃあとりあえずsinの微分から。まず微分の定義は

20071230eq3

だ。xをちっと増やしたときの式の値の変化量。まずはこの中にそのままsin xを突っ込む。

20071230eq4

加法定理を使って

20071230eq5

右辺sin xの項をまとめ、ついでに分母を分けて整理すると

20071230eq6

ここで最初に提示したsin x/xに関する極限の式を思い出してもらえば、右辺第2項目がcos xとなることはすぐ分かるだろう。

問題は1項目だ。cos 0 = 1だから分子は0になるだろうが、分母も0になる。0/0、不定形。こいつは0に収束するはずだ、 ということは分かっていたけれども、どうも上手い変形が思いつかなかった。まあここは大人しく参考書に頼ろう。

とりあえず見やすくするために問題の1項目だけを取り出す。

20071230eq7

まず、分子分母に(cos h + 1)を掛ける。

20071230eq8

ここで分子の1という数字に注目する。単位円を考えると、斜辺は1、底辺はcosθ、高さはsinθなので、 三平方の定理よりsin2θ+cos2θ=1が導かれる。よって、 1=sin2 h + cos2 hという変換ができる。この変換を分母の1に対して施せば、

20071230eq9

さあ、やっとsin h/hが外に出てきた。極限をとろう。

まず、sin h/hが1となることは最初に示した極限の式より。

次に分母のsin hは0、分子のcos hは1となるので、答えは

20071230eq10

0だ。やっと0になった。しかし上手いやり方だ。特別難しいことはしていないのに差が積になってしまった。 自分で思いつけたら楽しかったろうなぁ。

さあ最初の式に戻ろう。

20071230eq6

ここで、最初に示した式極限のより右辺第2項がcos x となることは明らかであり、第1項が0となることは今示した。 よって

20071230eq11

sineの微分はcosine!

さて、続けてcosineの微分。微分の定義の式

20071230eq3

へcos xを突っ込めば

20071230eq12

加法定理より

20071230eq13

cos x の項をまとめ、ついでに分母を分ければ

20071230eq14

ここで右辺第1項目は先ほどと同様の変形を施すことにより、0に収束することが分かる。第2項目についてはsin h/hが1へ収束するので、

20071230eq15.

cosineの微分は-sine!

つーことはもっかい微分すると-cos xになり、さらにもっかい微分するとsin xになる。要するに4階微分で元に戻るわけだ。

三角関数が四つの式を回るわけですよ。

参考

posted by Rion778 at 23:50 | Comment(4) | TrackBack(0) | 勉強ノート | このブログの読者になる | 更新情報をチェックする
2007年12月29日

半分はやさしさでできています

つまりやさしさだけじゃ何も出来ないということだな。

微妙に風邪ひいて頭と喉が痛いので風邪薬を探すも、バファリンしか入っていない薬箱。

風邪を引けばだるくて薬屋には行かないし、元気なときには面倒だから買い置きなんかしない。 いざというときの自分に対するやさしさが足りない。まあ頭痛は治まった気がするからいいや。

多少元気なうちに実家に帰ろうかとも思ったけど、実家はすごく寒い(暖房器具の質的な意味で)ことをかろうじて思い出す。やっぱやめ。 もう今年の残りは寝て過ごそう。

posted by Rion778 at 19:24 | Comment(0) | TrackBack(0) | diary | このブログの読者になる | 更新情報をチェックする

普通が一番ですよ

外食ばっかしてたら卵の消費期限がヤバイことになってたのでこんな時間に玉子焼き。

作ってる間の数十秒ヒマなのでパッケージを読んでたら

「とうもろこし、よもぎと、木酢液といった大地の恵みで云々かんぬん」

みたいな文言が。

トウモロコシはいいとして、ヨモギ?

と思ったけど、ぐーぐる先生に聞いてみたらヨモギは結構飼料として使われてるようだ。

まあそこまではいいとして、木酢液?

これもぐーぐる先生に聞いてみると、どうやら添加試験は結構されてるらしい。 斜め読みしたところでは卵のにおいがなくなるとかビタミン含量が増えるとかなんとか。なので飼料としてもそれなりに使われているようだ。

まあとにかく良い効果があるのは分かった。

でも前面に押し出していいの?木酢液。アルカロイドやら発ガン性物質やらの塊みたいなもんじゃないのか。昨今は「植物由来」 だからって無条件にホイホイと有り難がるような人間も減ってきただろうに(多分)。黙って使えばいいのに…

というか待て。

トウモロコシもヨモギも木酢液も飼料として珍しくないのだから、要するにそれって普通の卵って事じゃ

posted by Rion778 at 01:26 | Comment(0) | TrackBack(0) | diary | このブログの読者になる | 更新情報をチェックする
2007年12月26日

電話をかけて

プププっと聞こえたら

「ソフトバンクか…」

というサービス

どうなんだそれ。

posted by Rion778 at 21:45 | Comment(0) | TrackBack(0) | diary | このブログの読者になる | 更新情報をチェックする
2007年12月25日

本場では今日から年明けまでクリスマスなわけですが

これまでの人生において他に例を見ない高さから生卵落っことした。

なぜキャッチせずはたき上げてしまったのか。悔やまれるところです。

posted by Rion778 at 19:22 | Comment(0) | TrackBack(0) | diary | このブログの読者になる | 更新情報をチェックする

溶液を混合する際の濃度誤差について

2種の溶液を混合するとしよう。

この2種の溶液を溶液A、溶液Bと名づけ、混合する際の体積をそれぞれa、b と表記する。

このとき、混合液中の溶液Aの濃度を体積比で表せば次のようになる。

20071225eq1

これを目的の濃度とする。

だが、完全な精度で溶液を量りとることは不可能である。すなわち、a、bには常にいくらかの誤差が含まれることになる。

誤差を例えばεで表せば、

20071225eq2

が実際の溶液Aの濃度ということになる。

ここで誤差の大きさの程度を表すような指数を作るため、両者の比を取ることにする。

20071225eq3

つまり

20071225eq4

という式によって表される値を誤差の程度を表す指標としよう。式を見れば分かるように、 この値は目的の溶液Aの濃度に対し実際の溶液Aの濃度が「何倍に」なっているのかを表し、εa、 εbがそれぞれ0のときに1となる。

この式はεa、εbの関数と見ることもできるので、

20071225eq5

と表そう。ここで知りたいのは、誤差εの変化がこの指数にどの程度の影響を与えるのかという情報だ。

しかし、この式は分子に変数があるため変化量を正確に計算するのは結構面倒くさい。

が、そもそもそれほど厳密にこの指数の変化量を知る必要は無いだろう。

なぜなら対象としているのは誤差なのだから、εの値はいずれもa、bに比べてかなり小さい値だろう(というかそうじゃないといけない) 。なので、単位を上手く設定すればεはほとんど0の周囲となるはずだ。

0の周囲といえば…0の周囲でのテイラー展開、 マクローリン展開!

ところで、マクローリン展開するには上の式を微分しないといけないのだが、そのまま馬鹿正直に偏微分すると結構面倒なことになる。

20071225eq6

20071225eq7

別にこいつを使って計算を続けてもいいのだけど、ちょっとごちゃごちゃしすぎている。

そこで、 εaを考えるときはεbを0に、 εbを考えるときはεaを0に固定してしまおう。それから偏微分しよう。 εは0に近いはずだから多分問題はない。多分。あったらどうしよう…

まあとにかくそう考えることにすると、誤差の指標の式

20071225eq4

はεaの変化による部分、

20071225eq8

εbの変化による部分、

20071225eq9

とに分割できる。

で、計算はなるべく簡単にしたいし、0の付近を見るだけだから1次までで十分だと思う。 なのでそれぞれ一次までマクローリン展開して近似すると

20071225eq10

20071225eq11

この式から 、εaの影響はa < bの時大きく、a > bの時小さい、 εbの影響は溶液の総和が大きいほど小さいということが分かる。まあ考えれば当たり前のことなんだけど。

それで、例えば溶液Aの濃度に対して5%以内の精度が欲しいのであれば、 上式において第2項目が±0.05以下になるようなεの範囲を目安に、それよりも誤差が小さい計量器具で溶液を量りとればいい。

要は、「比率の小さい方の溶液は正確に量りとれ」、「少量を正確に量りとれないなら多目の量で混合しろ」と。

ちなみに、希釈とかで2種の溶液体積の比率が非常に大きく、なおかつ体積が少ない方の溶液の濃度に注目するときは、 こっちの誤差

20071225eq8

を計算するときに分母の誤差の影響はほとんどないとみなし、

20071225eq12

を用いても大した問題はない。少しだけ溶液Aの誤差の影響が少しだけ過大評価されるが、大きく見積もられる分には問題ないだろう。

うーん。しかしこう考えていくと学生実験はスゲーいい加減だったなぁと思ったり思わなかったり。

posted by Rion778 at 01:53 | Comment(0) | TrackBack(0) | 勉強ノート | このブログの読者になる | 更新情報をチェックする
2007年12月24日

フーリエ級数展開

イブなのでフーリエ級数展開の復習でもしよう。

何?神?いるに決まってんだろ。godなんだから。

…えーと、

テイラー展開においては、 級数展開された式

071221

のxへ0を代入→両辺を微分→xへ0を代入→両辺を微分…

という操作を繰り返すことにより、係数aを順次求めていった。

ところで、目的とする関数が周期関数であったなら、 同じように周期関数であるsineやcosineの級数として展開したいと思うところ(とりあえず周期は0から2π、 もしくは-πからπまでの周期関数を考える)。

20071224eq1

こんな感じに展開したいわけだ(もちろんcos0=1、 sin0=0なのでa0の項はa0になるし、b0の項は0になる)。 この形の展開をフーリエ級数展開と呼ぶ。狽使って書けば、

20071224eq6

となる。

テイラー展開では微分して0を代入すれば余計な項が消せた。

では、フーリエ級数展開において「微分して0を代入」に相当する操作は何なのか。

それは「cos nx (またはsin nx)を掛けて0から2πまで(もしくは-πからπまで)積分する」という操作だ。

この操作を行うと、cos nx (またはsin nx)以外の項がすべて0になるので、 テイラー展開のときの様に目的の項の係数を計算できる。

とりあえずはsin(x)とcos(x)についてそれぞれの0から2πまでの積分をグラフで見てみよう。

20071224graph1

20071224graph2

積分は、グレーで示した部分の面積の総和に相当する。いずれについても0>yの部分と0<yの部分の面積が等しいので、 計算結果が0になることが分かるだろう。

次に「cos nx (またはsin nx)を掛けて0から2πまで(もしくは-πからπまで)積分する」 のグラフをいくつか見てみよう。

…っとその前にグラフを書くのに使ったRのスクリプトの雛形はこんなの。 sineとcosineの違う組み合わせが見たかったら1行目をいじると見られる。

f <- function(x) sin(1*x)*sin(1*x)
plot(f, 0, 2*pi,ylim=c(-1,1))
xvals <- seq(0, 2*pi, length=100)
dvals <- f(xvals)
polygon(c(xvals,rev(xvals)),
 c(rep(0,100),rev(dvals)),col="gray")
grid()

それじゃグラフを

20071224graph3

20071224graph8

20071224graph4

20071224graph5

20071224graph6

20071224graph7

これらを見ると分かるように同じものを掛けて積分したときだけ値がゼロではない値となる。ちなみにπだ。 証明は良く分からないので面倒なので省略。要は直交しているベクトルの内積がゼロになるのと同じことだ。たしか。 関数が直交してるんだ。関数が。

…えーとつまり、両辺にcos nxを掛けて0から2πまで積分すると係数anにπの掛かったものが残り、 両辺にsin nxを掛けて0から2πまで積分すると係数bnにπの掛かったものが残るということだ。

式で書けば

20071224eq2

20071224eq3

係数を求める形にすれば、

20071224eq4

20071224eq5

つまりこれが係数an、bnを求めるための関数、要は一般項。

あとは、この操作によって得られた係数an、bnを順次

20071224eq6

へ入れていけば良いのだけれども、a0の項とb0の項は注意しないといけない。

まず、b0の項が恒常的に0となるのは言わずもがな。sin0は0だから。

a0の項だけど、この項だけを残すために上の式にcos0、つまり1を掛けて両辺を0から2πまで積分すると、

 20071224eq7

となり、a0に2πの掛かった状態となる。a0を求める状態にすると

20071224eq8

となる。要するにさっきのanを求める式の半分だ。なので、さっきのanを求める式へn=0を入れると、 実際のa0の倍の値が帰ってきてしまう。

よって、先ほどの式を使ってan、bnを求めるということにこだわるのであれば、 a0は倍の値が帰ってきてしまうということ、b0は0であるということを意識して、

20071224eq6

を次のように書き直す。

20071224eq9

つまり、この式の中にあるa、bという係数を

20071224eq4

20071224eq5

という2つの式によって求め、代入していけばフーリエ級数展開の完成というわけだ。

ところで、この係数an、bnというのは基準となる周期 (ここでは0から2π)のn分の1の周期を持ったsine、cosineの「大きさ」に他ならない。

どうせなら2つまとめて「大きさ」を評価したい。

sineとcosineをまとめるといえば…「オイラーの公式」。

きっと次回はそのへんから。

オイラーの公式を知らない人は、cos(x)、sin(x)、 eixをそれぞれテイラー展開して比べてみると幸せになれるかもしれない。iはもちろん虚数単位のi。 二乗すると-1になるやつな。それじゃ今すぐヤレ。面倒とか言わず。今すぐだ。今すぐ。

posted by Rion778 at 00:10 | Comment(0) | TrackBack(0) | 勉強ノート | このブログの読者になる | 更新情報をチェックする
2007年12月23日

テイラー展開(2)

何らかの具体的な課題を頭におかずに方法を探す人は,たいてい失敗の憂き目に会う.

──ヒルベルト

引き続き気分が晴れないので指数関数でもテイラー展開しよう。

ふつう単に指数関数といったら

20071223eq1

を指すことが多い。この関数は何回微分しても形が変わらない。というかそういう風にeという数字が決まっている。 eは自然対数の底とかネイピア数とか言い、πと同じ超越数。

微分しても形が変わらないということの意味は、exという関数のどこで接線を取ってみても、 その接線はxが1進む間にf(x)分だけ大きくなるような傾きを持っているということ(参考:アイとイーとパイのお話私事歳時記@ はてな)。

とにかく、微分しても形が変わらないのだからテイラー展開はしやすいだろう。それに、 eは無理数だけれどxにゼロを代入すれば1になる。

そこで、x=0の周りでテイラー展開、(つまりマクローリン展開)してやると、

20071223eq2

狽使って書き直せば

20071223eq3

さて、以前書いたフーリエ級数展開式をプロットするスクリプトを流用してグラフを書いてみよう。 スクリプトの雛形は以下の通り。

yx <- function(n,x) x^n/factorial(n)
my.Taylor <- function(m,k,o=0){
f <- function(m, x, a, o){
res <- 0
for(n in 0:a) res <- res + m(n,x)
return( o + res )
}
y <-seq(-3, 3, by=4*pi/1000)
return(plot(y,f(m,y,k,o),type="l",xlab="x",ylab="y",
 xlim=c(-3,3), ylim=c(-3,10), col="red"))
}
my.Taylor(yx, 1)
par(new=T)
f <- function(x) 1/(1-x)
plot(f, -3, 3, ylim=c(-3,10), col="blue",ylab="y",main="n=1")
grid()
legend(
 -3,10,
 legend=c("exp(x)", "Taylor polynomial"),
 lty = c(1,1),
 col = c("blue", "red"))

まずは0次近似(x=0における切片)。

20071223graph1

次に1次近似(切片+x=0における傾き)。

20071223graph2

2次近似

20071223graph3

3次

20071223graph4

4次

20071223graph5

5 次

20071223graph6

これ以上はこの範囲で見ていても大差ないのでやめよう。見てのとおり高次の項まで使うほど精度が良くなる。

一方で、x=0の付近であれば3次程度の近似でも非常に精度が良いことが見て取れるだろう。

ちなみに、指数関数の場合は次数を無限大に増やせば多項式はx=0の付近以外のあらゆるxについても指数関数へ収束する。

「指数関数の場合は」とわざわざ断るからには、収束しないような関数もある。

例えば次の関数

20071223eq4

何のことは無いただの反比例の式。グラフは

20071223graph7

Rでプロットするとx=1のところに線が入ってしまうが、式を見ると分かるようにx=1は定義されない。 その上右から近づいたときと左から近づいたときで発散の仕方も異なる。

それで、この式はx=0の付近でテイラー展開すると

20071223eq5

とかなりシンプルな形になるが、そのグラフは

20071223graph8

20071223graph9

20071223graph10

 

20071223graph11

20071223graph12

と、|x|<1の範囲以外では何項使おうが収束しない。というか発散する。また、nが偶数か奇数かで発散の仕方も異なる。 このときの|x|<1を収束半径とか呼ぶ。収束半径の外ではテイラー級数は収束しない上、有限項まで取ってみても近似にすらならない。

他にはsinとかcosとかの収束の様子を観察すると面白いのだけれど、もうグラフ描くのがめんどいのでパス。

ニコニコが見られる人にはこういう動画も。

この先は教科書的にフーリエ級数展開、フーリエ変換へと進むべきだろうか。

posted by Rion778 at 19:23 | Comment(0) | TrackBack(0) | 勉強ノート | このブログの読者になる | 更新情報をチェックする
2007年12月22日

テイラー展開

御伽噺の中ではカメは勝った。

でも、ウサギが勤勉なウサギだったら?

もうなんか鬱なのでテイラー展開の復習でもしよう。

以前フーリエ級数展開なんてのを見た。 あれはサインとコサインを上手い事調整しながら足し合わせていくことで目的の関数と同じものを作り上げるというものだった。 目的の関数をサインとコサインにバラすと言ってもいい。

テイラー展開は目的の関数を多項式にバラす。

071221

a0, a1, a2 ...は適当な係数。 こいつらを上手く決めてやれば、きっと右側の式は左側の式と同じになる気がするような、しないような。

が、このままでは係数はどれ一つ決められない。そこで、例えばxに0を入れてみよう。もちろんf(0)の値は得られるものとして。

そうすると、最初の係数a0以外はxが掛かっているものだからみんな消えてしまう。よって

071221eq2

一つ目の係数の値が分かった。というかまあ切片が分かっただけなんだけど。ちなみに、別にx=0でなくともxが0に近いなら

071221eq3

ということは言っても良いだろう。まあ近似といえば近似。

しかしもう少し情報を足したい。傾きくらいまねてみたい。傾きを求めるには微分だ。

というわけでさっきの多項式の両辺を微分してみよう。もちろんf(x)は微分できるとして。

071221eq4

a0は定数なので消えたし、a1についてたxも消えた。 これにもx=0を代入するとa1以外はみんな消えるので

071221eq5

これがf(x)のx=0付近における傾き。なので

071221eq6

これもやっぱ0の付近ではf(x)の近似になってるだろう。それも傾いてる分さっきより精度がいいはず。

この要領でf(x)を繰り返し微分し、x=0を代入していったら係数aはみんな求められる気がする。ためしにもう一回微分し、 そこへx=0を代入してみると、

071221eq7

係数a2について整理すれば

071221eq8

同様にa3も求めると

temp1

a4

071221eq9

と、いうことはan

071221eq10

で求められるだろう。

で、こいつらを最初の多項式に戻してやると

071221eq11

狽使って書き直すと

071221eq12

ところで、何回も微分することで係数が求められたのはx=0としたとき項が一つを残して全て0になるからだった。ということは、 xの周囲を少しいじったら、xに0以外の値を代入してもこれまでと同じような作業ができるだろう。具体的には、 最初の多項式のxの部分を(x-a)とする。

071221eq13

で、こうすると、x=aとしたときに(x-a)の部分が全て0になるので、さっきまでと同じような作業が可能となる。

途中すっとばして狽使った式を示すと、

071221eq14

で、こいつがテイラー級数。目的の関数をこの形にするのがテイラー展開。

さっきやったx=0の付近でのテイラー展開は特別にマクローリン展開とか言う。

ところで、テイラー展開をするにあたっては何回も、というか無限回の微分をする。 ということは目的の関数f(x)は無限回微分してもなお何かが残っているようなものでないといけない。そういう関数としては例えばsin、 cos、exなんかがある。

テイラー展開の何が便利かと言えば、こういう三角関数や指数関数を多項式という比較的計算の簡単なものに変換できる点にある。それに、 別に多項式を全て使わなくとも一次か二次程度の項まで使うだけでも多くの場合実用上差し支えのない精度が得られる。つまり

071221eq6

であってもxが0からそれほど離れていないところでのf(x)を計算するには多くの場合有効であるということ。

あと、関数が無限回微分できるということのほかに気をつけないといけないことがある。「x=0の付近では」 とかわざわざ言っているのは、xが0からいくらか離れてしまうとf(x)を表すことができない場合もあるということだ。 足す項を増やせばいいというわけでもない。足す項を無限に増やしていけば完全にf(x)へ収束するようなものもあれば、 いくら項を足していってもある範囲でしか収束しないようなものもある。

また次回くらいに実際に適当な関数をテイラー展開し、Rでプロットでもして確認してみよう。

posted by Rion778 at 01:24 | Comment(0) | TrackBack(0) | 勉強ノート | このブログの読者になる | 更新情報をチェックする
2007年12月15日

PCを使った2階微分方程式の数値解法

前回2階微分方程式をがんばって解いたんだけど, 一般解なんかどうでもいいからとにかく数値解がほしい,そんなときはどうすべきだろうか?

話が1階の微分方程式なら簡単だ.以前やったようにEuler法なんかが使える.

が,2階ではそうは行かない.

1階で上手くできたのは,「xがある値のときのxの増え方」を簡単に求められたからだ.しかし,2階の場合は 「xがある値のときのxの増え方の増え方」を示しているのだ.こんなややこしいものをそのまま使うことは難しい.

そこで,こいつをバラしてやる必要があるんだけれど,その前に前回の微分方程式を持ってこよう.

071215eq1

これが前回扱った単振動の微分方程式だった(項は移動してある).

この式が意味するのは,「zの増え方の増え方」が「z」により定義されるということ.ところで,位置の2階微分, つまり加速度は速度を微分したものだから,速度をvとして次のように書いてもいいだろう.

071215eq2

そして,速度vは位置を微分したものなのだから,

071215eq3

ここで今作った2つの式に注目すると,どちらも一階の微分方程式になっていることがわかる.2元連立一階微分方程式だ.

上の式を使うと,ある時点での速度,位置からちょっと時間が経った後の速度が計算できる.

下の式を使うと,ある時点での速度からちょっと時間が経った後の位置が計算できる.

速度と位置,この式の中で時間に伴い変化する2つの値は,どちらもPCを使って計算できる状態になったわけだ.

それでは例によってRでスクリプトを書こう.

t <- 0
z <- 0
v <- 1
Kdivm <- 1 deltat <- 0.01 i <- 1 Z <- numeric(0) time <- numeric(0) while(t <= 15){ ov <- v Z[i] <- z time[i] <- t v <- v - Kdivm * z * deltat z <- z + ov * deltat t <- t + deltat i <- i + 1 } plot(time,Z,type="l") grid()

t(時間),z(つりあい位置からの変位),v(速度),Kdivm(K/m)の中へ初期値を代入し, deltatへどれだけ短い幅に時刻を設定して計算するのかを入れる.

そうしてスクリプトを実行すると,次のような出力が得られる.

071215fig1

これは当然時間幅を短く取れば正確に,長く取れば不正確になる. 上のグラフでは時刻15までの重りの位置の変化をほぼ正確に(少し振れ幅が大きくなっている)表しているが, 例えば時間間隔を10倍の0.1に設定して計算しなおすと,

 071215fig2

このようにかなり誤差が大きくなってしまう.このあたりはやはりマシンパワーとの相談. Rではこの100倍程度の計算量でもかなり重くなってしまうことがあるので注意が必要.

まあ,あくまで勉強用のスクリプトなので.

 

posted by Rion778 at 22:15 | Comment(0) | TrackBack(0) | 勉強ノート | このブログの読者になる | 更新情報をチェックする
2007年12月14日

2階微分方程式の解法(単振動の運動方程式)

ばねにぶらさがり,上下にゆれる重りの運動を考えよう.つりあい位置からの変位をzとすると,重りに加わる力Fは

F=-Kz

で表すことができる(Kはばね定数).つりあい位置から遠いほど,つりあい位置へ向かう力が大きくなる.まあ当たり前の話だ.

次に,ニュートンの運動方程式を思い出そう.力Fは質量mと加速度aの積として次のように表現できる.

F=ma

加速度というのは速度を微分したものであり,速度というのは位置の変化量を微分したものだ.だから,加速度は位置の2階微分だ! ということが良く分かるように,次のような書き方をしよう.

F=mz''

zの上の2つの点で2回の微分をしていることを表している.

ここでこいつを最初の式とくっつけよう.

mz''=-Kz

両辺をmで割ってからちょいちょいと整理すると

z''+Kdivmz=0

2階微分方程式.こいつを解けば単振動が出てくるだろう.

「解く」というのは上の式におけるzを求めること.zは重りの位置なのだから運動していれば時間により変化する. つまり求めるのは時刻tの関数となる.

で,何を入れれば良いのか.例えばsinなんかは2階微分すると-sinになるので都合が良いかもしれない.次の式を入れてみる.

z=sinOmegat

ここでΩはまだ決めていない定数.こいつはあとで上手いこと調整してやらないといけない.とりあえずそのまま代入してみると

071215eq7

''は2回の微分を表す.sinを微分するとcosになり,cosを微分すると-sinになる.つまり2階微分で-sinになる.で, 定数は微分するたびに前に出てくるので,

071215eq8

sinΩtでくくると, 

071215eq9

これでΩの中身はK/mの平方根じゃないといけないことが分かった.代入する前の式のΩへk/mを代入してやると,

071215eq10

これでzの中身が分かった.変位zは時間と共にこのように変化するのだ.これでめでたしめでたし…ではない.

同じように適当に式を入れて考えていくと,全体を定数倍したもの,sinをcosにしたものも解となることが分かる。

その上,この微分方程式は線形微分方程式と呼ばれるものに属しているので, 定数倍したそれぞれの解を足し合わせたものもやはり解となる.つまり,

071215eq11

も解である(代入すれば確認できる).

これ以上の解は無いので(定数を上手くいじれば他の解はみんな出てくる),これが完全な解,「一般解」である. これでようやく終わった.ちなみに,ここで現れてきた定数a,bは「積分定数」と呼ばれる.

あとは初期条件などから定数を決定してやれば,おもりの運動が分かるというわけだ.

次回は2階微分方程式の数値解を力技で求める方法なんかを.

ところでようやくdvipngの使い方が分かった. コマンドプロンプトでやるときはシングルクオーテーションをダブルクオーテーションにしなければいけない,ただそれだけの事だった. どうして今まで気付かなかったんだろう.とりあえずマトモな数式が書けるようになって良かった.

参考

posted by Rion778 at 20:11 | Comment(2) | TrackBack(0) | 勉強ノート | このブログの読者になる | 更新情報をチェックする
2007年12月10日

PDFをPSへ変換する

powerdotでPDFを取り込もうとしてもどうも上手くいかない。PS形式だったら簡単なのに。

とりあえずGSviewを使ってその場しのぎ。以下、GSview4.8、Ghostscript8.54使用。

  1. GSviewを開く
  2. File→Convertで開くファイルを聞いてくるので、目的のPDFを指定
  3. Deviceにpswriteを選択
  4. Resolutionに600を選択
  5. OK。保存先を聞いてくるので適当に指定。拡張子忘れずに。

これでとりあえずPSになる。

が、スゲーカクカクだ。ビットマップにしてから変換してんのか。

OOoでPS出力できれば…無理か。PDFで出力できるだけマシか。

そういやOOoバージョンアップ忘れてたな。

どうもここのところ珍しく忙しい…気がする。気がするだけだな。やることが沢山あるだけだ。やることが沢山あるのは悪いことじゃない。

posted by Rion778 at 00:37 | Comment(0) | TrackBack(0) | PC関連。HTMLとか,Linuxとか,Rとか | このブログの読者になる | 更新情報をチェックする
2007年12月01日

通風乾湿計のための計算式をRで

乾湿計、と聞いてパッと姿が思い浮かべられるだろうか。

小学校の教室や保健室によくおいてある、普通の温度計と、ガーゼで湿らせて計る温度計の二つがセットになった温度計のことだ。

ガーゼで湿らせたほうの温度計は、湿度を計測するためについている。

空気が乾いているほど水は良く乾く。水は乾くとき、気化熱として熱を奪い去る。だから、どの程度熱が奪われたかを測定すれば、 空気の湿り具合が分かるというわけだ。

それで、実際にはどう使うのか。それには、乾湿計に書かれている表を使う。表には、乾球温度という数字が縦に、 乾湿示差という数字が横に書いてある。まずは乾球(湿らせてない温度計)の温度を読み取り、それに該当する列を表から探す。次に湿球 (湿らせてある温度計)の温度を読み取り、乾湿示差(乾球との温度差)を計算し、さっきの列からその値に該当する部分の数字を探し、 読み取る。その値が湿度(相対湿度)である。

ちなみに、学校とかによくあるタイプの乾湿計は簡易乾湿計と呼ばれるもので、「簡易」の名からも分かるように精度は良くない。 なぜなら、ガーゼから水が蒸発することにより周辺の空気の湿度が高くなってしまうためである。 気流があればその空気はすぐに新しいものと交換されるので問題は無いが、室内程度の気流ではそれは期待できない。

そこで、湿球の周りに気流を人為的に作ってやり、 より正確な湿度を計測しようというのが通風乾湿計である。例えば、アスマン通風乾湿計と呼ばれるものがその一種。

通風乾湿計の精度は非常に高く、例えば乾球と湿球の表示温度差に±0.2℃の誤差があるとすると、相対湿度の誤差は±1〜2%である。

ところで、相対湿度を10秒おきに調べたいというような場合には、乾球と湿球の温度を10秒おきに記録するといった作業が必要になる。 これはそのための機械さえあればさして難しいことではない。

が、測定された温度をいざ相対湿度へ変換しようというとき、さっきのようにいちいち表を調べていては莫大な手間がかかる。

そこで、乾球温度、湿球温度から相対湿度を計算する計算式が必要となる。当然そういったものは存在する。 簡易乾湿計についている表だって計算式によって作られている。

その前に相対湿度の定義を確認しておこう。相対湿度は次の式で定義される。

hr

Uを相対湿度として、eは空気中の水蒸気圧、esはその空気における飽和水蒸気圧である。 蒸気圧とその空気に含まれる分子数には比例関係がある。飽和水蒸気圧というのはその空気が含める「精一杯」の水分子数なので、 つまり相対湿度とは、いま空気に含まれている水蒸気分子の数が、その空気の「精一杯」に対して何割であるか、を示した指標なのである。

ここで求めるべきものは2つある。一つは空気中の水蒸気圧e、もう一つは飽和水蒸気圧es

まずは空気中の水蒸気圧だが、これは次のSprungの式がよく使われる(ところで、 手元の資料では文献の発表年が1855年になっている。一方でWikipediaなどには "Adolf Sprung 1848-1909" という記述が。一体いくつのときにこの式を作ったのだろう)。

sprung

ここでeswは湿球温度における飽和水蒸気圧、tは乾球温度、twは湿球温度、Pは大気圧である。Aは定数で、 湿球が氷結していないときは0.000662、氷結しているときは0.000583の値を入れる。

ここで問題なのは湿球温度における飽和水蒸気圧というやつである。相対湿度を計算するには、 乾球温度における飽和水蒸気圧も求めないといけない。幸いなことに、温度から飽和水蒸気圧を計算する式というのも存在する。

これにはいくつか種類が存在し、どれも近似式だが常温付近なら精度に問題は無い。以下はJISで採用されているSonntagの式。

sonntag

Tには絶対温度を入れる。長い上に返ってくるのは対数だが、PCで使うのなら大した問題ではない。なお、この式が適応できるのは「水」 の飽和水蒸気圧である(つまり0〜100℃が適用範囲)。

さて、これで水蒸気圧、飽和水蒸気圧を求める道具がそろった。あとは計算するだけ。例によってR用の関数を定義しよう。

#温度t(セ氏温度)における飽和水蒸気圧(Pa)
saturation <- function(t){ 
 T <- t+273.15
 es <- exp((-6096.9385/T) + 21.2409642 - 0.02711193*T +0.00001673952*T*T +2.433502*log(T))
 return(es)
 }
#乾球温度t、湿球温度twにおける相対湿度(%)
tuhu <- function(t,tw){
 esw <- saturation(tw)
 A <- 0.000662
 P <- 101325
 e <- esw - A*P*(t-tw)
 es <- saturation(t)
 hr <- e/es
 return(hr*100)
 }

使い方は「tuhu(乾球温度, 湿球温度)」。共にセ氏温度で与える。

ためしになんか適当に数字を入れてみよう。たとえば乾球温度20℃、湿球温度15℃だったとすると、

> tuhu(20,15)
[1] 58.5798

相対湿度はおよそ59%。

ちなみにベクトルで与えることも可能なので、何千、何万のデータでもok。エクセルで扱えないほど多くてもきっと大丈夫。

参考

posted by Rion778 at 19:29 | Comment(2) | TrackBack(0) | 勉強ノート | このブログの読者になる | 更新情報をチェックする

広告


この広告は60日以上更新がないブログに表示がされております。

以下のいずれかの方法で非表示にすることが可能です。

・記事の投稿、編集をおこなう
・マイブログの【設定】 > 【広告設定】 より、「60日間更新が無い場合」 の 「広告を表示しない」にチェックを入れて保存する。


×

この広告は1年以上新しい記事の投稿がないブログに表示されております。