ドキュメンテーション センター

  • 評価版
  • 製品アップデート

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

偏微分方程式

関数の概要

PDE ソルバー

ここでは MATLAB® の PDE ソルバーを説明します。

PDE 初期値および境界値問題のソルバー

説明

pdepe

1 つの空間変数と時間変数から構成される放物型と楕円型偏微分方程式 (PDE) の初期値および境界値問題を解きます。

PDE 補助関数

PDE 補助関数

説明

pdeval

pdepe の出力を使用して、PDE の数値解を計算します。

初期値問題

pdepe は 1 つの空間変数 x と 時間 t をもつ以下の形式の放物型および楕円型 PDE 系を解きます。

(10-5)

この PDE は t0 ≤ t ≤ tf かつ a ≤ x ≤ b です。区間 [a, b] は有限でなければなりません。m は 0, 1, 2 のいずれかで、それぞれスラブ平面、円筒、球面対称に対応します。m > 0 の場合、a ≥ 0 でなければなりません。

式 10-5 で、f(x,t,u,∂u/∂x) は流量の項、s(x,t,u,∂u/∂x) はソース項です。流量の項は、∂u/∂x に依存することになります。時間に関連する偏導関数の干渉は、対角行列 c(x,t,u,∂u/∂x) との乗算に制限されます。この行列の対角要素は、ゼロまたは正になります。ゼロの要素は楕円方程式に対応し、その他は放物線方程式に対応します。少なくとも 1 つの放物型方程式は存在しなければなりません。放物線方程式に対応する要素 c は、メッシュ点に孤立点 x がある場合、その点での値をゼロにすることができます。物質の境界による不連続 c や s は、メッシュ点が各境界に置かれている場合に使われます。

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

(10-6)

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

(10-7)

q(x, t) は、ゼロまたは非ゼロの要素をもつ対角行列です。この境界条件は、x ∂u/∂x に関連した偏導関数 u ではなく、f の項に関して表現されたものであることに注意してください。また、2 つの係数のうち p のみが u に依存します。

PDE ソルバー

PDE ソルバー

MATLAB の PDE ソルバー pdepe は、1 つの空間変数 x と時間 t で表される放物型または楕円型の連立微分方程式において、初期値と境界値問題を解きます。連立の中には、少なくとも放物線方程式が 1 つ含まれている必要があります。

pdepe ソルバーは、ユーザーが指定したノードをベースにした 2 次精度の空間離散化を使用して、偏微分方程式系を常微分方程式系に変換します。離散化手法は、[9]に記述されています。時間積分は、ode15s を使用して行います。pdepe ソルバーは、式 10-5 が楕円方程式を含む場合に生じる微分代数方程式を解くためと、指定されたスパース パターンをもつヤコビを取り扱うために、ode15s の機能を使います。ode15s は、タイム ステップと公式の両方を動的に変化させます。

離散化後、楕円方程式から代数方程式が生成されます。楕円方程式に対応した初期条件の要素が、離散化と整合性を保っていない場合、pdepe は時間積分を始める前にそれらの調整を試みます。このため最初の時刻に返される解は、他の時間のものと同程度の離散化誤差をもっています。メッシュが十分細かい場合、pdepe は与えられた初期条件に近い整合性のある初期条件を検出することができます。pdepe が整合性のある初期条件を見つけることが困難であるというメッセージを表示される場合は、メッシュをより細かく分割してください。放物線方程式に対応する初期条件ベクトルの要素の調整は必要ありません。

PDE ソルバー構文

ソルバーの基本的な構文は、次の型をしています。

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

入力引数は以下になります。

m

問題の対称性を指定します。m0 の場合にはスラブ平面対称、1 の場合には円筒対称、2 の場合は球面対称です。 = 式 10-5m に対応します。

pdefun

PDE の要素を定義する関数です。この関数は、式 10-5 で項 c、f および s を計算し、次の型をもちます。

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

ここで、xt はスカラー、ududx はベクトルで、解 u と x に関するその偏導関数を近似します。cf および s は列ベクトルです。c は行列 c の対角要素を格納します。

icfun

初期条件を計算する関数です。次の型をもちます。

u = icfun(x)

引数 x と共に呼び出すと、icfun は列ベクトル ux における解要素の初期値を計算し、返します。

bcfun

境界条件の項 p および q を計算する関数です。次の型をもちます。

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

ここで ul は左境界 xl = a の近似解であり、ur は右境界 xr = b の近似解です。plqlxl で計算された p に対応する列ベクトルと q の対角要素です。同様に、prqrxr に対応します。m > 0 かつ a = 0 のとき、x = 0 の近傍での解の制約条件は、流速 f が a = 0 でゼロになる必要があります。pdepe は、この境界条件を自動的に適用し、plql の中に戻る値を無視します。

xmesh

tspan のすべての値への数値解が要求される点を設定するベクトル [x0, x1, ..., xn] です。x0xn は、それぞれ a と b に対応します。

解の 2 次近似は xmesh で指定したメッシュ上で行われます。一般に、解の変化が激しい部分では、密なメッシュ点を使うことを推奨します。pdepe は x でメッシュの選択を自動的には選択 "しません"。ユーザーは固定された適切なメッシュを xmesh に与える必要があります。コストは、xmesh の長さに強く依存します。m > 0 の場合、座標特異点という理由で x = 0 の近くで細かいメッシュを使用する必要はありません。

xmesh の要素は、x0 < x1 < ... < xn を満たさなければなりません。xmesh の長さは 3以上 ( ≥ 3) でなければなりません。

tspan

xmesh のすべての値への解が要求される点を設定するベクトル [t0, t1, ..., tf] です。t0tf は、それぞれ t0 と tf に対応します。

pdepe は ODE ソルバーを使って、時間積分を行うときに、動的にタイム ステップと式を選択します。tspan で設定した点での解は、積分の定式を連続的に拡張する手法を使用して得られます。tspan の要素は、ユーザーが参照する位置を指定するのみのもので、tspan の長さによっては若干コストに差が出ます。

tspan の要素は、t0 < t1 < ... < tf を満たさなければなりません。tspan の長さは、3 以上 (≥ 3) である必要があります。

出力引数 sol は、以下の条件を満たす 3 次元配列になります。

  • sol(:,:,k) が u の解の k 番目の要素を近似する。

  • sol(i,:,k) が、時刻 tspan(i)、メッシュ点 xmesh(:) で、解の k 番目の要素を近似する。

  • sol(i,j,k) が、時刻 tspan(i)、メッシュ点 xmesh(j) で、解の k 番目の要素を近似する。

PDE ソルバー オプション

高度な応用として、入力引数に PDE ソルバーに渡すオプションと付加的なパラメーターを設定できます。

options

既定の積分プロパティを変更するオプション パラメーターの構造体。これは、7 番目の入力引数として設定します。

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

詳細は、「積分器のオプション」を参照してください。

積分器のオプション

MATLAB の PDE ソルバーの既定積分プロパティは、一般的な問題を扱えるように選択されています。場合によってはこれらの既定値を変更することで、ソルバーのパフォーマンスを向上させることができます。これを行うには、options 構造体に含まれる単一または複数のプロパティ値をソルバー pdepe に渡します。

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

関数 odeset を使って、options 構造体を作成します。pdepe で使用できるオプションは、基礎となる ODE ソルバーがもっているオプションのうち、次の表に示すもののみになります。一般的には入力引数 options を使わない既定設定で十分です。構造体の作成やプロパティの記述方法は、「積分器のオプション」を参照してください。

PDE プロパティ

プロパティ カテゴリ

プロパティ名

誤差制御

RelTol, AbsTol, NormControl

ステップ サイズ

InitialStep, MaxStep

単一の PDE

方程式を解く-  次の例は、以下に示すシンプルな PDE の簡単な定式化、解の計算およびプロットを説明します。

この方程式は、時間 t ≥ 0 に対して、区間 0 ≤ x ≤ 1 で成り立ちます。t = 0 で、解は以下の初期条件を満たします。

x = 0 および x = 1 で、解は以下の境界条件を満たします。

    メモ:   ファイル pdex1.mpdex1.m には、この例のコードがすべて含まれています。このファイルでは、必要な関数がすべてローカル関数としてコード化されています。エディターにコードを表示するには、コマンド ラインで「edit pdex1」と入力します。これを実行する場合は、コマンド ラインで「pdex1」と入力します。詳細は、「PDE ソルバー構文」を参照してください。

  1. PDE を書き換える。次の型に PDE を書き換えます。

    これは、式 10-5 で示される型で pdepe で使用できます。詳細は、「初期値問題」を参照してください。この例では、結果の方程式は以下になります。

    パラメーター m = 0 と以下の項が付随します。

  2. PDE をコード化する。上に示した型に PDE を書き直し (式 10-5)、各項を設定すると、pdepe が使用可能な関数として、PDE をコード化できます。関数は次の形式になります。

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

    ここで、cf および s は、c、f および s 項に対応します。次のコードは、例の c, f, s を計算します。

    function [c,f,s] = pdex1pde(x,t,u,DuDx)
    c = pi^2;
    f = DuDx;
    s = 0;
    
  3. 初期条件関数をコード化します。初期条件は、以下の形式の関数にコード化する必要があります。

    u = icfun(x)
    

    次のコードは、MATLAB の関数 pdex1ic の初期条件を表しています。

    function u0 = pdex1ic(x)
    u0 = sin(pi*x);
    
  4. 境界条件関数をコード化します。境界条件は、以下の形式の関数にコード化する必要があります。

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

    式 10-7 と同じ型で記述された境界条件は、次のようになります。

    および

    次のコードは、境界条件の要素 p(x,t,u) および q(x,t) を関数 pdex1bc で計算します。

    function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
    pl = ul;
    ql = 0;
    pr = pi * exp(-t);
    qr = 1;
    

    関数 pdex1bc で、plql は左の境界条件 (x = 0) に対応し、prqr は右の境界条件 (x = 1) に対応します。

  5. 解のメッシュ点を選択する。MATLAB PDE ソルバーを使用する前に、メッシュ点 (t,x) を、pdepe が解を計算する点として指定する必要があります。これらは、ベクトル tx に設定します。

    ベクトル tx は、ソルバーの中で異なる役割をもちます (「PDE ソルバー」を参照)。特に計算のコストと精度は、ベクトル x の長さに強く依存します。しかし、ベクトル t の値にはあまり影響を受けません。

    この例は、空間区間 [0,1] を 20 等分のメッシュに分割し、時間間隔 [0,2] を 5 等分した t に対して解を計算します。

    x = linspace(0,1,20);
    t = linspace(0,2,5);
    
  6. PDE ソルバーを適用する。この例は m = 0 として、関数 pdex1pdepdex1ic および pdex1bcpdepe が解を計算する位置である xt で定義されるメッシュによって pdepe を呼び出します。関数 pdepe は 3 次元配列 sol に数値解を返します。ここで sol(i,j,k) は解の k 番目の要素 (ukt(i)x(j) で計算) を近似します。

    m = 0;
    sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
    

    この例は、@ を使って pdex1pdepdex1icpdex1bc を、pdepe への関数ハンドルとして渡します。

      メモ:   関数ハンドルの詳細は、MATLAB のプログラミングの基礎ドキュメンテーションの function_handle (@)、func2strstr2func のリファレンス ページ、および「@」を参照してください。

  7. 結果を表示する。例の最後に、結果を表示します。

    1. 最初の解の成分を抽出し、表示します。この例では、解 u の成分は 1 つのみですが、説明のため 3 次元配列からそれを抽出します。表面プロットは解の動作を示します。

          u = sol(:,:,1);
          
          surf(x,t,u)    
          title('Numerical solution computed with 20 mesh points')
          xlabel('Distance x')
          ylabel('Time t')
      

    2. t の最終値 tf で、解のプロファイルを表示します。この例では、tf  = t = 2 になります。

          figure
          plot(x,u(end,:))
          title('Solution at t = 2')
          xlabel('Distance x')
          ylabel('u(x,2)')
      

解の計算-  上記のように解を得てプロットした後、特定の t 値での解の計算結果や、特定の点 x での解の時間変化に興味がある場合があります。(ステップ 7 で抽出した解の) k 番目の列 u(:,k) は、x(k) での解の時系列履歴を含んでいます。j 番目の行 u(j,:) は、t(j) での解のプロファイルをもちます。

ベクトル xu(j,:)、および補助関数 pdeval を使用して、解 u とその微分 ∂u/∂x を任意の点のセット xout で求められます。

[uout,DuoutDx] = pdeval(m,x,u(j,:),xout)

例「pdex3pdex3」は、xout = 0 における解の導関数を関数 pdeval を使用して計算します。詳細は関数 pdeval を参照してください。

PDE のシステム

この例では、偏微分方程式系の解法を示します。問題は、電気力学の分野から採用しています。この問題は、積分区間の両端で境界層をもち、解は次の微小な値 t で急激に変化します。

PDE は以下になります。

ここで、F(y) = exp(5.73y) – exp(–11.46y) です。これらの方程式は、時間 t ≥ 0 に対して、区間 0 ≤ x ≤ 1 で成り立ちます。

解 u は以下の初期条件を満たします。

および以下の境界条件を満たします。

    メモ:   ファイル pdex4.mpdex4.m には、この例のコードがすべて含まれています。このファイルでは、必要な関数がすべてローカル関数としてコード化されています。エディターにコードを表示するには、コマンド ラインで「edit pdex4」と入力します。これを実行する場合は、コマンド ラインで「pdex4」と入力します。

  1. PDE を書き換える。pdepe で想定される型において、方程式は次のようになります。

    u の偏導関数の境界条件は、流速の項で記述される必要があります。pdepe で想定される型において、左境界条件は次のようになります。

    右境界条件は次のようになります。

  2. PDE をコード化する。上記で示した型への PDE を書き換えると、pdepe で使用できる関数としてコード化できます。関数は次の形式になります。

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

    ここで、cf および s は、式 10-5 の c、f および s 項に対応します。

    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];  
    
  3. 境界条件関数をコード化する。初期条件関数は、次の型になります。

    u = icfun(x)
    

    次のコードは、関数 pdex4ic として初期条件を表しています。

    function u0 = pdex4ic(x);
    u0 = [1; 0]; 
    
  4. 境界条件関数をコード化する。境界条件関数は次の型になる必要があります。

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

    次のコードは、境界条件の要素 p(x,t,u) および q(x,t) (式 10-7) を関数 pdex4bc で計算するものです。

    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];     
    
  5. 解のメッシュ点を選択する。解は、値の小さい t で急激に変化します。プログラムは、この急激な変化に追従するようにタイム ステップ サイズを選択します。ただし、プロットでもこの動作を見るには、出力時間をそれに応じて選択しなければなりません。この解では [0,1] の両端において境界層があるので、メッシュ点はこの急激な変化を解像できるように配置する必要があります。解の動作を表すメッシュを選択するには、何回かテストをする必要があります。

    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];
    
  6. PDE ソルバーを適用する。この例は m = 0 をもつ pdepe、関数 pdex4pde、関数 pdex4ic、関数 pdex4bc を呼び出し、pdepe が解を計算する位置を x と t で定義されるメッシュによって指定します。関数 pdepe は 3 次元配列 sol に数値解を返します。ここで sol(i,j,k) は解の k 番目の要素 (μkt(i)x(j) で計算) を近似します。

    m = 0;
    sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t);
    
  7. 結果を表示する。表面プロットは解の動作を示します。

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

その他の例

以下の例を使用できます。次のように入力します。

edit examplename

を入力してコードを表示し、

examplename

で例を実行します。

例の名前

説明

pdex1

シンプルな PDE を使用して、簡単な定式化、計算、解のプロット手法を示します。

pdex2

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

pdex3

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

pdex4

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

pdex5

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

この情報は役に立ちましたか?