Main Content

このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

偏微分方程式の求解

"偏微分方程式" (PDE) では、解かれる関数がいくつかの変数に依存しており、それぞれの変数に対する偏導関数を微分方程式に含めることができます。偏微分方程式は、波形、熱流量、流体分散、および経時的に変化する空間的挙動を伴う他の現象をモデル化するのに便利です。

MATLAB で解くことができる PDE のタイプ

MATLAB® PDE ソルバー pdepe は、1 つの空間変数 x と時間 t をもつ PDE 系の初期-境界値問題を解きます。これらは、経時的にも変化する 1 変数の ODE と考えることができます。

pdepe は、求解する 1 次元の方程式に略式の分類を使用します。

  • 時間微分をもつ方程式は "放物型"。例として、熱方程式があります。ut=2ux2

  • 時間微分のない方程式は "楕円型"。例として、ラプラス方程式があります。2ux2=0

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 の系を解きます。

c(x,t,u,ux)ut=xmx(xmf(x,t,u,ux))+s(x,t,u,ux).

方程式には次の特性があります。

  • これらの PDE は t0 ≤ t ≤ tf かつ a ≤ x ≤ b です。

  • 空間区間 [a, b] は有限でなければならない。

  • m は 0、1、2 のいずれかで、それぞれ"スラブ平面"、"円柱"、"球面" での対称に対応する。m > 0 の場合、a ≥ 0 でなければならない。

  • 係数 f(x,t,u,ux) は流束項、s(x,t,u,ux) はソース項である。

  • 流束項は偏導関数 ∂u/∂x に依存していなければならない。

時間に関連する偏導関数の結合は、対角行列 c(x,t,u,ux) による乗算に限定されます。この行列の対角要素は、ゼロまたは正になります。ゼロの要素は楕円型方程式に対応し、その他の要素は放物型方程式に対応します。少なくとも 1 つの放物型方程式は存在しなければなりません。放物型方程式に対応する c の要素は、x の孤立値がメッシュ点 (解が評価される点) である場合、その値でゼロにすることができます。物質の境界による c や s の不連続は、メッシュ点が各境界に置かれている場合に使われます。

解法プロセス

pdepe を使って PDE を解くには、c、f、および s に対する方程式の係数、初期条件、境界での解の動作、および解を評価するメッシュ点を定義しなければなりません。関数呼び出し sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan) はこの情報を使用して、指定したメッシュでの解を計算します。

  • m は対称定数。

  • pdefun は求解される方程式を定義。

  • icfun は初期条件を定義。

  • bcfun は境界条件を定義。

  • xmesh は x の空間値のベクトル。

  • tspan は t の時間値のベクトル。

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 は方程式 ut=2ux2 を定義します。複数の方程式がある場合、c、f、および s は、その各要素が 1 つの方程式に対応するベクトルです。

初期条件

初期時刻 t = t0 において、すべての x に関する解要素は、次の形式の初期条件を満たします。

u(x,t0)=u0(x).

MATLAB では、次の形式の関数で初期条件をコード化できます。

function u0 = icfun(x)
u0 = 1;
end
この場合、u0 = 1u0(x,t0) = 1 という初期条件を定義します。複数の方程式がある場合、u0 は各要素が 1 つの方程式の初期条件を定義するベクトルとなります。

境界条件

境界 x = a または x = b において、すべての t に対し、解要素は次の形式の境界条件を満たします。

p(x,t,u)+q(x,t)f(x,t,u,ux)=0.

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
pLqL は左境界の係数で、pRqR は右境界の係数です。この場合、bcfun が境界条件を定義します。

uL(xL,t)=0uR(xR,t)=1

複数の方程式がある場合、出力 pLqLpR、および qR は、その各要素が 1 つの方程式の境界条件を定義するベクトルです。

積分のオプション

MATLAB の PDE ソルバーの既定積分プロパティは、一般的な問題を扱えるように選択されています。場合によっては、これらの既定値をオーバーライドすることでソルバーのパフォーマンスを向上させることができます。これを行うには、odeset を使用して options 構造体を作成します。その後、構造体を最後の入力引数として pdepe に渡します。

sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options)

基礎となる ODE ソルバー ode15s のオプションのうち、次の表に示すもののみを pdepe で使用できます。

カテゴリ

オプション名

誤差制御

RelTol, AbsTol, NormControl

ステップ サイズ

InitialStep, MaxStep

イベント ログ

Events

解の計算

pdepe で方程式を解いた後、MATLAB は解を 3 次元配列 sol として返します。ここで、sol(i,j,k)t(i)x(j) で評価された解の k 番目の要素を含みます。一般に、k 番目の解要素は、コマンド u = sol(:,:,k) を使って抽出できます。

指定する時間メッシュは出力目的にのみ使用され、ソルバーの内部タイム ステップには影響しません。しかし、指定する空間メッシュは解の質と速度に影響する場合があります。方程式を解いた後、pdeval を使用して、pdepe から返された解の構造体を異なる空間メッシュで評価することができます。

例: 熱方程式

放物型 PDE の例として、1 次元の熱方程式があります。

ut=2ux2.

この方程式は、0xL かつ t0 のときの熱散逸を記述します。目標は、温度 u(x,t) を求めることです。温度の初期値は非ゼロの定数であるため、初期条件は次のようになります。

u(x,0)=T0.

また、温度は左境界ではゼロ、右境界では非ゼロであるため、境界条件は次のようになります。

u(0,t)=0,

u(L,t)=1.

この方程式を MATLAB® で解くには、方程式、初期条件、境界条件をコード化し、適切な解のメッシュを選択してからソルバー pdepe を呼び出す必要があります。(この例のように) 必要な関数をファイルの最後にローカル関数として含めることも、あるいは個別の名前付きファイルとして MATLAB パス上のディレクトリに保存することもできます。

方程式のコード化

方程式をコード化する前に、それが pdepe ソルバーで想定される形式であることを確認する必要があります。

c(x,t,u,ux)ut=x-mx(xmf(x,t,u,ux))+s(x,t,u,ux).

この形式では、熱方程式は次のようになります。

1ut=x0x(x0ux)+0.

したがって、係数の値は以下のようになります。

  • m=0

  • c=1

  • f=ux

  • s=0

m の値は引数として pdepe に渡されますが、他の係数は方程式の関数で次のようにエンコードされます。

function [c,f,s] = heatpde(x,t,u,dudx)
c = 1;
f = dudx;
s = 0;
end

(メモ: 関数はすべて例の最後にローカル関数として含まれます。)

初期条件のコード化

熱方程式の初期条件関数は定数値を u0 に代入します。この関数は、x への入力を使用しなくても受け入れなければなりません。

function u0 = heatic(x)
u0 = 0.5;
end

境界条件のコード化

以下は pdepe ソルバーで想定される境界条件の標準形式です。

p(x,t,u)+q(x,t)f(x,t,u,ux)=0.

この形式で記述すると、この問題の境界条件は次のようになります。

u(0,t)+(0f)=0,

(u(L,t)-1)+(0f)=0.

したがって、pq の値は次のようになります。

  • pL=uL, qL=0.

  • pR=uR-1, qR=0.

対応する関数は次のようになります。

function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t)
pl = ul;
ql = 0;
pr = ur - 1;
qr = 0;
end

解のメッシュの選択

点 20 個の空間メッシュと点 30 個の時間メッシュを使用します。解は急速に定常状態に達するため、出力でのこの動作を捉えるために t=0 付近の時間点の間隔はより狭くなります。

L = 1;
x = linspace(0,L,20);
t = [linspace(0,0.05,20), linspace(0.5,5,10)];

方程式の求解

最後に、対称性 m、PDE 方程式、初期条件、境界条件、および xt のメッシュを使用して方程式を解きます。

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')

Figure contains an axes object. The axes object with title Heat Equation for 0 less equals x less equals 1 and 0 less equals t less equals 5, xlabel Distance x, ylabel Time t contains an object of type surface.

ローカル関数

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 サンプル ファイルのリストが記載されています。

サンプル ファイル

説明

リンクの例

pdex1

単純な PDE を使って定式化、計算、解のプロットを説明します。

単一の PDE の求解

pdex2

不連続を含む問題を解きます。

不連続性をもつ PDE の求解

pdex3

偏導関数の値の計算が必要な問題を解きます。

PDE の求解と偏導関数の計算

pdex4

積分区間の両端で境界層をもち、小さな t で急激に変化するような解をもつ、2 式から構成される偏微分方程式系を解きます。

PDE 系の求解

pdex5

初期条件としてステップ関数をもつ PDE 系を解きます。

初期条件のステップ関数をもつ 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.

参考

| | | |

関連するトピック