不連続性をもつ PDE の求解

この例では、物質と境界を接する PDE を解く方法を説明します。この問題では $x = 0.5$ において物質の境界により不連続性が作成され、初期条件には右の境界 $x = 1$ に不連続性があります。

区分的な PDE を考えます。

${\begin{matrix} \frac{\partial u}{\partial t} = x^{- 2} \frac{\partial}{\partial x} (x^{2} 5 \frac{\partial u}{\partial x}) - 1000 e^{u} & (0 \leq x \leq 0.5) \\ \frac{\partial u}{\partial t} = x^{- 2} \frac{\partial}{\partial x} (x^{2} \frac{\partial u}{\partial x}) - e^{u} & (0.5 \leq x \leq 1) \end{matrix}$

初期条件は次のとおりです。

$\begin{array}{l} u (x, 0) = 0 (0 \leq x < 1), \\ u (1, 0) = 1 (x = 1) . \end{array}$

境界条件は次のとおりです。

$\begin{array}{l} \frac{\partial u}{\partial x} = 0 (x = 0), \\ u (1, t) = 1 (x = 1) . \end{array}$

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

方程式のコード化

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

$c (x, t, u, \frac{\partial u}{\partial x}) \frac{\partial u}{\partial t} = x^{- m} \frac{\partial}{\partial x} (x^{m} f (x, t, u, \frac{\partial u}{\partial x})) + s (x, t, u, \frac{\partial u}{\partial x})$ .

この場合、PDE は適切な形式であるため、係数の値を読み取ることができます。

流束項 $f (x, t, u, \frac{\partial u}{\partial x})$ とソース項 $s (x . t, u, \frac{\partial u}{\partial x})$ の値は、 $x$ の値に応じて変化します。係数は次のようになります。

$m = 2$

$c (x, t, u, \frac{\partial u}{\partial x}) = 1$

${\begin{matrix} f (x, t, u, \frac{\partial u}{\partial x}) = 5 \frac{\partial u}{\partial x} & (0 \leq x \leq 0.5) \\ f (x, t, u, \frac{\partial u}{\partial x}) = \frac{\partial u}{\partial x} & (0.5 \leq x \leq 1) \end{matrix}$

${\begin{matrix} s (x, t, u, \frac{\partial u}{\partial x}) = - 1000 e^{u} & (0 \leq x \leq 0.5) \\ s (x, t, u, \frac{\partial u}{\partial x}) = - e^{u} & (0.5 \leq x \leq 1) \end{matrix}$

これで、方程式をコード化する関数を作成できるようになりました。この関数にはシグネチャ [c,f,s] = pdex2pde(x,t,u,dudx) がなければなりません。

x は独立空間変数。
t は独立時間変数。
u は、x および t について微分される従属変数。
dudx は空間偏導関数 $\partial u / \partial x$ 。
出力 c、f、および s は、pdepe で想定される標準 PDE 方程式形式の係数に対応。これらの係数は、入力変数 x、t、u、および dudx によってコード化されます。

結果として、この例の方程式は次の関数で表すことができます。

function [c,f,s] = pdex2pde(x,t,u,dudx)
c = 1;
if x <= 0.5
    f = 5*dudx;
    s = -1000*exp(u);
else
    f = dudx;
    s = -exp(u);
end
end

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

初期条件のコード化

次に、初期条件を返す関数を記述します。初期条件は最初の時間値において適用され、 $u (x, t_{0})$ の値を任意の $x$ の値に対して提供します。関数シグネチャ u0 = pdex2ic(x) を使用して関数を記述します。

初期条件は次のとおりです。

$\begin{array}{l} u (x, 0) = 0 (0 \leq x < 1), \\ u (1, 0) = 1 (x = 1) . \end{array}$

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

function u0 = pdex2ic(x)
if x < 1
    u0 = 0;
else
    u0 = 1;
end
end

境界条件のコード化

次に、境界条件を評価する関数を記述します。区間 $a \leq x \leq b$ で提示される問題について、境界条件はすべての t、および $x = a$ または $x = b$ のいずれかに適用されます。ソルバーで想定される境界条件の標準形式は次のとおりです。

$p (x, t, u) + q (x, t) f (x, t, u, \frac{\partial u}{\partial x}) = 0 .$

この例には球面対称性 ( $m = 2$ ) があるため、pdepe ソルバーは左境界条件を自動的に強制して解を原点に拘束し、境界関数内で左境界に指定された条件はすべて無視します。したがって、左境界条件には $p_{L} = q_{L} = 0$ を指定できます。右境界条件については、境界条件を標準形式で書き直し、 $p_{R}$ と $q_{R}$ の係数値を読み取ることができます。

$x = 1$ に対し、方程式は $u (1, t) = 1 \to (u - 1) + 0 \cdot \frac{\partial u}{\partial x} = 0$ になります。係数は次のようになります。

$p_{R} (1, t, u) = u - 1$
$q_{R} (1, t) = 0$

境界関数は関数シグネチャ [pl,ql,pr,qr] = pdex2bc(xl,ul,xr,ur,t) を使用しなければなりません。

入力 xl および ul は左境界の u および x に対応。
入力 xr および ur は右境界の u および x に対応。
t は独立時間変数。
出力 pl および ql は左境界 (この問題では $x = 0$ ) の $p_{L} (x, t, u)$ および $q_{L} (x, t)$ に対応。
出力 pr および qr は右境界 (この問題では $x = 1$ ) の $p_{R} (x, t, u)$ および $q_{R} (x, t)$ に対応。

この例の境界条件は次の関数で表されます。

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

解のメッシュの選択

空間メッシュは、不連続な境界を考慮するために $x = 0.5$ 近傍の値をいくつか含んでいなければならず、また、 $x = 1$ 近傍の点を含んでいなければなりません。これは、その点での初期値 ( $u (1, 0) = 1$ ) と境界値 ( $u (1, t) = 0$ ) が一致しないためです。小さい $t$ に対して解は急速に変化するため、この急激な変化を解決できるタイムステップを使用します。

x = [0 0.1 0.2 0.3 0.4 0.45 0.475 0.5 0.525 0.55 0.6 0.7 0.8 0.9 0.95 0.975 0.99 1];
t = [0 0.001 0.005 0.01 0.05 0.1 0.5 1];

方程式の求解

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

m = 2;
sol = pdepe(m,@pdex2pde,@pdex2ic,@pdex2bc,x,t);

pdepe は 3 次元配列 sol に解を返します。ここで、sol(i,j,k) は t(i) と x(j) で評価した解の k 番目の要素 $u_{k}$ を近似します。u0 が各解要素の初期条件を指定するため、sol のサイズは length(t) x length(x) x length(u0) となります。この問題では、u の要素は 1 つのみであるため、sol は 8 行 18 列の行列になりますが、通常はコマンド u = sol(:,:,k) を使用して k 番目の解要素を抽出できます。

最初の解要素を sol から抽出します。

u = sol(:,:,1);

解のプロット

$x$ と $t$ の選択されたメッシュ点でプロットされた解 $u$ の表面プロットを作成します。 $m = 2$ であるため、問題は球面対称性をもつ球面幾何学で提示されます。したがって、解は半径 $x$ の方向にのみ変化します。

surf(x,t,u)
title('Numerical solution with nonuniform mesh')
xlabel('Distance x')
ylabel('Time t')
zlabel('Solution u')

Figure contains an axes object. The axes object with title Numerical solution with nonuniform mesh, xlabel Distance x, ylabel Time t contains an object of type surface.

次に、 $x$ と $u$ のみをプロットして、表面プロットの等高線の側面図を取得します。 $x = 0.5$ にラインを追加して、物質の境界の影響を強調表示します。

plot(x,u,x,u,'*')
line([0.5 0.5], [-3 1], 'Color', 'k')
xlabel('Distance x')
ylabel('Solution u')
title('Solution profiles at several times')

Figure contains an axes object. The axes object with title Solution profiles at several times, xlabel Distance x, ylabel Solution u contains 17 objects of type line. One or more of the lines displays its values using only markers

ローカル関数

ここでは、PDE ソルバー pdepe が解を計算するために呼び出すローカル補助関数を紹介しています。あるいは、これらの関数を独自のファイルとして MATLAB パスのディレクトリに保存することもできます。

function [c,f,s] = pdex2pde(x,t,u,dudx) % Equation to solve
c = 1;
if x <= 0.5
    f = 5*dudx;
    s = -1000*exp(u);
else
    f = dudx;
    s = -exp(u);
end
end
%----------------------------------------------
function u0 = pdex2ic(x) %Initial conditions
if x < 1
    u0 = 0;
else
    u0 = 1;
end
end
%----------------------------------------------
function [pl,ql,pr,qr] = pdex2bc(xl,ul,xr,ur,t) % Boundary conditions
pl = 0;
ql = 0;
pr = ur - 1;
qr = 0;
end
%----------------------------------------------

参考

pdepe