微分方程式の境界値問題を解く
１。微分方程式を解くって何だっけ？
このワークシートはMath by Codeの一部です。
前回まで、微分方程式の型や解法、それで得られた関数や多項式を見てきた。
初心に帰って、そもそも微分方程式を解くって何だっけ？
から考えてみよう。
＜解の種類＞
微分方程式には一般解と特殊解があったね。
独立変数をｘとして、ｙがｘの関数のとき、ｘの値に関係なく微分方程式を満たす解が一般解だった。
それに対して微分方程式のｘなどのパラメータに条件をつけたものが特殊解だった。
そもそも、微分方程式は独立変数を時間ｔとして、ｔの関数としての運動を記述するイメージから始まった。そのときの条件というのは、ｔ＝０のときの状態ｙがどうかという初期値だった。
その初期条件から未来のｔの状態ｙを予測するというイメージだったね。これを初期値問題という。
微分方程式ではもう一つ大切な使い方がある。それが、境界値問題だ。
独立変数ｘが空間で位置を表すときに、
ｘの変域の境界での状態ｙはいくつか？どう変化するか？という境界条件を与えるときに、
内部の状態ｙがどんな分布に落ち着くのか、どう安定するのか？
などの問いが生まれてきた。
言い換えると、時間の条件で解くのが初期値問題で、空間の条件で解くのが境界値問題だね。
＜境界値問題を解く＞
xの区間I=[a,b]で微分方程式を考える。
区間の両端x=a,x=bでのyまたはy'が境界条件を満たす解と求めることを境界値問題という。
（例）
たとえば、微分方程式y''=exp(x) -2の解の境界値問題を考えてみよう。3つの境界条件がある。
境界条件A：y(0)=0, y(1)=0
境界条件B：y'(0)=0, y(1)=0
境界条件C：y’(0)=0,y'(1)=0
まず一般解。変数が最初から分離しているので、両辺の積分を繰り返す。
y'=exp(x)- 2x+c , y = exp(x) - x2 +cx + d(これが一般解)
A：y(0)=1-0+d=0 ,y(1)=e-1+c+d=0 これから、d=-1, c=2-eとなる。解はy=exp(x)-x2+(2-e)x-1
B：y'(0)=1 + c =0 , y(1)=e-1+c+d=0 これから、c=-1, d=2-eとなる。解はy=exp(x)-x2-x+(2-e)
C：c=-1, y'(1)=e-2+c=0 これから、c,d は決まらず。解なし。
＜境界条件の３パターン＞
境界条件BC(Boundary Condition) を３つのパターンに分けることがよくある。
xの区間I=[a,b]で微分方程式を考える。
区間の両端のｙの値を与えるもので、ディリクレ条件（固定端条件）という。
区間の両端のｙ’の値を与えるもので、ノイマン条件（自由端条件）という。
これらの単純形以外は混合型とか第３種条件ということが多い。
上の例ではAがディレクレ型、Bが混合型、Cがノイマン型だね。
（例）
常微分方程式 y''-xy'-y=0
境界条件 y(0)=1, y'(0)=0
ｘ＝０のまわりのべき級数解ｙを求めよう。
y=a_0+a_1x+a_2x^2+ a_3x^3+ ........... + a_kx^k+..... =Σ(a_m xm) [ m for 0 to ∞]
y'=a_1 +2a_2x + 3a_3x^2 + .......... + ka_kx^(k-1) +......1 =Σ(m a_m xm-1) [ m for 1 to ∞]
y''=2 a_2 + 3・2 a_3 x + ........ + k(k-1) a_k x^(k-2) + ...... Σ(m(m-1) a_m xm-2) [ m for 2 to ∞]
定数項和2a_2-a_0=0から、 a_2=1/2a_0
xの係数和 3・2 a_3 -a_1-a_1=0 から、 a_3=1/3a_1
x2の係数和 4・3 a_4 -2a_2-a_2=0 から。a_4=1/4a_2
xkの係数和 (k+2)(k+1) a_(k+2) -ka_k-a_k=0 から。a_(k+2)=1/(k+2)a_3xk+1
x_even=a_2s=
x_odd= a_2s+1=
y0(x)=1+1/2x2+1/4!! x4+ 1/6!! x6+.....,
y1(x)=x+1/3!! x3+1/5!! x5+......
として、y(x)=a_0y0(x)+a_1y1(x)が一般解だ。次に境界値問題を解こう。
特別解は、
y(0)=a_0*1+a_1*0=1, y'(0)= a_0*0+a_1*1=0から、a_0=1, a_1=0から、
y(x)=1+1/2x2+1/4!! x4+ 1/6!! x6+.....=
２．熱伝導方程式の境界値問題
1次元熱伝導方程式（拡散方程式）uは温度で、位置x,時間tの関数としよう。
∂u/∂t＝c2 ∂2u/∂x2 または表式は、 ut = c2 uxx
初期条件u(x,0)=f(x) ：０とLの間でのt=0の温度分布fはxの関数。
境界条件u(0,t)=u(L,t)=0 と,x=0,x=Lで断熱されているとすると、f(0)=f(L)=0となる。
u(x,t)=g(x)h(t)とすると、境界条件からg(0)=g(L)=0
ut=∂(g(x)h(t))/∂t=g(x)h'(t),
c2 uxx=c2 ∂(g(x)h(t))/∂x=c2 h(t) g''(x)
g(x)h'(t)= c2 h(t) g''(x)となるので、
変数分離してg''(x)/g(x)=h'(t)/(c2h(t))=kとおこう。
p=g'(x),Q=g'(x)とするとP=g(x),q=g''(x)で、g''=kgだから、
部分積分によって、[x for 0 to L] ∫(g')2dx=[ gg']-∫gg''=0-0- k∫g2dxが非負なので、kは負。
なので、k=-p2とおく。
・2階線形常微分方程式X：g''(x)+p2 g(x)=0
・1階線形常微分方程式T：h'(t)+p2 c2h(t)=0
最初の熱伝導方程式は、この２つの常微分方程式に分解されたということだね。
・方程式Xの一般解は、特性方程式t2+p2=0は複素数解t=0+pi, q=0 -piを持つので、
g(x)=e0x(A cospx + B sinp x)= A cos px + B sin pxとなるね。
境界条件から、g(0)=Acos0+Bsin0=A=0 , g(L)=0+ B sin pL = 0 pL=nπ p=nπ/L
以上から、B=1として、
gn(x)= sin (nπ/L x) [n for 1 to ∞1]
・方程式Tの一般解は、から、dh/dt=-(qn)2hn qn=cnπ/Lとする。
変数分離して積分 ∫1/hn dh=-(qn)2 ∫ 1 dt log|hn|= -(qn)2 t +c
hn(t)=Cn exp( -(qn)2 t )
だから、un(x,t) = gn(x)hn(t) ＝Cn sin (nπ/L x) exp( -(qn)2 t )
そこで、u(x,t) = Σ Cn sin (nπ/L x) exp( -(qn)2 t )
[n for 1 to ∞](qn=cnπ/L)を一般解としよう。
初期条件からu(x,0)=f(x)は上式のt=0だから、 ここで、フーリエsin級数、f(t)=Σ an sin[ ] ならan=2/T ∫ f(t) sin[ ] dtを参考にすると、
Cn= とすればよいね。
このとき、u(x,t) = 。
（例）
f(x)=if(x≦L/2, x, L-x)とする。
x=0, x=Lでは断熱されていて初期温度はｘ＝L/２で最高になる三角形の分布f(x)=u(x,0)
Cn=
カンタンのためにsin(nπ/L x)=S(x), cos(nπ/L x)=C(x)とおく。
x={0,L/2,L}に対して、
ｎ=(1,2,3,0) mod 4 のとき、s=(1,0,-1,0), c=(0,-1,0,1), cc=(-1,1,-1,1) とおくと、
S(x)={sin0, sin(nπ/2), sin(nπ)} ={0, s, 0} ,
C(x)={cos0, cos(nπ/2), cos(nπ)}={1, c, cc}
また、部分積分によって∫x S(x)dx
=- L/nπ x C(x) + L/nπ∫C(x)dx =- L/nπ x C(x) + (L/nπ)2S(x) より、
Cn=
これより、
Cn= ｎ=(1,2,3,0) mod 4 のとき、s=(1,0,-1,0)
つまり、nが偶数ならCn=0,
n=1,5,9,....,ならCn , n=3,7,....なら、Cn
u(x,t) =
1次元熱伝導のモデル
３．長方形領域のラプラス方程式のディリクレ問題
すごく長いタイトルですが、
・ラプラス方程式とはラプラシアンをしてゼロ。つまり、u=u(x,y)について、uxx+uyy=0
・たてｙ＝ｂ、よこｘ＝aの長方形領域の中での状態uをみるということ。
・ディリクレ問題というのは、境界条件が固定端ということだった。
u(0,y)=u(a,y)=0 長方形の左右は断熱。
u(x,b)=0,u(x,0)=f(x) 長方形の上端は断熱。下端はfで分布。
・熱伝導方程式と同じように1変数関数の積で表し、変数分離しよう。
u(x,y)=g(x)h(y)。uxx=g''(x)h(y) , uyy=g(x)h''(y)より、uxx+uyy=g''(x)h(y) +g(x)h''(y)=0
g''(x)/g(x) =-h''(y)/h(y) =k とおくと、g''(x)=kg(x), h''(y)=-kh(y)。
・2階線形常微分方程式X：g''(x)- k g(x)=0
・2階線形常微分方程式X：h''(y)+k h(y)=0
・方程式Xの一般解は、熱伝導のときと同様にkが負だから、g(x)= A cos√-k x + B sin√-k xとなるね。
境界条件から、g(0)=Acos0+Bsin0=A=0 , g(a)=0+ B sin√-k a = 0 √-k a=nπ √-k=nπ/a
k=-(nπ/a)2 以上から、gn(x)= sin (nπ/a x) [n for 1 to ∞]とおける。
・方程式Yの一般解は、熱伝導のときと反対に-kが正だから、特性方程式は２つの実数解±√-kをもつ。
だから、hn(y)= An exp(p y) +Bn exp(-p y) [n for 1 to ∞]とおける。簡単のためにnπ/a＝ｐ。
ここで、双曲線関数cosh(x)=[e^x+e^(-x)]/2, sinh(x)=[e^x-e^(-x)]/2を連想すると、
e^xはcoshとsinhの和であり、e^(-x)はcohとsinhの差で表せることに気づくね。式の変形しよう。
hn(y)=An(cosh(py)+sinh(p y))+Bn(cosh(p y)-sinh(py))
hn(y)=(An+Bn)cosh(py)+(An-Bn)cosh(p y)
=Cn cosh(nπ/a y)+Dnsinh(nπ/a y)
境界条件から、u(x,b)=hn(b)=Cn cosh(nπb/a )+Dnsinh(nπb/a )=0 となり、
Cn=-Dn sinh(nπb/a)/cosh(nπb/a)。
hn(y)= Cn cosh(nπ/a y)+Dnsinh(nπ/a y)
=-Dn[cosh(nπ/a y)sinh(nπb/a)ーsinh(nπ/a y)cosh(nπb/a) ]cosh(nπb/a)
=En sinh(nπ/a bーnπy/a） 加法定理sh(x±y)=sh(x)ch(y)±ch(x)sh(y)から
=En sinh(nπ (bーy)/a）
以上から、重ね合わせをして、
u(x,y)=Σgn(x)hn(y)[n for 1 to ∞]
= ΣEn sin (nπ x/a) sinh(nπ (bーy)/a）[n for 1 to ∞]
sinフーリエ級数an=2/T ∫ f(t) sin[ 2π n t/T ] dtとすると、
f(t)=Σ an sin[ 2π n t/T] ) =2/T Σ sin[ 2π n t/T] ) ∫ f(t) sin[ 2π n t/T ] dt だから、
f(0)=f(a)=0から、[-a,a]でf(x)が奇関数で2a周期で無限につなぐために、
f(x)=2/aΣ sin(nπx/a) ∫f(k)sin(nπk/a) dk [k for 0,a] としよう。
境界条件から、u(x,0)=ΣEn sin (nπ x/a) sinh(nπ (b)/a）=Σ sin(nπx/a) 2/a ∫f(k)sin[nπk/a] dkとなり、 En=2/[a sinh(nπ b/a）]∫f(x)sin[nπx/a] dxだね。
u(x,y)=Σgn(x)hn(y)[n for 1 to ∞]
= ΣEn sin (nπ x/a) sinh(nπ (bーy)/a）[n for 1 to ∞]
=
4．波動方程式
・波動方程式は熱伝導方程式の左辺を時間で2階偏微分にしたものだ。
弦の振動を表したりするのに向いているね。
つまり、u=u(x,y)について、∂2u/∂t2＝c2 ∂2u/∂x2 または表式は、 utt = c2 uxx
初期条件 初期変位u(x,0)=f(x) 、初期速度ut(x,0)=v(x)。
境界条件u(0,t)=u(L,t)=0 から、f(0)=f(L)=0となる。
u(x,t)=g(x)h(t)とすると、境界条件からg(0)=g(L)=0
utt=∂2(g(x)h(t))/∂t2=g(x)h''(t),
c2 uxx=c2 ∂(g(x)h(t))/∂x=c2 h(t) g''(x)
g(x)h''(t)= c2 h(t) g''(x)となるので、
変数分離してg''(x)/g(x)=h''(t)/(c2h(t))=kとおこう。
p=g'(x),Q=g'(x)とするとP=g(x),q=g''(x)で、g''=kgだから、
部分積分によって、[x for 0 to L] ∫(g')2dx=[ gg']-∫gg''=0-0- k∫g2dxが非負なので、kは負。
なので、k=-p2とおく。
・2階線形常微分方程式X：g''(x)+p2 g(x)=0
・2階線形常微分方程式T：h''(t)+p2 c2h(t)=0
最初の熱伝導方程式は、この２つの常微分方程式に分解されたということだね。
・方程式Xの一般解は、特性方程式t2+p2=0は複素数解t=0+pi, q=0 -piを持つので、
g(x)=e0x(A cospx + B sinpx)= A cos px + B sin pxとなるね。
境界条件から、g(0)=Acos0+Bsin0=A=0 , g(L)=0+ B sin pL = 0 pL=nπ p=nπ/L
以上から、B=1として、
gn(x)= sin (nπ/L x) [n for 1 to ∞1]
・方程式Tの一般解は、特性方程式t2+(pc)2=0は複素数解t=0+pci, q=0 -pciを持つので、
h(t)=e0x(C cospct + D sinpct)= C cos pct + D sin pctとなるね。pc=cnπ/L=qnとして、
hn(t)=Cncos qn t + Dn sin qn tとおく。
・以上から、一般解は、
un(x,t)=Σgn(x)hn(t)[n for 1 to ∞]= Σ (Cn cos qn t + Dn sin qn t ) sin (nπ/L x)[n for 1 to ∞]
qn=cnπ/L を固有値、それに対するunを固有関数という。{q1, q2,.....}をスペクトル。
qn/2π= cn/2Lが振動数。
初期条件u(x,0)=f(x)から、u(x,0)=Σ qn(Cn cos0 + Dn sin 0 ) sin (nπ/L x)[n for 1 to ∞]
=Σ Cn sin (nπ/L x) [n for 1 to ∞]=f(x)
sinフーリエ級数an=2/T ∫ f(t) sin[ ] dtとすると、f(t)=Σ an sin[ ] だから、
Cn=2/L ∫0L f(x) sin (nπ/L x) dx [n for 1 to ∞]
初期条件ut(x,0)=v(x)から、ut(x,0)=Σ qn(-Cn sin0 + Dn cos 0 ) sin (nπ/L x)[n for 1 to ∞]
=Σ Dn qn sin (nπ/L x)[n for 1 to ∞] = v(x)
sinフーリエ級数an=2/T ∫ f(t) sin[ ] dtとすると、f(t)=Σ an sin[ ] だから、
Dn qn=2/L ∫0L v(x) sin (nπ/L x) dx [n for 1 to ∞]（cnπ/L=qn）
Dn =2/cnπ ∫0L v(x) sin (nπ/L x) dx [n for 1 to ∞]
つまり、
[n for 1 to ∞]
u(x,t) = Σ (Cn cos qn t + Dn sin qn t ) sin (nπx /L)（ncπ/L=qn）
Cn = 2/L ∫0L f(x) sin (nπ/L x) dx
Dn = 2/cnπ ∫0L v(x) sin (nπ/L x) dx
（例）
初期変位f(x)=if(x≦L/2, x, L-x)とすると、初期速度v(x)=0となり、Dn=0。
u(x,t)=Σ Cn cos(ncπ/L t) sin (nπ/L x)（ncπ/L=qn）
x=0, x=Lは固定されていて変位はｘ＝L/２で最高になる三角形の分布f(x)=u(x,0)
Cn=
このあとCnの求め方は、上の２．熱伝導方程式の境界値問題の例と同じなので、
途中省略
n=1,5,9,....,ならCn , n=3,7,....なら、Cn
u(x,t) =