pdeCoefficients

偏微分方程式の係数の抽出

R2021a 以降

構文

coeffs = pdeCoefficients(pdeeq,u)

symCoeffs = pdeCoefficients(pdeeq,u,'Symbolic',true)

説明

coeffs = pdeCoefficients(pdeeq,u) は、偏微分方程式 (PDE) の係数を倍精度数と関数ハンドルの構造体として抽出します。この構造体は、Partial Differential Equation Toolbox™ で関数 specifyCoefficients の入力として使用できます。

pdeeq は、スカラーの PDE またはシンボリック形式の PDE 系で、u の関数です。関数 pdeCoefficients は、pdeeq を

$$m\frac{{\partial }^{2}u}{\partial t^2}+d\frac{\partial u}{\partial t}-\nabla ·\left(c\otimes \nabla u\right)+au=f$$

の形式の連立方程式に変換し、係数 m、d、c、a、および f を含む構造体 coeffs を返します。詳細については、Equations You Can Solve Using Partial Differential Equation Toolbox (Partial Differential Equation Toolbox)を参照してください。

pdeCoefficients は、PDE を上記の発散形式に変換できない場合、警告メッセージを表示し、残りのすべての勾配を係数 f に書き込みます。Partial Differential Equation Toolbox は、このタイプの PDE を解くことができません。

symCoeffs = pdeCoefficients(pdeeq,u,'Symbolic',true) は、PDE の係数をシンボリック式の構造体として返します。

例

シンボリック PDE からの係数の抽出

変数 u(x,y) のラプラシアンを表すシンボリック PDE を作成します。

syms u(x,y)
pdeeq = laplacian(u,[x y]) == -3

pdeeq(x, y) = 
$$\frac{{\partial }^2}{\partial x^2}\mathrm{ }u\left(x,y\right)+\frac{{\partial }^2}{\partial y^2}\mathrm{ }u\left(x,y\right)=-3$$

PDE の係数を抽出します。

coeffs = pdeCoefficients(pdeeq,u)

coeffs = struct with fields:
    a: 0
    c: [4×1 double]
    m: 0
    d: 0
    f: -3

pdeCoefficients は、シンボリック PDE を以下の形式のスカラー PDE 方程式に変換します。

$$m\frac{{\partial }^{2}u}{\partial t^2}+d\frac{\partial u}{\partial t}-\nabla \cdot (c\nabla u)+au=f$$

また、係数 m、d、c、a、および f を抽出し、構造体 coeffs に格納します。2 次元の連立 $$N$$ 次方程式の場合、c は 4$$N^2$$ 個の要素をもつテンソルとなります。詳細については、c Coefficient for specifyCoefficients (Partial Differential Equation Toolbox)を参照してください。この例では $$N=1$$ であるため、係数 c は、c$${}_{xx}$$、c$${}_{xy}$$、c$${}_{yx}$$、および c$${}_{yy}$$ を表す 4 行 1 列の列ベクトルとなります。

coeffs.c

ans = 4×1

    -1
     0
     0
    -1

斉次ラプラス方程式の求解

単位円で囲まれた領域における 2 次元の斉次ラプラス方程式を解きます。関数 pdeCoefficients を使用して、ラプラス方程式の係数を抽出します。

PDE モデルを作成し、このモデルに幾何形状を含めます。

model = createpde;
geometryFromEdges(model,@circleg);

幾何形状をプロットし、境界条件の定義に使用する辺のラベルを表示します。

figure
pdegplot(model,'EdgeLabels','on')
axis equal

ラプラス方程式を表すシンボリック式 pdeeq を作成します。

syms u(x,y)
pdeeq = laplacian(u,[x y])

pdeeq(x, y) = 
$$\frac{{\partial }^2}{\partial x^2}\mathrm{ }u\left(x,y\right)+\frac{{\partial }^2}{\partial y^2}\mathrm{ }u\left(x,y\right)$$

ラプラス方程式の係数を抽出します。

coeffs = pdeCoefficients(pdeeq,u)

coeffs = struct with fields:
    a: 0
    c: [4×1 double]
    m: 0
    d: 0
    f: 0

PDE モデルの係数を指定します。

specifyCoefficients(model,'m',coeffs.m,'d',coeffs.d, ...
    'c',coeffs.c,'a',coeffs.a,'f',coeffs.f);

円を形成する 4 辺すべてにディリクレ境界条件 $$u(x,y)=x^4+y^4$$ を適用します。

applyBoundaryCondition(model,'dirichlet','Edge',1:4,'u',@(region,state) region.x.^4 + region.y.^4);

この幾何形状に対する既定のメッシュを生成します。

generateMesh(model);

PDE を解いて解をプロットします。

results = solvepde(model);
pdeplot(model,'XYData',results.NodalSolution)

非定数係数 f をもつポアソン方程式の求解

単位円で囲まれた領域における 2 次元のポアソン方程式を解きます。この PDE の発散形式は、非定数係数 f をもちます。この PDE を解くには、pdeCoefficients を使用して係数を抽出し、specifyCoefficients を使用して PDE モデルの係数を指定します。

PDE モデルを作成し、このモデルに幾何形状を含めます。

model = createpde;
geometryFromEdges(model,@circleg);

ポアソン方程式を表すシンボリック式 pdeeq を作成します。

syms u(x,y)
pdeeq = diff(u,x,x) + diff(u,y,y) - 1/(x^2 + y^2)

pdeeq(x, y) = 
$$\frac{{\partial }^2}{\partial x^2}\mathrm{ }u\left(x,y\right)-\frac{1}{x^2+y^2}+\frac{{\partial }^2}{\partial y^2}\mathrm{ }u\left(x,y\right)$$

方程式の係数を抽出します。

coeffs = pdeCoefficients(pdeeq,u)

coeffs = struct with fields:
    a: 0
    c: [4×1 double]
    m: 0
    d: 0
    f: @makeCoefficient/coefficientFunction

係数 f は定数ではありません。これは @makeCoefficient/coefficientFunction として表示されます。これは、この係数が内部関数のワークスペースから返されたことを表します。関数ハンドルの背後にある式を表示するには、coeffs.f('show') を使用します。

coeffs.f('show')

ans = 
$$\frac{1}{x^2+y^2}$$

PDE モデルの係数を指定します。

specifyCoefficients(model,'m',coeffs.m,'d',coeffs.d, ...
    'c',coeffs.c,'a',coeffs.a,'f',coeffs.f);

円を形成する 4 辺すべてにディリクレ境界条件 $$u(x,y)=0$$ を適用します。

applyBoundaryCondition(model,'dirichlet','Edge',1:4,'u',0);

この幾何形状に対する既定のメッシュを生成します。

generateMesh(model);

PDE を解きます。関数 pdeplot のオプション 'XYData' を使用して、ノード解をプロットします。

results = solvepde(model);
pdeplot(model,'XYData',results.NodalSolution)

オプション 'FlowData' を使用して、ノード位置における解の勾配をプロットします。

pdeplot(model,'FlowData',[results.XGradients results.YGradients])

非発散形式の PDE

非発散形式で PDE を構築します。

syms u(x,y)
pdeeq = diff(u,x,x) + cos(x+y)/4*diff(u,x,y) + 1/2*diff(u,y,y)

pdeeq(x, y) = 
$$\frac{{\partial }^2}{\partial x^2}\mathrm{ }u\left(x,y\right)+\frac{\cos \left(x+y\right)\hspace{0.17em}\frac{\partial }{\partial y}\mathrm{ }\frac{\partial }{\partial x}\mathrm{ }u\left(x,y\right)}{4}+\frac{\frac{{\partial }^2}{\partial y^2}\mathrm{ }u\left(x,y\right)}{2}$$

係数を抽出します。pdeCoefficients は、PDE を以下の発散形式に変換できない場合があります。

$$m\frac{{\partial }^{2}u}{\partial t^2}+d\frac{\partial u}{\partial t}-\nabla \cdot (c\otimes \nabla u)+au=f$$,

この場合、警告メッセージを表示しますが、結果も出力します。

coeffs = pdeCoefficients(pdeeq,u)

Warning: After extracting m, d, and c, some gradients remain. Writing all remaining terms to f.

coeffs = struct with fields:
    a: 0
    c: @makeCoefficient/coefficientFunction
    m: 0
    d: 0
    f: @makeCoefficient/coefficientFunction

抽出した係数 c および f の関数ハンドルを表示するには、'show' オプションを使用します。PDE に含まれる残りの勾配は、すべて係数 f に書き込まれます。

coeffs.c('show')

ans = 
$$\left(\begin{array}{c}-1\\ \frac{1}{8}-\frac{\cos \left(x+y\right)}{4}\\ \frac{1}{8}-\frac{\cos \left(x+y\right)}{4}\\ -\frac{1}{2}\end{array}\right)$$

coeffs.f('show')

ans = 
$$\frac{\cos \left(x+y\right)\hspace{0.17em}\frac{\partial }{\partial y}\mathrm{ }\frac{\partial }{\partial x}\mathrm{ }u\left(x,y\right)}{4}-\frac{\sin \left(x+y\right)\hspace{0.17em}\frac{\partial }{\partial y}\mathrm{ }u\left(x,y\right)}{4}-\frac{\sin \left(x+y\right)\hspace{0.17em}\frac{\partial }{\partial x}\mathrm{ }u\left(x,y\right)}{4}-\frac{\frac{\partial }{\partial y}\mathrm{ }\frac{\partial }{\partial x}\mathrm{ }u\left(x,y\right)}{4}$$

この PDE には PDE Toolbox が必要とする発散形式が含まれていないため、PDE Toolbox はこの PDE を解くことができません。

時間依存波動方程式の求解

単位円で囲まれた領域における時間依存波動方程式を解きます。この波動方程式は、時間 $$t$$ と座標 $$x$$ および $$y$$ に依存するスカラー関数 $$u(t,x,y)$$ をもつ PDE です。

PDE モデルを作成し、このモデルに幾何形状を含めます。

model = createpde(1);
geometryFromEdges(model,@circleg);

波動方程式を表すシンボリック PDE を作成します。

syms u(t,x,y)
pdeeq = diff(u,t,t) - laplacian(u,[x y])

pdeeq(t, x, y) = 
$$\frac{{\partial }^2}{\partial t^2}\mathrm{ }u\left(t,x,y\right)-\frac{{\partial }^2}{\partial x^2}\mathrm{ }u\left(t,x,y\right)-\frac{{\partial }^2}{\partial y^2}\mathrm{ }u\left(t,x,y\right)$$

PDE の係数を抽出します。

coeffs = pdeCoefficients(pdeeq,u)

coeffs = struct with fields:
    a: 0
    c: [4×1 double]
    m: 1
    d: 0
    f: 0

PDE モデルの係数を指定します。

specifyCoefficients(model,'m',coeffs.m,'d',coeffs.d, ...
    'c',coeffs.c,'a',coeffs.a,'f',coeffs.f);

幾何形状全体に関する時間依存問題の初期条件を設定します。

setInitialConditions(model,0,1);

円を形成する 4 辺すべてにディリクレ境界条件 $$u(x,y)=0$$ を適用します。

applyBoundaryCondition(model,'dirichlet','Edge',1:4,'u',0);

この幾何形状に対する既定のメッシュを生成します。

generateMesh(model);

時間範囲 0 秒 ～ 50 秒、時間間隔 2 秒の条件で、時間依存 PDE の解を求めます。

results = solvepde(model,linspace(0,2,50));

5 秒間隔で波動方程式の解をプロットします。

figure;
for k = 1:10
    subplot(5,2,k);
    pdeplot(model,'XYData',results.NodalSolution(:,k*5))
    axis equal
end  

複数の PDE からなる方程式

2 つの 2 階 PDE からなる方程式を解きます。この PDE 系を解くには、pdeCoefficients を使用して PDE の係数をシンボリックに抽出し、pdeCoefficientsToDouble を使用してその係数を倍精度数に変換し、specifyCoefficients を使用して PDE モデルの係数を指定します。

この PDE 系は、固定された構造板に一様な圧力負荷をかけたときの構造板のたわみを表します。従属変数 $$u_1$$ および $$u_2$$ をもつ PDE 系は、次の式で与えられます。

$$-{\nabla }^2{u}_1+u_2=0$$,

$$-D{\nabla }^2{u}_2=p$$,

ここで、$$D$$ は板のたわみ剛性で、次の式で与えられます。

$$D=\frac{E{h}^3}{12(1-{\nu }^2)}$$,

また、$$E$$ は弾性率、$$\nu$$ はポアソン比、$$h$$ は板の厚さ、$$u_1$$ は板の横方向のたわみ、$$p$$ は圧力負荷です。

2 つの連立方程式の PDE モデルを作成します。

model = createpde(2);

矩形形状を作成します。矩形の側面の長さを指定します。次に、この幾何形状を PDE モデルに組み込みます。

len = 10.0;
gdm = [3 4 0 len len 0 0 0 len len]';
g = decsg(gdm,'S1',('S1')');
geometryFromEdges(model,g);

方程式の物理パラメーターの値を指定します。外部の圧力 $$p$$ は、任意の値を取ることができるシンボリック変数 pres とします。

E = 1.0e6;
h_thick = 0.1;
nu = 0.3;
D = E*h_thick^3/(12*(1 - nu^2));
syms pres

PDE 系をシンボリック方程式として宣言します。PDE の係数を抽出し、その係数をシンボリック形式で返します。

syms u1(x,y) u2(x,y)
pdeeq = [-laplacian(u1) + u2; -D*laplacian(u2) - pres];
symCoeffs = pdeCoefficients(pdeeq,[u1 u2],'Symbolic',true)

symCoeffs = struct with fields:
    m: 0
    a: [2×2 sym]
    c: [4×4 sym]
    f: [2×1 sym]
    d: 0

係数 m、a、c、f、および d を表示します。

structfun(@disp,symCoeffs)

$$0$$

$$\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right)$$

$$\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& \frac{25000}{273}& 0\\ 0& 0& 0& \frac{25000}{273}\end{array}\right)$$

$$\left(\begin{array}{c}0\\ \mathrm{pres}\end{array}\right)$$

$$0$$

関数 subs を使用して pres の値を代入します。subs の出力はシンボリック オブジェクトであるため、Partial Differential Equation Toolbox への有効な入力となるように、係数を pdeCoefficientsToDouble 関数を使用して double データ型に変換します。

symCoeffs = subs(symCoeffs,pres,2);
coeffs = pdeCoefficientsToDouble(symCoeffs)

coeffs = struct with fields:
    a: [4×1 double]
    c: [16×1 double]
    m: 0
    d: 0
    f: [2×1 double]

PDE モデルの PDE 係数を指定します。

specifyCoefficients(model,'m',coeffs.m,'d',coeffs.d, ...
    'c',coeffs.c,'a',coeffs.a,'f',coeffs.f);

バネの剛性を指定します。4 つの辺すべてで分布するバネを定義し、境界条件を指定します。

k = 1e7;
bOuter = applyBoundaryCondition(model,'neumann','Edge',(1:4), ...
    'g',[0 0],'q',[0 0; k 0]);

幾何形状のメッシュ サイズを指定し、PDE モデルのメッシュを生成します。

hmax = len/20;
generateMesh(model,'Hmax',hmax);

モデルを解きます。

res = solvepde(model);

ノード位置における解にアクセスします。

u = res.NodalSolution;

板の横方向のたわみをプロットします。

numNodes = size(model.Mesh.Nodes,2);
figure;
pdeplot(model,'XYData',u(1:numNodes),'contour','on')
title 'Transverse Deflection'

板の中央における横方向のたわみを求めます。

wMax = min(u(1:numNodes,1))

wMax = 
-0.2763

この結果を、解析的に計算した、板の中央における横方向のたわみと比較します。

pres = 2;
wMax = -.0138*pres*len^4/(E*h_thick^3)

wMax = 
-0.2760

入力引数

pdeeq — シンボリック形式の PDE

シンボリック方程式 | シンボリック式 | シンボリック ベクトル

シンボリック形式の PDE。シンボリック方程式、シンボリック式、またはシンボリック ベクトルとして指定します。

u — PDE の従属変数

シンボリック関数

PDE の従属変数。シンボリック関数として指定します。u には、2 次元または 3 次元の定常変数または時間依存変数が含まれていなければなりません。たとえば、次のいずれかのステートメントを使用して変数 u を作成します。

• syms u(x,y)

• syms u(t,x,y)

• syms u(x,y,z)

• syms u(t,x,y,z)

出力引数

coeffs — double で演算を行う PDE の係数

倍精度数と関数ハンドルの構造体

double で演算を行う PDE の係数。関数 specifyCoefficients が必要とする倍精度数と関数ハンドルの構造体として返されます。この構造体のフィールドは、a、c、m、d、および f です。Partial Differential Equation Toolbox が必要とする形式で係数を解釈する方法については、以下を参照してください。

• c Coefficient for specifyCoefficients (Partial Differential Equation Toolbox)

• m, d, or a Coefficient for specifyCoefficients (Partial Differential Equation Toolbox)

• f Coefficient for specifyCoefficients (Partial Differential Equation Toolbox)

pdeCoefficients が関数ハンドルとして係数を返す際、関数ハンドルは入力引数として location と state の 2 つの構造体を取り、倍精度出力を返します。関数ハンドルは、コマンド ウィンドウで @makeCoefficient/coefficientFunction として表示されます。coeffs.f の関数ハンドルの式をシンボリック形式で表示するには、coeffs.f('show') を使用します。

生成されたすべての係数を specifyCoefficients で使用できない場合もあります。たとえば、m が非ゼロの場合、係数 d は特別な行列形式を取らなければなりません。詳細については、d Coefficient When m Is Nonzero (Partial Differential Equation Toolbox)を参照してください。

symCoeffs — シンボリック形式の PDE 係数

シンボリック式の構造体

シンボリック形式の PDE 係数。シンボリック式の構造体として返されます。

バージョン履歴

R2021a で導入

参考

syms | diff | laplacian | pdeCoefficientsToDouble | specifyCoefficients (Partial Differential Equation Toolbox)

トピック

• 非線形熱伝導に関する偏微分方程式の求解