このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。
単一の PDE の求解
この例では、単一の PDE の定式化、計算、および解のプロット方法を説明します。
次の偏微分方程式を考えます。
この方程式は、時間 について区間 で定義されます。 で、解は次の初期条件を満たします。
また、 および で、解は以下の境界条件を満たします。
この方程式を MATLAB® で解くには、方程式、初期条件、および境界条件をコード化し、適切な解のメッシュを選択してから、ソルバー pdepe
を呼び出す必要があります。(ここで行ったように) 必要な関数をファイルの最後にローカル関数として含めることも、あるいは個別の名前付きファイルとして MATLAB パスのディレクトリに保存することもできます。
方程式のコード化
方程式は、コード化する前に pdepe
ソルバーで想定される形式に書き換える必要があります。pdepe
で想定される標準形式は次のとおりです。
.
この形式で記述すると、PDE は次のようになります。
.
適切な形式の方程式を使用すると、関連する項を読み取ることができます。
これで、方程式をコード化する関数を作成できるようになりました。この関数にはシグネチャ [c,f,s] = pdex1pde(x,t,u,dudx)
がなければなりません。
x
は独立空間変数。t
は独立時間変数。u
は、x
およびt
について微分される従属変数。dudx
は空間偏導関数 。出力
c
、f
、およびs
は、pdepe
で想定される標準 PDE 方程式形式の係数に対応。これらの係数は、入力変数x
、t
、u
、およびdudx
によってコード化されます。
結果として、この例の方程式は次の関数で表すことができます。
function [c,f,s] = pdex1pde(x,t,u,dudx) c = pi^2; f = dudx; s = 0; end
(メモ: 関数はすべて例の最後にローカル関数として含まれます。)
初期条件のコード化
次に、初期条件を返す関数を記述します。初期条件は、最初の時間値 tspan(1)
で適用されます。この関数にはシグネチャ u0 = pdex1ic(x)
がなければなりません。
対応する関数は次のようになります。
function u0 = pdex1ic(x) u0 = sin(pi*x); end
境界条件のコード化
次に、境界条件を評価する関数を記述します。区間 で提示される問題について、境界条件はすべての 、および または のいずれかに適用されます。ソルバーで想定される境界条件の標準形式は次のとおりです。
この標準形式で境界条件を書き直し、係数値を読み取ります。
の場合、方程式は次のようになります。
係数は次のようになります。
の場合、方程式は次のようになります。
係数は次のようになります。
境界条件関数は によって表現され、この項は既にメインの PDE 関数で定義されているため、方程式のこの部分を境界条件関数内で指定する必要はありません。 および の値のみを各境界で指定する必要があります。
境界関数は関数シグネチャ [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
を使用しなければなりません。
入力
xl
およびul
は左境界の および に対応。入力
xr
およびur
は右境界の および に対応。t
は独立時間変数。出力
pl
およびql
は左境界 (この問題では ) の および に対応。出力
pr
およびqr
は右境界 (この問題では ) の および に対応。
この例の境界条件は次の関数で表されます。
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t) pl = ul; ql = 0; pr = pi * exp(-t); qr = 1; end
解のメッシュの選択
方程式を解く前に、メッシュ点 を、pdepe
が解を評価する点として指定する必要があります。これらは、ベクトル t
と x
に設定します。ベクトル t
と x
は、ソルバーの中で異なる役割をもちます。特に求解のコストと精度は、ベクトル x
の長さに強く依存します。しかし、ベクトル t
の値にはあまり影響を受けません。
この問題では、空間区間 [0,1] 内で等間隔に配置された 20 個の点をもつメッシュと、時間区間 [0,2] からの t
の 5 つの値を使用します。
x = linspace(0,1,20); t = linspace(0,2,5);
方程式の求解
最後に、対称性 m
、PDE 方程式、初期条件、境界条件、および x
と t
のメッシュを使用して方程式を解きます。
m = 0; sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
pdepe
は 3 次元配列 sol
に解を返します。ここで、sol(i,j,k)
は t(i)
と x(j)
で評価した解の k
番目の要素 を近似します。u0
が各解要素の初期条件を指定するため、sol
のサイズは length(t)
x length(x)
x length(u0)
となります。この問題では、u
の要素は 1 つのみであるため、sol
は 5 行 20 列の行列になりますが、通常はコマンド u = sol(:,:,k)
を使用して k
番目の解要素を抽出できます。
最初の解要素を sol
から抽出します。
u = sol(:,:,1);
解のプロット
解の表面プロットを作成します。
surf(x,t,u) title('Numerical solution computed with 20 mesh points') xlabel('Distance x') ylabel('Time t')
この問題の初期条件と境界条件は、次で与えられる解析解となるように選択されました。
同じメッシュ点で解析解をプロットします。
surf(x,t,exp(-t)'*sin(pi*x)) title('True solution plotted with 20 mesh points') xlabel('Distance x') ylabel('Time t')
次に、数値解と解析解を の最終値である において比較します。この例では、 となります。
plot(x,u(end,:),'o',x,exp(-t(end))*sin(pi*x)) title('Solution at t = 2') legend('Numerical, 20 mesh points','Analytical','Location','South') xlabel('Distance x') ylabel('u(x,2)')
ローカル関数
ここでは、PDE ソルバー pdepe が解を計算するために呼び出すローカル補助関数を紹介しています。あるいは、これらの関数を独自のファイルとして MATLAB パスのディレクトリに保存することもできます。
function [c,f,s] = pdex1pde(x,t,u,dudx) % Equation to solve c = pi^2; f = dudx; s = 0; end %---------------------------------------------- function u0 = pdex1ic(x) % Initial conditions u0 = sin(pi*x); end %---------------------------------------------- function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t) % Boundary conditions pl = ul; ql = 0; pr = pi * exp(-t); qr = 1; end %----------------------------------------------