常微分方程式の解法
Mathematicaはすべての常微分方程式を解けるわけではないが、かなり
多くの問題をとくことができる。
まず次の例からはじめる。
DSolve[y'[x] == a y[x],y[x],x]
DSolve[y''[x] == a y[x],y,x]
DSolve[y''[x] == a y[x],y[x],x]
ここで、微分方程式を解く命令はDSolveである。
この命令に続いて、方程式
y'[x] == a y[x]、未知関数y[x]
あるいはy、独立変数xを順に書く。
未知関数の書き方は二通りあるが、いずれにしても解はもとまる。
ただし、関数として、解を求めたいのであれば、yを用いるのがよい。
y[x]を未知関数として指定すれば、関数y[x]をもとめるが、yと指定すれば純関数の形で
解を求める。この純関数の概念は数学上の関数の概念に近い。
y'[x]はy[x]の一階微分をあらわす。それ以外の２階微分についても同様である。
この時、解は通常任意定数C[1], C[2]などを含んだ形の一般解が
もとまる。
基本的にはこの形で解けるものはすべて解の一般解をあたえる。
常微分方程式の初期値問題
次の例はいわゆる初期値問題といわれるものである。 これはある時刻(x=0とすることが多い)に未知関数の値とその微分係数をあたえてその後の 解を決定するものである。この時、括弧の中に方程式と初期値が与えられている。
DSolve[{y'[x] == a y[x],y[0] == 1},y[x],x]
DSolve[{y'[x] == - y[x],y[0] == 1},y[x],x]
DSolve[{x^2 y''[x] - 3 x y'[x] + 4 y[x] == 0,y[0] == 0,
y'[0]==1},y,x]
以下でこれ以外の例をあたえる。各自確認して欲しい。
DSolve[y''[x] - 2y'[x] + 10y[x] == 0,y,x]
DSolve[y''[x] - 2y'[x] + y[x] == 0,y,x]
DSolve[x^2 y''[x] - 3 x y'[x] + 4y[x] == 0,y,x]
DSolve[y'''[x] - y[x] == 0,y,x]
DSolve[y''''[x] + 3y''[x] - y[x] == 0,y,x]
DSolve[x^3 y'''[x] - 3 x^2 y''[x] - y[x] == 0,y,x]
DSolve[x^3 y'''[x] + x^2 y''[x]-2x y'[x] +2 y[x] == 0,y,x]
DSolve[y'''[x] - 2 y''[x]- y'[x] + 2 y[x] == 2x^2 - 6x + 4,y,x]
DSolve[y''[x] - y'[x] + 2 y[x] == 2Exp[x],y,x]
DSolve[y''[x] - y'[x] + 2 y[x] == 2Exp[x] + Sin[x],y,x]
DSolve[y''[x] + 2y'[x] + 5 y[x] == 16Exp[x] + Sin[2x],y,x]
DSolve[y'''[x] - y''[x]- 4y'[x] +4 y[x] == 6 Exp[-x] ,y,x]
DSolve[y'''[x] - y''[x]- 4y'[x] +4 y[x] == 6 Exp[-x]+Sin[x] ,y,x]
DSolve[y''[x] - 4y'[x] + 4 y[x] == Exp[-2x]/x^2 ,y,x]
DSolve[y''[x] + 4y'[x] + 4 y[x] == Exp[-2x]/x^2 ,y,x]
DSolve[x^2 y''[x]- x y'[x] == 2 Exp[x] x^3 ,y,x]
数値解法
方程式の一般解あるいは初期値問題の解が関数の形で求まらない時にもその 数値的な解を求めることはできる。これは解のグラフを描いてその性質を調べたりするのに役立つ。 数値的に微分方程式を解くためにはNDSolve命令を用いる。次の例はさらにグラフを描いている。 ここで/.は置換命令である。すなわち、y[x]をsolutionで置換する。 Evaluate命令はy[x]を評価する。
solution =NDSolve[{y''''[x] + 3y''[x] - y[x] == 0, y[0]==0, y'[0]==0, y''[0]==1, y'''[0]==0},y,{x,-5,5}];
Plot[Evaluate[y[x] /. solution], {x,-5,5} ]
この方法を用いて上の方程式の解のグラフを描くことができる。
solution =NDSolve[{y'''[x] - y''[x]- 4y'[x] +4 y[x] == 6 Exp[-x]+Sin[x] ,
y[0]==0, y'[0]==0, y''[0]==1},y,{x,-5,5}];
Plot[Evaluate[y[x] /. solution], {x,-5,5} ]
方程式が特異性を持っているときはxの範囲に注意する必要がある。たとえば次の例では x=0は特異性なので解の範囲からのぞいておく必要がある。
solution =NDSolve[{x^2 y''[x]- x y'[x] == 2 Exp[x] x^3, y[1]==0, y'[1]==1},y,{x,1,2}];
Plot[Evaluate[y[x] /. solution], {x,1,2} ]
非線型方程式を解く
非線型方程式の厳密解を求めるのは可積分の場合を除いて難しいが、数値的に解くのは 可能である場合がある。このようないくつかの場合について述べる。 まず、調和振動子の場合を考える。
solution =NDSolve[{y''[x] +y[x] == 0, y[0]==0, y'[0]==1},y,{x,-6,6}];
Plot[Evaluate[y[x] /. solution], {x,-6,6} ]
このグラフは単振動である。これに非線形の外力項を付け加える。するとグラフは次のようになる。
solution =NDSolve[{y''[x] +y[x] -y[x]^2 == 0, y[0]==0, y'[0]==1},y,{x,-4,3.5}];
Plot[Evaluate[y[x] /. solution], {x,-4,3.5} ]
ここで注意することは解の存在範囲はかならずしもすべてのxを指定することはできない。
いまのばあい、x=4付近では解が急速の増加していることがわかる。すなわち、解は爆発すると
思われる。
これに対して、非線形項の符号を逆にしてみる。このときは次のようになる。
solution =NDSolve[{y''[x] +y[x] + y[x]^2 == 0, y[0]==0, y'[0]==1},y,{x,-3.5,4}];
Plot[Evaluate[y[x] /. solution], {x,-3.5,4} ]
この場合もxの絶対値が大きくなると解は爆発することがわかる。このように非線型性は解の爆発を起こすことがある。
次の例をみる。
solution1 =NDSolve[{y''''[x] + 3y''[x] - y[x] == 0, y[0]==0, y'[0]==0, y''[0]==1, y'''[0]==0},y,{x,-5,5}];
Plot[Evaluate[y[x] /. solution1], {x,-5,5} ]
これに非線型性を与えてみる。グラフの形に大きな変化はみられない。
solution2 =NDSolve[{y''''[x] + 3y''[x] - y[x]^2 == 0, y[0]==0, y'[0]==0, y''[0]==1, y'''[0]==0},y,{x,-5,5}];
Plot[Evaluate[{y[x] /. solution2, u(x) /. solution1}], {x,-5,5} ]
xの範囲を大きくとってみる。解の様子はかなり異なっている。解の爆発が観察できる。
solution1 =
NDSolve[{y''''[x] + 3y''[x] - y[x] == 0, y[0] == 0, y'[0] == 0,
y''[0] == 1, y'''[0] == 0}, y, {x, -10, 10}];
Plot[Evaluate[y[x] /. solution1], {x, -10, 10}];
solution2 =
NDSolve[{y''''[x] + 3y''[x] - y[x]^2 == 0, y[0] == 0, y'[0] == 0,
y''[0] == 1, y'''[0] == 0}, y, {x, -10, 10}];
Plot[Evaluate[y[x] /. solution2], {x, -10, 10}]
別の非線型性を与えてみよう。この場合にははっきりとした爆発現象は観察できない。
solution1 =
NDSolve[{y''''[x] + 3y''[x] - y[x] == 0, y[0] == 0, y'[0] == 0,
y''[0] == 1, y'''[0] == 0}, y, {x, -10, 10}];
Plot[Evaluate[y[x] /. solution1], {x, -10, 10}];
solution2 =
NDSolve[{y''''[x] + 3y''[x] - y[x] +y'[x]^2 == 0, y[0] == 0, y'[0] == 0,
y''[0] == 1, y'''[0] == 0}, y, {x, -10, 10}];
Plot[Evaluate[y[x] /. solution2], {x, -10, 10}]
そこでxのとりうる範囲をもっと広げてみる。 この場合にも解の爆発がおきている。
solution1 =
NDSolve[{y''''[x] + 3y''[x] - y[x] == 0, y[0] == 0, y'[0] == 0,
y''[0] == 1, y'''[0] == 0}, y, {x, -10, 10}];
Plot[Evaluate[y[x] /. solution1], {x, -10, 10}];
solution2 =
NDSolve[{y''''[x] + 3y''[x] - y[x] +y'[x]^2 == 0, y[0] == 0, y'[0] == 0,
y''[0] == 1, y'''[0] == 0}, y, {x, -14, 14}];
Plot[Evaluate[y[x] /. solution2], {x, -14, 14}]
非線型性の次数を変えてみる。するとつぎのようになる。
solution1 =
NDSolve[{y''''[x] + 3y''[x] - y[x] == 0, y[0] == 0, y'[0] == 0,
y''[0] == 1, y'''[0] == 0}, y, {x, -10, 10}];
Plot[Evaluate[y[x] /. solution1], {x, -10, 10}];
solution2 =
NDSolve[{y''''[x] + 3y''[x] - y[x] + Sqrt[1+ y'[x]^2] == 0, y[0] == 0, y'[0] == 0,
y''[0] == 1, y'''[0] == 0}, y, {x, -14, 14}];
Plot[Evaluate[y[x] /. solution2], {x, -14, 14}]