連立微分方程式
実際に扱う微分方程式の問題の多くは、運動方程式など2階の微分方程式や連立微分方程式になります。そこで今回は、高階の微分方程式の解法の演習を行います。高階の微分方程式は、連立微分方程式に帰着できますので、今回の演習は、連立微分方程式の演習になります。
次のような連立微分方程式を考えます。
$\displaystyle{\frac{dy_{1}(x)}{dx}}=f_1(x,y_1,\cdots,y_n)$
$\displaystyle{\frac{dy_{2}(x)}{dx}}=f_2(x,y_1,\cdots,y_n)$
$\cdots\cdots \cdots \cdots \cdots \cdots $
$\displaystyle{\frac{dy_{n}(x)}{dx}}=f_n(x,y_1,\cdots,y_n)$
$y_1(x^{(0)})=y_1^{(0)}, y_2(x^{(0)})=y_2^{(0)}, \cdots, y_n(x^{(0)})=y_n^{(0)}$
これをベクトルを用いて表現すると、次のようになります。
$\displaystyle{\frac{d{\bf y}(x)}{dx}}={\bf f}(x,{\bf y}(x))$
${\bf y}(x^{(0)})= {\bf y}^{(0)}$
ここで
${\bf y}(x)=\left(\begin{array}{c} y_1(x) \\ y_2(x) \\ \vdots \\ y_n(x) \end{array}\right),\ \ \ {\bf f}(x, {\bf y}(x))=\left(\begin{array}{c} f_1(x, {\bf y}(x)) \\ f_2(x, {\bf y}(x)) \\ \vdots \\ f_n(x, {\bf y}(x)) \\ \end{array}\right),\ \ \ {\bf y}^{(0)}=\left(\begin{array}{c} y_1^{(0)} \\ y_2^{(0)} \\ \vdots \\ y_n^{(0)} \\ \end{array}\right)$
です。
この連立常微分方程式を解く方法は、基本的には連立ではない常微分方程式の解法と同じで、4次のルンゲ・クッタ法の公式は、ベクトルを用いて以下のように表わされます。
${\bf y}^{(j)}={\bf y}^{(j-1)}+\displaystyle{\frac{1}{6}}({\bf k}_1+2{\bf k}_2+2{\bf k}_3+{\bf k}_4)$
${\bf k}_1=h{\bf f}(x^{(j-1)}, {\bf y}^{(j-1)})$
${\bf k}_2=h{\bf f}(x^{(j-1)} +\displaystyle{\frac{h}{2}},\ {\bf y}^{(j-1)}+\displaystyle{\frac{{\bf k}_1}{2}})$
${\bf k}_3=h{\bf f}(x^{(j-1)} +\displaystyle{\frac{h}{2}},\ {\bf y}^{(j-1)} +\displaystyle{\frac{{\bf k}_2}{2}})$
${\bf k}_2=h{\bf f}(x^{(j-1)} +h,\ {\bf y}^{(j-1)} +{\bf k}_3)$
ここで、$j$は$x$軸方向のグリッドの番号(ルンゲ・クッタ法による数値的解法の計算ステップ数)を表しています。
それでは、次の連立常微分方程式を4次のルンゲ・クッタ法で解いてみましょう。
$\displaystyle{\frac{dy_{1}(x)}{dx}}=y_2$
$\displaystyle{\frac{dy_{2}(x)}{dx}}=-1$
ただし
$y_1(0)=0,\ \ \ y_2(0)=1$
とします。この連立微分方程式の厳密解は
$y_1=x-\displaystyle{\frac{1}{2}}x^2$
$y_2=1-x$
ですので、数値解と厳密解を比較して、正しく計算できているか確認しましょう。
1 program runge4_sim
2 implicit none
3
4 integer,parameter:: n = 2
5 integer jmax,j
6 double precision x,y(n),h,f
7 external f
8
9 x = 0.0
10 y(1) = 0.0
11 y(2) = 1.0
12
13 h = 0.1
14 jmax = 2.0/h
15
16 write(*,'("x=",1pe12.4," y1=",2e12.4," y2=",2e12.4)') &
17 x,y(1),x-0.5d0*x*x,y(2),1.0d0-x
18
19 do j=1,jmax
20 call runge4( x, y, n, h, f )
21 x = x + h
22 write(*,'("x=",1pe12.4," y1=",2e12.4," y2=",2e12.4)') &
23 x,y(1),x-0.5d0*x*x,y(2),1.0d0-x
24 end do
25
26 end program runge4_sim
27
28 function f( i, x, y, n )
29 implicit none
30 integer i,n
31 double precision f,x,y(n)
32
33 if( i .eq. 1 )then
34 f = y(2)
35 else if( i .eq. 2 )then
36 f = -1.0d0
37 end if
38
39 end function
40
41 subroutine runge4( x, y, n, h, func )
42 implicit none
43 integer i,n
44 double precision x,y(n),h,func
45 double precision,allocatable,dimension(:)::k1,k2,k3,k4,work
46
47 allocate(k1(n),k2(n),k3(n),k4(n),work(n))
48
49 do i = 1, n
50 k1(i) = h * func( i, x, y, n )
51 end do
52 do i = 1, n
53 work(i) = y(i) + k1(i) / 2.0d0
54 end do
55
56 do i = 1, n
57 k2(i) = h * func( i, x+h/2.0d0, work, n )
58 end do
59 do i = 1, n
60 work(i) = y(i) + k2(i) / 2.0d0
61 end do
62
63 do i = 1, n
64 k3(i) = h * func( i, x+h/2.0d0, work, n )
65 end do
66 do i = 1, n
67 work(i) = y(i) + k3(i)
68 end do
69
70 do i = 1, n
71 k4(i) = h * func( i, x+h, work, n )
72 end do
73
74 do i = 1, n
75 y(i) = y(i) + (k1(i) + 2.0d0*k2(i) + 2.0d0*k3(i) + k4(i)) / 6.0d0
76 end do
77
78 deallocate(k1,k2,k3,k4,work)
79
80 end subroutine runge4
allocationの使い方については、Fortran90の基礎の "6.3. 配列の動的割付け"を参照してください。
次に、高階の微分方程式の例として、一様重力場における物体の運動方程式(2階の微分方程式)を考えてみましょう。今、重力加速度の大きさを$g$とすると、垂直方向($z$方向)の運動方程式は、
$\displaystyle{\frac{d^2 z}{dt^2}}=-g$
となります。ここで、
$\displaystyle{\frac{dz}{dt}}=v$
とおくと、この運動方程式は次の1階の連立微分方程式と等価になります。
$\displaystyle{\frac{dz}{dt}}=v$
$\displaystyle{\frac{dv}{dt}}=-g$
このように、2階の微分方程式である運動方程式を解くことは、連立微分方程式を解くことに等しくなるわけです。ここで上式において、
$z\rightarrow y_1,\ \ \ v\rightarrow y_2,\ \ \ t\rightarrow x,\ \ \ g\rightarrow 1$
とすると、プログラムで解いた問題と同じになりますね。つまり先のプログラムは、加速度-1、初速度1の等加速度運動を解く問題だったと考えることができます。
〇課題○ 次の連立常微分方程式を解きなさい。
$\displaystyle{\frac{dy_{1}(x)}{dx}}=y_2$
$\displaystyle{\frac{dy_{2}(x)}{dx}}=-y_1$
ただし
$y_1(0)=0,\ \ \ y_2(0)=1$
とする。また、この連立微分方程式の厳密解を求め、数値計算の結果と比較しなさい。
〇課題○ 次の微分方程式を解きなさい。
$\displaystyle{\frac{d^2 y}{dx^2}}=x\displaystyle{\frac{dy}{dx}}+y$
ただし
$y(0)=0,\ \ \ y’(0)=1$
とする。
〇課題○ 次の運動方程式を解きなさい。
$\displaystyle{\frac{d^2 x}{dt^2}}=-\displaystyle{\frac{x}{(x^2+y^2)^{3/2}}}$
$\displaystyle{\frac{d^2 y}{dt^2}}=-\displaystyle{\frac{y}{(x^2+y^2)^{3/2}}}$
ただし
$x(0)=1,\ \ \ y(0)=0,\ \ \ \dot{x}(0)=0,\ \ \ \dot{y}(0)=1$
とする。この運動方程式は、どのような運動を表す方程式だろうか。また、上記の初期条件を与えると、物体はどのような運動をするだろうか。
＊ヒント＊ 上式は、以下のようなベクトル方程式に書き直すことができます。
$\displaystyle{\frac{d^2 {\bf r}}{dt^2}} =-\displaystyle{\frac{1}{r^2}}\displaystyle{\frac{{\bf r}}{r}}$
ただし、${\bf r}=(x,y)$、$r=(x^2+y^2)^{1/2}$です。
上式は、$x=y_1,\ \ y=y_2,\ \ \dot{x}=y_3,\ \ \dot{y}=y_4$と置き換えると、4次の連立常微分方程式になります。