偏微分方程式の求解
"偏微分方程式" (PDE) では、解かれる関数がいくつかの変数に依存しており、それぞれの変数に対する偏導関数を微分方程式に含めることができます。偏微分方程式は、波形、熱流量、流体分散、および経時的に変化する空間的挙動を伴う他の現象をモデル化するのに便利です。
MATLAB で解くことができる PDE のタイプ
MATLAB® PDE ソルバー pdepe
は、1 つの空間変数 x と時間 t をもつ PDE 系の初期-境界値問題を解きます。これらは、経時的にも変化する 1 変数の ODE と考えることができます。
pdepe
は、求解する 1 次元の方程式に略式の分類を使用します。
時間微分をもつ方程式は "放物型"。例として、熱方程式があります。
時間微分のない方程式は "楕円型"。例として、ラプラス方程式があります。
pdepe
では、系の中に放物型方程式が少なくとも 1 つ必要です。つまり、系の方程式の少なくとも 1 つに時間微分が含まれていなければなりません。
pdepe
はまた、角の対称性により 1 次元問題に低次元化される、特定の 2 次元および 3 次元の問題も解きます (詳細については、対称定数 m
の引数の説明を参照してください)。
Partial Differential Equation Toolbox™ は、ディリクレとノイマンの境界条件によって、この機能を 2 次元および 3 次元の一般化問題へと拡張します。
1 次元 PDE の求解
1 次元 PDE には、時間 t と 1 つの空間変数 x に依存する関数 u(x,t) が含まれます。MATLAB PDE ソルバー pdepe
は、次の形式の 1 次元放物型および楕円型 PDE の系を解きます。
方程式には次の特性があります。
これらの PDE は t0 ≤ t ≤ tf かつ a ≤ x ≤ b です。
空間区間 [a, b] は有限でなければならない。
m
は 0、1、2 のいずれかで、それぞれ"スラブ平面"、"円柱"、"球面" での対称に対応する。m > 0 の場合、a ≥ 0 でなければならない。係数 は流束項、 はソース項である。
流束項は偏導関数 ∂u/∂x に依存していなければならない。
時間に関連する偏導関数の結合は、対角行列 による乗算に限定されます。この行列の対角要素は、ゼロまたは正になります。ゼロの要素は楕円型方程式に対応し、その他の要素は放物型方程式に対応します。少なくとも 1 つの放物型方程式は存在しなければなりません。放物型方程式に対応する c の要素は、x の孤立値がメッシュ点 (解が評価される点) である場合、その値でゼロにすることができます。物質の境界による c や s の不連続は、メッシュ点が各境界に置かれている場合に使われます。
解法プロセス
pdepe
を使って PDE を解くには、c、f、および s に対する方程式の係数、初期条件、境界での解の動作、および解を評価するメッシュ点を定義しなければなりません。関数呼び出し sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan)
はこの情報を使用して、指定したメッシュでの解を計算します。
xmesh
および tspan
ベクトルはともに 2 次元グリッドを形成し、そこで pdepe
により解が評価されます。
方程式
PDE は、pdepe
で想定される標準形式で表現しなければなりません。この形式で記述すると、係数 c、f、および s の値を読み取ることができます。
MATLAB では方程式を次の形式の関数を使ってコード化できます。
function [c,f,s] = pdefun(x,t,u,dudx) c = 1; f = dudx; s = 0; end
pdefun
は方程式 を定義します。複数の方程式がある場合、c、f、および s は、その各要素が 1 つの方程式に対応するベクトルです。初期条件
初期時刻 t = t0 において、すべての x に関する解要素は、次の形式の初期条件を満たします。
MATLAB では、次の形式の関数で初期条件をコード化できます。
function u0 = icfun(x) u0 = 1; end
u0 = 1
は u0(x,t0) = 1 という初期条件を定義します。複数の方程式がある場合、u0
は各要素が 1 つの方程式の初期条件を定義するベクトルとなります。境界条件
境界 x = a または x = b において、すべての t に対し、解要素は次の形式の境界条件を満たします。
q(x,t) は、ゼロまたは非ゼロの要素をもつ対角行列です。この境界条件は、x に対する u の偏導関数ではなく、流束の項 f で表現さていることに注意してください。また、2 つの係数 p(x,t,u) と q(x,t) のうち、p のみが u に依存します。
MATLAB では、次の形式の関数で境界条件をコード化できます。
function [pL,qL,pR,qR] = bcfun(xL,uL,xR,uR,t) pL = uL; qL = 0; pR = uR - 1; qR = 0; end
pL
と qL
は左境界の係数で、pR
と qR
は右境界の係数です。この場合、bcfun
が境界条件を定義します。複数の方程式がある場合、出力 pL
、qL
、pR
、および qR
は、その各要素が 1 つの方程式の境界条件を定義するベクトルです。
積分のオプション
MATLAB の PDE ソルバーの既定積分プロパティは、一般的な問題を扱えるように選択されています。場合によっては、これらの既定値をオーバーライドすることでソルバーのパフォーマンスを向上させることができます。これを行うには、odeset
を使用して options
構造体を作成します。その後、構造体を最後の入力引数として pdepe
に渡します。
sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options)
基礎となる ODE ソルバー ode15s
のオプションのうち、次の表に示すもののみを pdepe
で使用できます。
カテゴリ | オプション名 |
---|---|
誤差制御 | |
ステップ サイズ | |
イベント ログ |
解の計算
pdepe
で方程式を解いた後、MATLAB は解を 3 次元配列 sol
として返します。ここで、sol(i,j,k)
は t(i)
と x(j)
で評価された解の k
番目の要素を含みます。一般に、k
番目の解要素は、コマンド u = sol(:,:,k)
を使って抽出できます。
指定する時間メッシュは出力目的にのみ使用され、ソルバーの内部タイム ステップには影響しません。しかし、指定する空間メッシュは解の質と速度に影響する場合があります。方程式を解いた後、pdeval
を使用して、pdepe
から返された解の構造体を異なる空間メッシュで評価することができます。
例: 熱方程式
放物型 PDE の例として、1 次元の熱方程式があります。
この方程式は、 かつ のときの熱散逸を記述します。目標は、温度 を求めることです。温度の初期値は非ゼロの定数であるため、初期条件は次のようになります。
また、温度は左境界ではゼロ、右境界では非ゼロであるため、境界条件は次のようになります。
この方程式を MATLAB® で解くには、方程式、初期条件、境界条件をコード化し、適切な解のメッシュを選択してからソルバー pdepe
を呼び出す必要があります。(この例のように) 必要な関数をファイルの最後にローカル関数として含めることも、あるいは個別の名前付きファイルとして MATLAB パス上のディレクトリに保存することもできます。
方程式のコード化
方程式をコード化する前に、それが pdepe
ソルバーで想定される形式であることを確認する必要があります。
この形式では、熱方程式は次のようになります。
したがって、係数の値は以下のようになります。
の値は引数として pdepe
に渡されますが、他の係数は方程式の関数で次のようにエンコードされます。
function [c,f,s] = heatpde(x,t,u,dudx) c = 1; f = dudx; s = 0; end
(メモ: 関数はすべて例の最後にローカル関数として含まれます。)
初期条件のコード化
熱方程式の初期条件関数は定数値を に代入します。この関数は、 への入力を使用しなくても受け入れなければなりません。
function u0 = heatic(x) u0 = 0.5; end
境界条件のコード化
以下は pdepe
ソルバーで想定される境界条件の標準形式です。
この形式で記述すると、この問題の境界条件は次のようになります。
したがって、 と の値は次のようになります。
.
対応する関数は次のようになります。
function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t) pl = ul; ql = 0; pr = ur - 1; qr = 0; end
解のメッシュの選択
点 20 個の空間メッシュと点 30 個の時間メッシュを使用します。解は急速に定常状態に達するため、出力でのこの動作を捉えるために 付近の時間点の間隔はより狭くなります。
L = 1; x = linspace(0,L,20); t = [linspace(0,0.05,20), linspace(0.5,5,10)];
方程式の求解
最後に、対称性 、PDE 方程式、初期条件、境界条件、および と のメッシュを使用して方程式を解きます。
m = 0; sol = pdepe(m,@heatpde,@heatic,@heatbc,x,t);
解のプロット
pcolor
を使用して解の行列を可視化します。
colormap hot pcolor(x,t,sol) colorbar xlabel('Distance x','interpreter','latex') ylabel('Time t','interpreter','latex') title('Heat Equation for $0 \le x \le 1$ and $0 \le t \le 5$','interpreter','latex')
ローカル関数
function [c,f,s] = heatpde(x,t,u,dudx) c = 1; f = dudx; s = 0; end function u0 = heatic(x) u0 = 0.5; end function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t) pl = ul; ql = 0; pr = ur - 1; qr = 0; end
PDE の例とファイル
最も一般的な 1 次元 PDE の問題に対し開始点として効果的に使用できるサンプル ファイルがいくつかあります。例を確認し実行するには、[微分方程式の例] アプリを使用します。このアプリを実行するには、次のように入力します。
odeexamples
個々のファイルを編集用に開くには、次のように入力します。
edit exampleFileName.m
例を実行するには、次のように入力します。
exampleFileName
次の表には、利用可能な PDE サンプル ファイルのリストが記載されています。
サンプル ファイル | 説明 | リンクの例 |
---|---|---|
| 単純な PDE を使って定式化、計算、解のプロットを説明します。 | |
| 不連続を含む問題を解きます。 | |
| 偏導関数の値の計算が必要な問題を解きます。 | |
| 積分区間の両端で境界層をもち、小さな t で急激に変化するような解をもつ、2 式から構成される偏微分方程式系を解きます。 | |
| 初期条件としてステップ関数をもつ PDE 系を解きます。 |
参照
[1] Skeel, R. D. and M. Berzins, "A Method for the Spatial Discretization of Parabolic Equations in One Space Variable," SIAM Journal on Scientific and Statistical Computing, Vol. 11, 1990, pp. 1–32.
参考
bvp4c
| ode45
| pdepe
| odeset
| pdeval