DSolveの使用:ユーザガイド
はじめに
このチュートリアルの目的は,DSolveでさまざまな種類の問題を解くための自己充足的作業ガイドを提供することである.
DSolveを使用する最初のステップは,問題を正確に設定することである.次はDSolveを使って解の式を得ることである.一旦解が見付かると,記号的あるいは数値的方法を使って検証することができる.また,Wolframシステム関数のPlot,Plot3D,ContourPlot等を使ってプロットすることもできる.プロットによって,閉形式の式では明確ではない解についての情報が分かることがよくある.
問題に境界条件が指定されていない場合,DSolveの出力は任意パラメータを含む一般解の形式となる.GeneratedParametersオプションを使うと,このような任意パラメータにラベルを付けることができる.
多くの分野では,ロジスティック方程式の成長率等,微分方程式に記号パラメータが含まれる.微分方程式には,直前の計算から生じた機械数のような非厳密な数値が含まれることもある.DSolveでは記号パラメータも非厳密な数値も使うことができるが,その存在を認識した上で解を正確に解釈するのがよい.
計算中にDSolveが仮定したり,不具合に遭遇したりすると,問題を描写する警告メッセージを発する.このメッセージは通常無視してもよいが,問題に対して与えられた答に重大な制約が存在することを指摘する場合もある.
問題の文は,曖昧さがないように分析するのがよい.つまり,問題が正しい形式であることを確認するのである.そうすることにより,DSolveから意味のある答が得られる.
問題の設定
DSolveに与える第1引数は微分方程式,第2引数は未知の関数であり,最後の引数で独立変数を識別する.
sol = DSolve[y'[x] + 5y[x] == 1, y[x], x]DSolveの出力は微分方程式の解のリストになる.方程式の中には複数解を持つものもあるので,別のリストが必要になる.ここでは,方程式は一階線形なので,解はy[x]->
+-5 x C[1]の1つしかない.初期条件が指定されていないので,解には1つの未知定数C[1]がある.部分指定を使うと,解のリストから解を抽出することができる.
m = sol[[1]]この解の形式は
自身を見付けるのには便利であるが,
の導関数およびある点における
の値を見付けるのには不便である.
y[x] /. my'[x] /. my[1] /. m解を今後の計算で使う場合,
ではなく
を使って未知関数を指定するのが一番よい.これにより,Function[x,expr]という形式の純関数を使って解が与えられる.
sol = DSolve[y'[x] + 5y[x] == 1, y, x]解が純関数の形式である場合,指定された点における
の導関数の式と
の値の式を見付けることができる.
m = sol[[1]]y[x] /. my'[x] /. my[1] /. m問題に複数の解がある場合,解のリストから個々の解を選ぶことも,直接そのリストを使うこともできる.
sol = DSolve[y'[x] ^ 2 == x + 11, y, x]y[x] /. sol[[1]]y[x] /. sol[[2]]y[x] /. sol方程式系を解くためには,DSolveの第1引数は方程式のリストに,第2引数は未知関数のリストにしなければならない.
sol = DSolve[{x'[t] == y[t] + z[t], y'[t] + z[t] - x[t] == 0, z'[t] + y[t] == x[t]}, {x, y, z}, t]系のそれぞれの解は,未知関数の置換規則のリストになっている.未知関数の式は,前の例と同様にして抽出することができる.
{x[t], y[t], z[t]} /. sol[[1]]//Simplify問題に初期条件が指定されている場合は,未知定数のいくつか,またはすべてが削除される.
DSolve[{y''[x] + y[x] == 5, y[0] == 1, y'[0] == 7}, y[x], x]DSolve[{y''[x] + y[x] == 5, y[0] == 1}, y[x], x]偏微分方程式では,DSolveの第3引数は方程式の独立変数のリストである.
DSolve[D[u[x, y], x] + Sin[x] * D[u[x, y], y] == 8, u, {x, y}]DSolve[{x'[t] + y[t] == Sin[t], x[t] + y[t] == 1, x[0] == 4}, {x, y}, t]問題の解が常に陽形式で与えられるとは限らない.その場合,解は未評価のSolveオブジェクト,あるいは,InverseFunctionを使って与えられる.
sol = DSolve[y'[x] + y[x] ^ 3 + y[x] ^ 2 == 1, y[x], x]sol[[1]]sol = DSolve[Derivative[2][y][x] + y[x] * Derivative[1][y][x] ^ 4 == 0, y, x]soly = y[x] /. sol[[1]]implicitequation = ( Head[soly][[1]][y[x]] == soly[[1]] )解の検証
DSolveで与えられた解は,さまざまな方法を使って検証することができる.その中でも最も簡単な方法は,解を方程式に代入することである.結果がTrueとなれば,解は有効である.
eqn = y'[x] == x;sol = DSolve[eqn, y, x]eqn /. soleqn = {y'[x] == x, y[0] == 1};sol = DSolve[eqn, y, x]eqn /. sol時として,代入の結果がTrue,Falseよりも複雑になることがある.そのような場合は,Simplifyを使って代入の結果を簡約し,検証することができる.簡約された結果がTrueであれば,解は有効である.
eqn = y''[x] + 5 * y'[x] + 6y[x] == 1;sol = DSolve[eqn, y, x]eqn /. solSimplify[eqn /. sol]PDE = D[u[x, y], x] + Sin[x] * D[u[x, y], y] == x ^ 2;sol = DSolve[PDE, u, {x, y}]Simplify[PDE /. sol]方程式に特殊関数が含まれている場合は,FullSimplifyを使って解を検証する必要があるかもしれない.
eqn = x * y''[x] + y'[x] - y[x] == 1;sol = DSolve[eqn, y, x]FullSimplify[eqn /. sol]解が大きい場合,あるいはSimplifyおよびFullSimplifyでは解が検証できない場合は,RandomRealまたはRandomComplexを使って問題の変数とパラメータすべての値を生成することにより,数値的チェックを行うことができる.その場合,乱数値の集合をいくつか使ってチェックを繰り返した方がよい.
Clear[a, y, x, r]eqn = y''[x] - (a * x ^ 6 + x ^ 2) * y[x];sol = DSolve[eqn == 0, y, x]Union[Flatten[Table[Chop[eqn /. sol /. {x -> RandomReal[], a -> RandomReal[], C[1] -> RandomReal[], C[2] -> RandomReal[]}], {i, 1, 10}]]]数値チェックでは,確実に解を検証することはできないが,より高い精度を使ったり,変数に複素数値を許可したりすることにより,さらに厳密なチェックを行うことができる.
Chop[eqn /. sol /. {x -> RandomReal[ {1, 2}, WorkingPrecision -> 20], a -> RandomReal[ {1, 2}, WorkingPrecision -> 20], C[1] -> RandomReal[ {1, 2}, WorkingPrecision -> 20], C[2] -> RandomReal[ {1, 2}, WorkingPrecision -> 20]}]r := RandomComplex[]Chop[N[eqn /. sol /. {x -> r, a -> r, C[1] -> r, C[2] -> r}]]上記の方法は,解が陽形式で求められるときだけ使用できる.最後の例では,一階常微分方程式の解が陰形式で与えられるときの検証方法を示す.
eqn = Derivative[1][y][x] + 2 * x y[x]^2 + y[x]^3;sol = DSolve[eqn == 0, y, x]sol[[1]]Simplify[Solve[D[sol[[1]], x], y'[x]]]解のプロット
DSolveで与えられる解のプロットは,解の性質(例えば,振動するものかどうか)についての有用な情報を与えることがある.また,グラフの形状が理論から,あるいは,微分方程式に関連するベクトル場のプロットから判明する場合は,プロットは解の検証方法としての役割を果たすことができる.異なるWolfram言語のグラフィックス関数を使用する例を以下に挙げる.
sol = DSolve[y'[x] + x * y[x] == Cos[x ^ 2], y, x]Plot[Evaluate[y[x] /. sol /. {C[1] -> 1}], {x, -7, 7}, PlotRange -> All]sol = DSolve[{y''[x] / y[x] == -4 * Exp[-x / 4], y[0] == 1, y'[0] == 1 / 2}, y, x]Plot[Evaluate[y[x] /. sol], {x, 0, 30}]sol = DSolve[{y'[x] ^ 2 == x - x ^ 3, y[0] == 1}, y, x]Plot[Evaluate[y[x] /. sol], {x, 1 / 100, 9 / 10}]sol = DSolve[Derivative[1][y][x] == -5 y[x]^2 - (y[x]^3/x), y[x], x]expr = C[1] /. Solve[sol[[1]], C[1]][[1]] /. {y[x] -> y}ContourPlot[expr, {x, -0.4, -0.1}, {y, 1., 1.8}]sol = DSolve[{x'[t] + t * y[t] == 0, 2y'[t] - 3x[t] == 0, x[0] == 1, y[0] == 3}, {x, y}, t]Plot[Evaluate[{x[t], y[t]} /. sol], {t, -1, 10}, WorkingPrecision -> 20]ParametricPlot[Evaluate[{x[t], y[t]} /. sol], {t, -1, 7}, WorkingPrecision -> 20]sol = DSolve[{x''[t] + 3 * y[t] == UnitStep[t], x[t] - 5 * y[t] == t ^ 2, x[0] == 1, x'[0] == 1}, {x, y}, t]Plot[Evaluate[{x[t], y[t], x[t] - 5 * y[t]} /. sol], {t, -2, 2}]sol = DSolve[D[u[x, y], x] + x ^ 2 * D[u[x, y], y] == Exp[x], u, {x, y}]Plot3D[Evaluate[u[x, y] /. sol /. {C[1][t_] -> Sin[3 * t]}[[1]]], {x, -2, 2}, {y, -2, 2}, Mesh -> False]生成されたパラメータ
微分方程式の一般解にはC[1],C[2]等のラベルの付く未定係数が含まれる.
DSolve[y'[x] + y[x] == 1, y[x], x]未定係数の名前を変更するためには,GeneratedParametersオプションを使う.
DSolve[y'[x] + y[x] == 1, y[x], x, GeneratedParameters -> P]パラメータCは,異なる定数C[i]を生成する指標の集合に作用する純関数と考えなければならない.
parameter = C[#]&;indexset = {1, 2, 3, 4};parameter[indexset[[1]]]parameter[indexset[[3]]]内部的に純関数を使用することにより,DSolveは高階常微分方程式および常微分方程式系にあるC[i]の引数 i を正確に増加させることができる.
DSolve[y''[x] + y[x] == 1, y[x], x, GeneratedParameters -> (const[#]&)]1(デフォルト)以外の任意の値でパラメータの指標付けを開始したい場合は,純関数を使うと特に便利である.
DSolve[y''[x] + y[x] == 1, y[x], x, GeneratedParameters -> (const[1 + #]&)]パラメータ指標に添え字付き変数やその他のスタイルを使って解を表示すると便利なことがある.
DSolve[y''[x] + y[x] == 1, y[x], x, GeneratedParameters -> (Subscript[c, #1]&)]Module変数を使うと,各DSolveの呼出しで一意であるパラメータ名を得ることができる.
DSolve[y''[x] + y[x] == 1, y[x], x, GeneratedParameters -> Module[{C}, C[#]&]]DSolve[y''[x] + y[x] == 1, y[x], x, GeneratedParameters -> Module[{C}, C[#]&]]ホロノミック解
多項式係数または有理関数係数を持つ線形常微分方程式の解はホロノミック関数DifferentialRootで表すことができる.
sol = DSolve[y'''[x] - x ^ 2 y'[x] - y[x] - Cos[x] == 0 && y[0] == 0 && y'[0] == 1 && y''[0] == 0, y, x]Plot[y[x] /. sol[[1]], {x, 0, 1}]オプションMethod"Holonomic"を使うと,DSolveが線形微分方程式にホロノミック解を返すよう強制する.
DSolve[{y''[x] + y[x] == 0, y[0] == 0, y'[0] == 1}, y[x], x]DSolve[{y''[x] + y[x] == 0, y[0] == 0, y'[0] == 1}, y[x], x, Method -> "Holonomic"]特異解
DSolveはデフォルトでは線形または非線形の常微分方程式に対して,任意のパラメータに依存する一般解を返す.非線形微分方程式によっては特異解を持つものもある.これらの特異解は特定の値を一般解の任意の定数に割り当てることで得ることはできないが,力学系等の研究では便利である.
DSolveのオプションIncludeSingularSolutionsTrueを使うと,非線形常微分方程式について一般解と一緒に特異解が返される.
DSolve[y'[x] ^ 2 - y'[x] * x + y[x] == 0, y[x], x]DSolve[y'[x] ^ 2 - y'[x] * x + y[x] == 0, y, x, IncludeSingularSolutions -> True]Plot[Table[x c - c^2, {c, -10, 10}]//Evaluate, {x, -50, 50}, PlotRange -> All]DSolve[y'[x] == y[x] (1 - y[x]), y[x], x]DSolve[y'[x] == y[x] (1 - y[x]), y, x, IncludeSingularSolutions -> True]Plot[Table[(E^x/E^x + E^c), {c, -10, 10}]//Evaluate, {x, -20, 20}, PlotRange -> All]仮定
常微分方程式は,パラメータ・変数の特定のタイプや範囲によって異なる解を持つ場合がある.DSolveのAssumptionsを使うと,パラメータおよび変数のタイプや範囲を指定して,必要な解を選ぶことができる.
DSolve[{y''[x] + λ y[x] == 0, y[0] == 0, y[π] == 0}, y[x], x]DSolve[{y''[x] + λ y[x] == 0, y[0] == 0, y[π] == 0}, y[x], x, Assumptions -> 0 < λ < 20]eqn = {y'[x] == Abs[x + Abs[x]] y[x], y[0] == 1};sol = DSolve[eqn, y, x, Assumptions -> x∈Reals]記号パラメータと非厳密な数値
- 唯一の変数が独立・従属変数である方程式.従って,DSolveの第1引数に現れる変数がすべて第2,第3引数にも現れる.
- 質量あるいはバネ定数等,他の記号的数量を含む方程式.この場合の解は,独立変数,従属変数,付加的な記号的パラメータに依存する.
DSolve[y''[x] - 8 * x * y[x] == 0, y, x]sol = DSolve[{y''[x] + k ^ 2 * y[x] == x ^ 2, y[0] == 3k + 1, y'[0] == k}, y, x]DSolveは上記のどちらのタイプの方程式も扱うことができる.2つ目のタイプの方程式では,可能なパラメータ値すべてについての解が使えるようにしておくと非常に便利である.
Plot[Evaluate[Table[y[x] /. sol /. {k -> i}, {i, 1, 5, 2}]], {x, -7, 7}]記号パラメータがあると,かなり複雑な出力になる可能性がある.
eqn = y''[x] - c * x ^ a * y[x] ;sol = DSolve[eqn == 0, y, x]sol1 = y[x] /. sol[[1]] /. {a -> 0, c -> 4}解がパラメータの大部分の値については有効であるが,すべての値については有効でないという場合がある.
sol = DSolve[{y'[x] == x, y[0] == 1 / a}, y, x]sol = DSolve[{y'[x] == x, y[0] == k}, y, x]つまり,記号パラメータを持つ微分方程式が解けるということは,DSolveのような記号ソルバのパワフルで不可欠な機能なのである.しかし,下記の点には注意しなければならない.
- 解が複雑になることがある.そのような計算には,かなりの時間とメモリが必要となることが多い.
- パラメータが例外的な値のときは,答が有効ではない場合がある.
- パラメータが特殊値のときは,解の検証を記号的に行うと簡単であるが,一般的には数値検証法の方がよい.
Wolfram言語の数には,無限精度,機械精度,任意精度の3種類がある.無限精度数は「厳密数」といわれ,残りの2つは不完全な情報を表すため「非厳密数」と呼ばれる.
exactpi = Pi;Precision[exactpi](machinepi = N[Pi])//InputFormPrecision[machinepi]highprecisionpi = N[Pi, 30]Precision[highprecisionpi]DSolveは記号ソルバなので,使用されるアルゴリズムは主に入力が厳密であるという前提に基づいている.しかし,DSolveは非厳密な数値を含む方程式も通常通りに扱う.
DSolve[x'[t] == 3. x[t], x, t]DSolve[x'[t] == 3.`40 * x[t], x, t]非厳密な入力が生じるのは,例えば,方程式の係数が前の計算から導かれており,近似的にのみ分かっている場合等である.そのような場合,方程式を厳密形式に変換することは,計算を著しく遅くするので,実践的ではない.
n = 5;
Subscript[x, 0][t_] := 0;
Subscript[x, n][t_] := 1;
eqns = Table[{Subscript[x, i]'[t] - ( Subscript[x, i + 1][t] - 2 Subscript[x, i][t] + Subscript[x, i - 1][t]), Subscript[x, i][0] - (1 / n)}, {i, n - 1}];
vars = Table[Subscript[x, i], {i, n - 1}];
sol = DSolve[Map[# == 0&, Flatten[eqns]], vars, t];//TimingLeafCount[sol]eqns /. sol /. {t -> RandomReal[ {0, 1}, WorkingPrecision -> 200]}//N//Chop//Flatten//Unionn = 5;
Subscript[x, 0][t_] := 0.;
Subscript[x, n][t_] := 1;
eqns = Table[{Subscript[x, i]'[t] - ( Subscript[x, i + 1][t] - 2 Subscript[x, i][t] + Subscript[x, i - 1][t]), Subscript[x, i][0] - (1 / n)}, {i, n - 1}];
vars = Table[Subscript[x, i], {i, n - 1}];
sol = DSolve[Map[# == 0&, Flatten[eqns]], vars, t];//TimingLeafCount[sol]eqns /. sol /. {t -> RandomReal[ {0, 1}] }//N[#]&//Chop//Flatten//Unionこのように,DSolveのような記号関数の中でも,非厳密数を使い続けたほうがよい場合がよくある.しかし,そのような場合に得られる解にはある程度の数値誤差が含まれるため,注意深くチェックしなければならない.問題が大きすぎない(方程式が5つ未満である等)場合は,関数Rationalizeを使って入力を厳密形に変換しなければならない.
DSolve[{x'[t] == 0.0001 * x[t], x[0] == 3.07}, x[t], t]DSolve[Rationalize[{x'[t] == 0.0001 * x[t], x[0] == 3.07}, 0], x[t], t]問題が適切であるかどうか
DSolveは初期条件か境界条件が指定されていない場合,問題の一般解を返す.
DSolve[y'[x] == 1 - y[x], y, x]しかし,初期条件か境界条件が指定されている場合は,DSolveの出力は内在する微分方程式と与えられた条件の両方を満足しなければならない.
eqns = {y'[x] == 1 - y[x], y[3] == 5};sol = DSolve[eqns, y, x]eqns /. solこのような場合,DSolveへ与える問題が妥当であるかどうかをチェックしてみるとよい.つまり,問題が正しい形式であるかどうかをチェックするのである.初期値あるいは境界値問題の解が有名な関数(正則関数等)に存在することが保証される場合,解が一意の場合,解がデータに連続的に依存する場合,問題は適切であるといえる.
階常微分方程式(あるいは
個の一階方程式からなる系)と
個の初期条件があるとすると,指定された条件の下で問題は適切であることを示す標準的な存在・一意性定理がある.前の例の一階線形常微分方程式の右辺は
の多項式なので,無限に微分可能である.これは,右辺がLipschitz連続であることしか必要としないピカール(Picard)の存在・一意性定理を適用するのに十分である.
実際に生じる問題は,正しい理論的原則から導かれているため,ほとんどが適切である.しかし,注意事項として,DSolveが問題を満足する解を見付けることが難しい例を以下に挙げる.
generalsol = DSolve[{y'[x] == 1 / y[x]}, y, x]Plot[Evaluate[y[x] /. generalsol /. {C[1] -> 1}], {x, -1, 3}]DSolve[{y'[x] == 1 / y[x], y[0] == 1}, y, x]DSolve[{y''[x] + y[x] == 0, y[0] == 1, y[Pi] == 3}, y[x], x]問題に解があっても,一般解が陰形式であったり高次の超越関数を含んでいたりするために,DSolveでは見付けることができないことがある.
DSolve[y'[x] == 1 / (x - y[x]) && y[0] == 1, y, x]DSolve[x'[y] == (x[y] - y) && x[1] == 0, x, y]DSolveValueを使う
DSolveとともに,関数DSolveValueも微分方程式を解くために使うことができる.DSolveValueは常微分方程式,偏微分方程式,微分代数方程式,遅延微分方程式,積分方程式,積分微分方程式,ハイブリッド微分方程式を解くことができる.
DSolveValueの出力は従属関数の形式
または
で制御される.
DSolveValue[y'[x] + y[x] == a Sin[x], y[x], x]sol = DSolveValue[{y'[x] + y[x] == a Sin[x], y[0] == 0}, y, x]{y'[x] + y[x] == a Sin[x], y[0] == 0} /. y -> sol//SimplifyPlot[Table[sol[x] /. {a -> i}, {i, 7}], {x, -2, 10}]DSolveとは対照的に,DSolveValueの第2引数は従属および独立変数の任意の式でよい.DSolveValue[eqn,expr,x]は独立変数
を持つ常微分方程式 eqn の記号解によって決定された expr の値を与える.
DSolveValue[{y''[x] + y[x] == 0, y[0] == 0, y'[0] == 1}, y[3], x]DSolveValue[{y''[x] + y[x] == 0, y[0] == 0, y'[0] == 1}, y'[x], x]微分方程式の解が複数の分岐を持つ場合,DSolveValueは警告メッセージを生成し解のうちの一つだけを返す.
DSolveValue[{y'[x] + y[x] == Sin[x], y'[0] ^ 2 == 1}, y[x], x]DSolve[{y'[x] + y[x] == Sin[x], y'[0] ^ 2 == 1}, y[x], x]以上でDSolveを効率的に使うための基本原則についての説明を終わる.DSolveの開発中やこのドキュメントの執筆中に役立った「参考文献」のリストがある.