individual_entry.blog.li.nu
日々外観・機能性が進化し続ける、内容不定の日記帳
published:
10
2010
12

偏微分方程式 [3] - 有限区間の 1 次元拡散方程式

Menu : back to the index / comment (0) / trackback
Info : posted at 09:43 PM, categorized as メモ帳 img 数学 img 偏微分方程式

 今回からは 2 階線形偏微分方程式についてです。最初は、2 階の中でも最も簡単な、有限区間の拡散方程式の一般解および特解の求め方についてまとめておきます。

1 次元拡散方程式

001.png
で示される 2 階線形同次の偏微分方程式を、1 次元の拡散方程式と呼びます。a は (一般に正の) 定数です。どうやってこの式が出てきたかは後に譲るとして、これをどうやって解いていけばいいのか考えます。
 この方程式は、熱の広がりがモデルとなっています。したがって、最終的にはそのモデルの細かな条件設定に合うような特解を一般解から求めることが要求されます。例えば、長さの決まった棒に分布している熱の広がりを知りたいかもしれませんし、あるいはほとんど無限といってもいいくらい広い範囲にある物質の拡散状況を知りたいかもしれません、また、熱が最初からどのように分布しているのか・・・など、実はこういう条件設定によってずいぶんと解法が変わってきます。要するに常微分方程式でいう初期条件的なものがあると特解の形も変わってくるというわけです。

初期条件と境界条件

 常微分方程式とは違い、偏微分方程式の条件設定には 2 つあります。まあ、変数が 2 個あるんだから x と t それぞれについての初期値みたいなのがあるんだろうなというイメージで大体分かりますね。それらは具体的には初期条件 (Initial Condition、略して I.C. と書かれる)、境界条件 (Boundary Condition、略して B.C. と書かれる) と呼ばれています。
 初期条件とは、偏微分方程式の未知関数 u が時間 t と 位置 x の関数 u(x,t) であるとすれば、常微分方程式同様 t=0、つまり「観測が始まった時点で未知関数はどんな風でしたか」、という条件です。拡散方程式でいえば、観測が始まったとき、熱は問題にしている棒状でどのような分布になっていましたか、ということですね。ここで、分布という言葉を使っているように、気をつけなければならないのは、常微分方程式は変数が一つしかないので初期条件は単なる一つの値、定数です。しかしながら、偏微分方程式の初期条件は、単に t=0 としても u(x,0)、すなわち、まだ x の関数ではあります。ここが常微分方程式との大きな違いになります。したがって、偏微分方程式の初期条件は、一般に x の関数で与えられることとなります。
 t の初期値に対する u の値が初期条件なら、もう一方の変数である x の値に対応する u が境界条件だということは容易に想像がつきますが、境界条件とは、「問題にしている領域の端点での熱は時間によってどう変化していますか」、という条件です。もし、無限に長い棒の分布を知りたいなら、境界条件は存在しないことになります (もっとも、この場合は「少なくとも発散はしない」というのが境界条件として与えられるようですが)。以下に、初期条件と境界条件のイメージを載せておきます。
002.png
図 : 初期条件 (I.C.)
003.png
図 : 境界条件 (B.C.) は 1 つだけじゃない時もある

有限区間拡散方程式の解法 Step.0 - 与えられた問題の確認

 実際に拡散方程式を解きます。今回は、以下のような条件が与えられていたとします。なお、定数 a は 1 としておきます。一番簡単な場合ですね。
004.png
この問題では、長さ 1 [m] の棒の熱の状態を知りたいとします。初期条件より、観測が始まった時点では、中心である x = 1/2 から左右 1/4 ずつに一様な熱の分布があり、残りには分布がない状態だというのが分かります。また、境界条件より、端点では常に熱はないままであるという条件も課されています。

有限区間拡散方程式の解法 Step.1 - 変数分離法を適用し、偏微分方程式を常微分方程式の組に変換する

 (同次) 偏微分方程式を解く上で非常に強力な道具となるのが変数分離法という方法です。常微分方程式でこの名前見たぞって思うかもしれませんが、それとは全く別のものです。偏微分方程式における変数分離法とは、解である u(x,t) を、それぞれ x だけの関数 F(x), t だけの関数 G(t) による積
005.png
で表されると仮定し、そのまま元の方程式に代入することで F(x), G(t) を確定させるという方法です。後で実際にやれば分かりますが、これは偏微分方程式を常微分方程式に変換してくれる神のツールです。いきなりこんなことをやってもいいのかって感じですが、ちゃんとこうやって仮定すれば必ず解が出てくるので、とりあえず使わせてもらいましょう (ちなみにこれは同次方程式を解くときの常套手段なので他の種類の方程式でも結構な頻度で見ることになります。この方法は拡散方程式だけのためのものではありません)。実際、これを用いると、元の方程式にあらわれている偏微分は
006.png
と、常微分に変換されたのが分かります。t の微分はプライム (') のかわりにドットを上につけるので、上の画像のように書くこともできます。分数形式の方がいいかどうかはお好みですね。実際元の方程式が常微分に分離される過程は、これから元の方程式にこれらを代入することで見えてきます。以上を代入すれば、
007.png
となりますが、常微分にくっついている係数に注目します。t の常微分には x だけの関数 F(x) が、x の常微分の方には t だけの関数 G(t) があります。ですから、両辺とも F(x)G(t) で割ってやれば、両辺とも x だけの式、t だけの式にしてやることが可能です。つまり、両辺それぞれを常微分方程式に出来るかもしれないってことですね。実際、両辺を割れば
008.png
ですね。でもこれでは両辺がそれぞれ x,t だけの式になったというだけで、別に常微分方程式になったわけではありません。さらにこれらが何か定数とイコールであってほしい訳です。都合がいいことに、これはさらに任意の定数 k と等値と見ることができます。すなわち、
009.png
です。なんでこんなことしていいの?って感じですが、よく考えると、左辺は t だけの式、右辺は x だけの式であり、それらが常に等しいのです。つまり、x, t の値に関係なく、恒等的に等しくなければならないのですね。x, t の値に無関係に等しいということは、当然それぞれの式は x, t に関係ない・・・そう、定数でなければなりませんよね。というわけで適当な任意定数 k と等しい、と置くことができるのです。これより、式が 2 つに分離され、偏微分方程式を解くことは 2 つの常微分方程式
010.png
を解くことと同じになります。これなら解けますね。もうちょっと見やすくすると、
011.png
k は定数ですからどっちも定数係数の同次線形常微分方程式、これを解くのは簡単ですね。

有限区間拡散方程式の解法 Step.2 - x に関する 2 階常微分方程式の一般解を求める

 x に関する方程式を解いてみます。2 階線形常微分方程式の解法は既に書いたのでいいと思います。手順通りにやっていくと、特性方程式は
012.png
ですが、k が一体いくつか分からないので、場合分けが必要です。 λ の値がルート含みだと表記上面倒なので、ここでは常に正の実数 p を用いて一般解を表記することにします。つまり、正の場合は
013.png
とし、λ=プラスマイナス p という風にして一般解を書きます。要するに λ をルートなしにしたいので、k のルート自身を新たに p とおきましたってだけのことですね。この場合は、
014.png
という解を得ます。
 次に、k = 0 の場合を考えます。このときは
015.png
ですね。
 最後に、k が負の時ですが、このときもルートを避けるためだけという理由で
016.png
とおいておけば、
017.png
を得ますね。
 以上より、
014.png
015.png
017.png
という、3 種類の一般解が考えられますが、この中で問題に適するのは 1 つしかありません。じゃあそれはどれなんだということで、この時点で問題に適する解はどれなのか吟味して、合わないものは切り捨てていきます。最終的に答えじゃないものまでいつまでも本当の答えと同格に扱っても面倒なだけですからね。
 これは、x の値に関する条件である境界条件から分かります。端点 x=0, 1 ではそれぞれ、変数分離法によって u = F(x)G(t) であるので、u(0,t) = u(1,t) = F(0)G(t) = F(1)G(t) = 0 でなければなりません。すなわち、あらゆる G(t) に対して恒等的にこれらが 0 であるためには、F(0) = F(1) = 0 でなければならないということです。
 次は、このことをそれぞれの場合にあてはめて、どれが問題に合うんだろうとやっていくことになります。まず、一番上の式にあてはめると
018.png
となります。しかし、これは C1 = C2 となるしかありません。Cramer の公式を頭に思い浮かべればすぐだと思いますが、右辺ベクトルが 0, 0 の時点でいずれの C も 0 になってしまいます。もちろん、下の式の両辺に e をかけるか e で割るかして上の式に代入して・・・と考えてもよいでしょう。それで、もし両方とも 0 で都合が悪いかということですが、この場合は F(x) = 0、つまり u は F(x)G(t) ですから、解自体が恒等的に 0 です。初期条件を見るに、全然そんなことはないですね。したがって、これは不適となります。
 真ん中の式はもっと分かりやすいです。同じように代入すれば
019.png
ですからね。したがってこれもダメです。
 では、最後の希望はどうかというと、
020.png
となります。C1 はもうどうしようもなく、0 となるしかありません。そして、一見これも C2 ともに 0 になるしかないように見えます。でも、変数分離法なら解けるっていったはずなのに全部ダメ?もう何も信じられない・・・と早まらずに、最後の式をよく見てみましょう。
021.png
のところです、これって別に C2 が 0 にならなくても 0 になると思いませんか?要するに sin の中身次第では C2 が 0 でなくても 0 になれるでしょう、ということです。だって、sin は π の整数倍なら常に 0 ですからね。というわけで、
022.png
とします (上の画像では 0 からとなっていますが、-1, -2 ・・・・も含めておいてください)。これなら全然 OK ですね。したがって、C2 が 0 とならないで、この式が 0 になるような p は全て境界条件、初期条件からいって問題ないといえますので、この m それぞれに対応する、そのときそのときの C を新たに Cm と置き直せば、結局、この x に関する微分方程式の一般解は
023.png
と得られました。長々と書きましたが慣れるともうほとんど一瞬で判断がつくはずなのでここまで丁寧に書かなくていいと思います。

有限区間拡散方程式の解法 Step.3 - t に関する 1 階常微分方程式の一般解を求める

 t の方は簡単です。1 階線形なら解の公式がありましたよね。解の公式については既に導出をしたことがあるため詳しくは省略し、直接あてはめますと、Q(t)=0 なので頭の中で暗算出来て
024.png
となります。丁寧に書きましたがこれはすぐに出せると思います。常微分方程式をそれなりに学習したなら覚えているはずですね。で、F(x) の決定時に
016.png
とおきましたから、これは要するに マイナス p の 2 乗、そしてさらに p は mπ だったので
025.png
を得ます。これで G(t) も決定することができました。

有限区間拡散方程式の解法 Step.4 - 解を重ね合わせ、偏微分方程式の一般解とする

 以上より、この偏微分方程式の解は
026.png
であることが分かりました。なお、G(t) の C と F(x) の Cm との積を改めて Cm とおきなおしています。こういうやり方が嫌であれば、Dm 等おいておくのが無難でしょう。
 偏微分方程式において、同次方程式の解は、重ね合わせが有効です。すなわち、上の整数 m がどんな値であろうが方程式の解となるため、これらを全部重ね合わせて一般解とする (汎用性を上げて、より多くの特解に対応できるようにする) 必要があります。となると、マイナス無限大から無限大までの総和での足しあわせになりますね。
027.png
m が整数っていうのも、良い条件です。そのまま総和記号として表現できますからね。もしこれが連続値だったら積分になります (無限区間の拡散方程式だと積分ですがこれは後の回でまとめます)。この時点でなんとなく想像がついたかもしれませんが、特解は、初期条件を見るに、どう見ても不連続です。こんな変な関数を表現しようとするなら、もうフーリエ級数展開するくらいしかないでしょう。そこで、上の式を見てください。無限級数 (マイナス側はあるけど) だし、初期条件に合うよう t=0 とすると e が消滅して、「ΣCmsin(mΠx) = t(0,x) = f(x)」 となりますよね。ということは、どうもこれフーリエ級数展開すればうまく特解を出せるっぽいです。というわけで、次はフーリエ級数展開できるよう、どうやってマイナス側を消して m を 1 から無限大の無限級数にできるかな、と考えます。
 ここで、負の場合を -m 番目と表現し、Σ の中で分離させます。すなわち、
028.png
です。さらに m=0 では sin が 0 なので、m=0 は除外して m=1 からとしても問題ありません。これでフーリエ級数展開の定義と同じ、1 から無限大の無限級数となりました (画像は m=0 となっていますが m=1 の間違いです。以下、無限級数が m=0 となっていますが全て m=1 の間違いになっています)。あとはどうにかしてこれを 1 つの項にまとめられれば美味しい、というわけです。e の累乗にくっついている m は、マイナス m としても 2 乗ですから m に戻ります。また、sin は奇関数なのでマイナスが入るとマイナスは外に出せます。ですから、上の式は
029.png
となりますね。ここまで来たら勝利も同然、なぜならば、任意定数同士の差、これはどっちも任意だからまたまとめて任意定数とおけばいいだけの話なのです。したがって、この C の差同士を改めて C とおけば、
030.png
という一般解を得ます。これであとは初期条件、境界条件に合う解を見つけ出すだけとなりました。

有限区間拡散方程式の解法 Step.5 - 初期条件・境界条件を考慮し、特解を求める

 最後に求めなければならないのは、m によって異なる係数 Cm の全てです。をれを決めるヒントは初期条件ですね。t=0 を代入すれば
031.png
であり、これは f(x) のフーリエ級数展開そのものとなります。cos の項はないですけどね。つまりこれは奇関数として展開しなさいということをいっています。でも、問題の初期条件は、図示してみると
032.png
のようであり、周期 1 と見ると、奇関数の要素がありません。仮に真ん中を中心として周期 1/2 にしても偶関数であって奇関数 とはなりえません。これは手詰まりなのでしょうかということですが、ここで、半区間展開というものを行います。要するに、関数を無理矢理奇関数にします。といっても形を崩したら元も子もないので、このパルス状の図形を奇関数的に、x 軸の負の領域に配置します。図的には以下のようです。
033.png
このように、単純に線対称として折り返して配置したものは偶関数に見えるので偶関数的拡張、点対称で折り返したものは奇関数的拡張といいます。今回は奇関数として展開することが求められているため、奇関数的拡張、図で言えば黄色い方を付け足します。これで、周期 2 の奇関数となりましたので、フーリエ級数展開が可能となります。こんなことやってもいいのかと思うかもしれませんが、このままフーリエ級数展開して得られたものは、この周期 2 のパルス状、エセ正弦波が延々と繰り返されているものとなります。で、このうち x=0 から 1 までを定義域にしておけば、何も問題ないのです。要らない領域は定義されてないからどうでもいい、と切り捨てればいいだけってことですね。
 これで、係数 Cm が求まります。これはフーリエ正弦級数ですから、公式によって
034.png
ですぐに求まります (上の式において sin の中身は mπx ではなく正しくは mπx/L でした)。この場合は周期 L = 1 なので、それを代入します。まあ、本当は上の画像における f(x) は奇関数拡張後のものなので g(x) とでもしておくのがいいかもしれませんがここでは省略してあります。で、f(x) は奇関数であり、sin も奇関数ですから、これは偶関数となって、実際の積分区間は半分で良くなります。
035.png
ここまで来て、ようやく計算らしい計算が始まりますね。
036.png
となります。これはこれでいいんですが、三角関数は幸運なことに和積公式とかいうのがありますから、これ、1 つの項にまとめられます。詳しくは既に書いたので省略することにして、和積公式とは
036.png
という公式でした。このうち、真ん中のものが適用出来ますよね。α が 2 ぶんの mπ、β が 4 ぶんの mπ です。したがって
037.png
と求まりました。結局これは m が偶数の時は左の方の sin のせいで無条件に 0 になって、生き残るのは奇数の時だけ、ということになります。m は有限の値ではなく無限大まで行くので、偶数の場合を切り捨てても同じことです。したがって、上を奇数の場合のみに書き換える、すなわち、m を 2m-1 と置き換えれば、
038.png
を得ました。これで係数が全て求まったため、元の一般解に代入し、最終的な解
039.png
を得ます。半端ない長さですね。

解の動きの確認

 解が実際にどうなっているかは、グラフをコンピュータに描かせるしかありません。総和まで計算してグラフを書いてくれるソフトは少ないため、ここではそれが出来る GCalc-Plus というソフトで描きます。まず、t=0 で、本当に初期状態とグラフの形状が一致しているか確認します。
040.png
本当は m = 無限大まで描かせるべきなのでしょうが、ここではそんなことをしたらコンピュータが死んでしまうので 999 番目までの総和で表現しています。なので、不連続点付近はちょっと歪みが見えますが、おおむね正しい形状になっているのが確認できます。これが、時間を進めていくことでどのように拡散していくのか見ていきます。
041.png 042.png
043.png 044.png
045.png 046.png
047.png 048.png
このように、たった 1 秒経過しただけであっという間に熱が失われてしまっているのが分かります。

例題

049.png
 例題にも挑戦します。定義域は 0050.png とおきます。すると、x,t の偏微分は常微分に変わり
051.png
となります。これを代入すれば、方程式は
052.png
ですが、両辺 FG で割ると、それぞれ x, t だけの式となり
053.png
となります。a はどっちにつけておくかですが、1 階より 2 階の方が圧倒的に解くのが面倒なので 1 階常微分方程式になる t の方に a をくっつけておきます。ここで、右辺左辺それぞれ x,t に関わらず恒等的に等しくなければならず、そのようになるにはこれらは何らかの定数 k でなければなりません。したがって、
054.png
を得ます。これより、偏微分方程式は二つの常微分方程式に分離され、
055.png
を解けばいいということになります。
 x の方程式の一般解は、k の符号ごとに場合分けされます。
056.png
境界条件を代入すると、u(x,t)=F(x)G(t) より u(0,t)=u(2π,t)=F(0)G(t)=F(2π)G(t)=0, あらゆる t に対して成立するためには F(0)=F(2π)=0 でなければならず、上 2 つはこれを満たすのは C1=C2=0 のときのみとなります。それは u(x,t) が恒等的に 0 であることを意味しますが、初期条件より明らかに解ではありません。したがって除外します。一番下に同じく代入すれば、C1=0 と、
057.png
という条件を得ます。C2 は 0 にならなくても、sin が 0 になるような条件ならば初期条件を満たすことになります。それは、sin の中身 2πp が π の整数倍になるということであり、
058.png
となります。これより、それぞれの m に対応する任意定数 C2 を新たに Cm とおけば
059.png
と定まります。
 G(t) の方は、解の公式にあてはめることで、直接
060.png
を得ます。ルートを避けるため、k はプラスマイナス ip の二乗としてあった点に注意です。
 以上から、u(x,t) は
061.png
となります。C の m と C プライム m の積は改めて Cm とおきなおしてあります。ここで、全ての整数 m に対してこの u(x,t) は解ですから、重ね合わせても解になります。したがって、
062.png
も解です。そして、この解はさらにまとめることができて、m が負の場合を m → -m としてシグマの中で分離させます、そうすれば
063.png
また、m=0 の場合は sin があるので 0 です。したがって、総和は m=1 からとすることができ、C_m と C_{-m} の差を改めて C_m とおきなおせば、一般解
064.png
を得ます。
 次に、初期条件へのあてはめを行います。t=0 のとき、
065.png
です。Cm を求めたいのですが、これは f(x) のフーリエ級数展開 (L=2π, 周期 2L=4π) の式そのものであり、Cm をフーリエ正弦級数として求めることができます。このフーリエ級数展開は sin の項しか出てないので奇関数でなければなりません。したがって f(x) を奇関数拡張して、
066.png
となります。あとはこの計算をすればよいだけです。面倒ではありますが難しくはないですね。難しいと面倒は違いますからね。
067.png
まあ長いだけでたいしたことないと思います。これは偶数の時 0 となり、奇数の時は 1-(-1) の m 乗は 2 になるので、奇数の時だけ生き残るように書き換えます。具体的には m → 2m-1 とすると
068.png
が得られました。したがって、特解は、上より奇数番号だけ生き残らせなければならないこと (つまり m → 2m-1 になってること) に注意すれば
069.png
が答えです。このように、特解の形が最初の例とは全く違う形になっているのが分かります。一般に、偏微分方程式の特解は初期条件と境界条件に大きく依存します。以下、おまけです。
070.png 071.png
072.png 073.png
074.png 075.png
076.png
これも、元の形状は違いますがあっという間に熱が逃げていってしまっているのが分かります。実際のモデルで考えれば、こんなすぐに熱がなくなってしまうことは考えにくいです。これは、定数 a=1 というのが相等熱の伝わりがいい時の値なんだなということになります。

trackback

Trackback URL : http://li.nu/mt/mt-tb.cgi/265

comment

投稿する

投稿状況

このウェブページにはまだコメントがありません。右側の記入欄に必要事項 (メールアドレス、URL は任意です) を記入して、投稿できます。この Blog にユーザー登録すると、この手間を省け、独自の画像を表示できたりします。

powered by

Movable Type 5.0
Movable Type