MATLAB Function Reference

pdepe

1空間変数と時間で表される放物線偏微分方程式と楕円偏微分方程式システム(PDEs)の初期-境界値問題を解きます。

表示

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

引数

m	問題の対称性に対応したパラメータ。m が、0 の場合、スラブ、1 の場合、円筒、2 の場合、球面です。
pdefun	PDE の要素を定義する関数
icfun	初期条件を定義する関数
bcfun	境界条件を定義する関数
xmesh	tspan の中のすべての値に対して、数値解が必要とされる点を指定するベクトル [x0, x1, ..., xn] です。xmesh の要素は、x0 < x1 < ... < xn を満足する必要があります。xmesh の長さは、 3 でなければなりません。
tspan	xmesh の中のすべての値に対して、解が必要とする点を指定するベクトル [t0, t1, ..., tf] です。tspan の要素は、t0 < t1 < ... < tf を満足する必要があります。 tspan の長さは、 3 でなければなりません。
options	pdepe で使用可能な ODE ソルバで使用するいくつかのオプション、RelTol, AbsTol, NormControl, InitialStep, MaxStep です。多くの場合、これらのオプションに対するデフォルト値は、満足の得られる解を与えます。詳細は、odeset を参照してください。
p1,p2,...	pdefun, icfun, bcfun に渡されるオプションパラメータ

詳細

sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan) は、1空間変数 と時間変数 からなる放物線 PDEs と楕円 PDEs に対する初期-境界値問題を解きます。空間的な離散化の結果求まる常微分方程式は、tspan の中に指定された時刻で、近似解を得るために、積分されます。関数 pdepe は、xmesh の中で与えられるメッシュ上の解の値を戻します。

pdepe は、つぎの型の PDEs を解きます。

(1-1)

PDEs は、 と で、成り立ちます。区間 は、有限で、 は、0、1、2のいずれかで、スラブ、円筒、球面対称に対応します。 の場合、 は、 0 である必要があります。

式 1-1 で、 は、流量の項、 は、ソースの項です。時間に関する偏微分の干渉は、対角行列 との乗算に制限されます。この行列の対角要素は、ゼロ、または、正数のいずれかです。ゼロの要素は、楕円方程式に対応し、その他は、放物線の方程式に対応します。少なくとも一つは、放物線方程式が存在する必要があります。放物線方程式に対応する の要素は、 のそれらの値が、メッシュ点上に存在すれば、 の孤立点でゼロになります。物質境界による や の中の不連続は、メッシュ点が、個々のインタフェースに配置されている場合には、邪魔になりません。

と、すべての に対して、解の要素は、つぎの型の初期条件を満足します。

(1-2)

すべての と、、または、 のどちらか一方と共に、解要素は、つぎの型の境界条件を満足します。

(1-3)

の要素は、ゼロ、または、非ゼロのどちらかです。境界条件が、 でなく、流量 の項で表されることに注意してください。2つの係数に関して、 だけ、 に依存します。

sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan) のコールの中で、

• m は、 に対応します。
• xmesh(1) と xmesh(end) は、 と に対応します。
• tspan(1) と tspan(end) は、 と に対応します。
• pdefun は、項 , , (式 1-1) を計算します。つぎの型をもちます。

[c,f,s] = pdefun(x,t,u,dudx)

入力引数は、スカラ x と t 、ベクトル u と dudx で、解 と に関する偏微係数を近似します。c, f, s は、列ベクトルです。c は、行列 の対角要素をストアします(式 1-1)。

• icfun は、初期条件を計算します。つぎの型をしています。

u = icfun(x)

引数 x と共にコールすると、icfun は、列ベクトル u に、x での解要素の初期値を計算し、戻します。

• bcfun は、境界条件の と の項を計算します(式 1-3)。つぎの型をしています。

[pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t)

ul は、左境界 xl = での近似解で、ur は、右境界 xr = での近似解です。pl と ql は、xl で計算した と に対応する列ベクトルです。同様に、pr と qr は、xr に対応します。 で、 の場合、 近傍での解の有界は、流量 が、 で、ゼロになることを必要とします。pdepe は、自動的に、この境界条件を含み、pl と ql に戻される値を無視します。

pdepe は、多次元配列 sol として、解を戻します。 = ui = sol(:,:,i) は、解ベクトル の i 番目の要素への近似になります。要素 ui(j,k) = sol(j,k,i) は、 = (tspan(j),xmesh(k)) で、 を近似します。

ui = sol(j,:,i) は、時刻 tspan(j) とメッシュ点 xmesh(:) で、解の成分 i を近似します。pdeval を使って、xmesh に含まれない点で、近似とその偏微係数 を計算します。詳細は、pdeval を参照してください。

sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options) は、デフォルトの積分パラメータの代わりに、関数 odeset で作成した引数 options に設定した値を使って、上述したようにして、解きます。ODE ソルバのオプションのいくつかは、pdepe の中で使用可能です。これは、RelTol, AbsTol, NormControl, InitialStep, MaxStep です。入力引数 options を設定しないで得られたデフォルトは、一般には、問題ありません。詳細は、odeset を参照してください。

sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options,p1,p2...) は、付加的なパラメータ p1, p2, ... を、関数 pdefun, icfun, bcfun に渡します。options を設定しない場合、options = [] を使ってください。

注意

• 配列 xmesh と tspan は、pdepe の中では、異なった役割をもっています。

tspan - 関数 pdepe は、時間ステップと公式をダイナミックに選択して、ODE ソルバを使って、時間積分を行います。tspan の要素は、単にユーザが答えたい場所を指定するだけで、コストは、tspan の長さに、ほとんど影響しません。

xmesh - 解への2次の近似は、xmesh に指定されたメッシュ上でなされます。一般的に、解の変化が激しい場合、密なメッシュ点を使用することがベストです。pdepe は、自動的に の中で、メッシュを選択しません。 ユーザは、xmesh に適切な固定メッシュを与える必要があります。コストは、xmesh の長さに強く依存します。 の場合、軸の特異性のために、 近くで、細かいメッシュを使う必要はありません。

• 時間積分は、ode15s を使って、行います。pdepe は、式 1-1 が、楕円方程式を含んでいる場合に生じる微分代数方程式を解くため、また、指定したスパースパターンをもつヤコビアンを取り扱うために、ode15s の機能を使います。
• 離散化の後、楕円方程式は、代数方程式になります。楕円方程式に対応する初期条件の要素が、離散化と矛盾している場合、pdepe は、時間積分の初めに、それらを調整しようとします。このために、最初の段階で戻される解は、他の時点で戻される解と比べて、離散化誤差をもっている可能性があります。メッシュが十分に細かい場合、pdepe は、与えられた条件に近い矛盾のない初期条件を見つけます。pdepe が、矛盾のない初期条件を見つけることが困難であるとのメッセージを表示した場合、メッシュをより細かくして見てください。

放物線方程式に対応する初期条件ベクトルの要素に関しては、調整する必要はありません。

例題

例題 1. つぎの例題は、直接的な定式化、計算、解のプロットを示すものです。

この方程式は、時刻 に対して、区間 で成り立っています。

PDE は、初期条件を満足させます。

それと共に、境界条件も満足させます。

一つの M-ファイルの中に、pdepe で必要とされるすべての関数をサブ関数を使って、設定すると便利です。

function pdex1

m = 0;
x = linspace(0,1,20);
t = linspace(0,2,5);

sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
% 最初の解要素を u として抽出
u = sol(:,:,1);

% サーフェスプロットは、解を調べるために、良い方法となる場合があります。
surf(x,t,u)    
title('Numerical solution computed with 20 mesh points.')
xlabel('Distance x')
ylabel('Time t')

% 解のプロファイルも表示することができます。
figure
plot(x,u(end,:))
title('Solution at t = 2')
xlabel('Distance x')
ylabel('u(x,2)')
% --------------------------------------------------------------
function [c,f,s] = pdex1pde(x,t,u,DuDx)
c = pi^2;
f = DuDx;
s = 0;
% --------------------------------------------------------------
function u0 = pdex1ic(x)
u0 = sin(pi*x);
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
pl = ul;
ql = 0;
pr = pi * exp(-t);
qr = 1;

この例題の中で、PDE、初期条件、境界条件は、サブ関数 pdex1pde, pdex1ic, pdex1bc にコード化されています。

サーフェスプロットは、解の挙動を示します。

つぎのプロットは、t (すなわち、t = 2) の最終値で、解のプロファイルを示しています。

例題 2. つぎの例題は、PDEs システムの解を示すものです。問題は、区間の両端で境界層をもっています。解は、小さな に対して、急激に変化します。

PDEs は、つぎのように表せます。

ここで、

.

です。この方程式は、時刻 で、区間 で成り立ちます。

PDE は、初期条件を満足します。

そして、境界条件も満足します。

pdepe で、使用可能な型として、方程式は、つぎのように変形されます。

に関する偏微係数に対する境界条件は、流量の項目で記述される必要があります。pdepe で、使用可能な型として、左境界条件は、つぎのように変形されます。

また、右境界条件は、つぎのように変形されます。

解は、小さな に対して、急激に変化します。プログラムは、この急激な変化に対応するために、時間に関するステップサイズを選択し、プロット上で、その挙動を調べます。従って、例題では、出力時間を選択する必要があります。[0,1] の両端で、解に境界層が存在しています。それで、例題は、これらの急激な変化を解くために、0 と 1 の近傍にメッシュ点を配置しています。解の挙動を表現するために、適切なメッシュ選択を行うために、数回試みる必要があります。

function pdex4
m = 0;
x = [0 0.005 0.01 0.05 0.1 0.2 0.5 0.7 0.9 0.95 0.99 0.995 1];
t = [0 0.005 0.01 0.05 0.1 0.5 1 1.5 2];

sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t);
u1 = sol(:,:,1);
u2 = sol(:,:,2);

figure
surf(x,t,u1)
title('u1(x,t)')
xlabel('Distance x')
ylabel('Time t')

figure
surf(x,t,u2)
title('u2(x,t)')
xlabel('Distance x')
ylabel('Time t')
% --------------------------------------------------------------
function [c,f,s] = pdex4pde(x,t,u,DuDx)
c = [1; 1];                                  
f = [0.024; 0.17] .* DuDx;                   
y = u(1) - u(2);
F = exp(5.73*y)-exp(-11.47*y);
s = [-F; F];                                                          
% --------------------------------------------------------------
function u0 = pdex4ic(x);
u0 = [1; 0];                                 
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex4bc(xl,ul,xr,ur,t)
pl = [0; ul(2)];                               
ql = [1; 0];                                  
pr = [ur(1)-1; 0];                            
qr = [0; 1];                                  

この例題で、PDEs、初期条件、境界条件は、サブ関数 pdex4pde, pdex4ic, pdex4bc にコード化されています。

サーフェスプロットは、解要素の挙動を示しています。

参考

function handle (function_handle), pdeval, ode15s, odeset, odeget

参考文献

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.

pcolor

pdeval