|
pdsolve/numeric で使用可能な PDE の数値解法
|
|
説明
|
|
|
•
|
このページには、解度のオーダーの情報、数値境界条件が必要とされ、差分スキームが使用できるような PDE 問題について、時間依存型問題に利用できる全ての数値 PDE 解法の記述が含まれています。PDE を解くための数値解法についての更なる詳細は pdsolve/numeric/method をご参照ください。
|
|
|
全ての計算法は、時間変数に関して一階の変数係数非同次 PDE のために適用することが可能となります。以下では、h は空間ステップとして使用され、k は時間ステップとして使用されます。さらに、空間変数について PDE の微分のオーダーは xord によって表示されます。
|
|
|
ForwardTime1Space[forward/backward]
|
|
|
•
|
時間微分項を前進差分で離散化し、空間微分項を前進差分、あるいは後退差分で離散化する方法は、一階の時間 / 空間 PDEs の解を求めることができるような、陽の一段階法となります。
|
|
•
|
この計算法の解度は O(h,k) となっています。
|
|
•
|
右進行波によって表現される一階 PDEs は、後退差分法 ( ForwardTime1Space[backward] として指定されます。 ) を使用し、左境界としての境界条件を指定する必要があります。あるいは、左進行波が前進差分法 ( ForwardTime1Space[forward] として指定されます。 ) を必要とする場合は、右境界として境界条件を指定する必要があります。数値境界条件は必須ではありません。
|
|
•
|
これは陽のスキームなので、安定性には問題によるいくつかの a のための k < a*h 形式の制限が必要となります。
|
|
•
|
メッシュ点 [1,i] で計算される値のために使用されるスキームは、点 (t,x[i]) について PDE に差分を適用するという結果であり、後退差分法のための PDE に以下の離散化式を代入することで得ることができます。
|
u[t] = (u[1,i]-u[0,i])/k,
u[x] = (u[0,i]-u[0,i-1])/h,
u = u[0,i]
|
|
前進差分法の場合は、単純に u[x] のための離散化式を以下と交換します。
|
u[x] = (u[0,i+1]-u[0,i])/h.
|
|
|
CenteredTime1Space[forward/backward]
|
|
|
•
|
時間微分項を中央差分で離散化し、空間微分項を前進差分、あるいは後退差分で離散化する方法は、導関数 u,u[x],u[t],u[xt] を含む PDE の解を求めるために使用できる、陰の一段階法となります。
|
|
•
|
この計算法の解度は O(h,k^2) となっています。
|
|
•
|
右進行波によって記述された PDEs は、後退差分法 ( CenteredTime1Space[backward] として指定されます。 ) を使用し、左境界としての境界条件を指定する必要があります。あるいは、左進行波が前進差分法 ( CenteredTime1Space[forward] として指定されます。 ) を必要とする場合は、右境界として境界条件を指定する必要があります。数値境界条件は必須ではありません。
|
|
•
|
これは陰のスキームであり、多くの問題について無条件に安定しています ( チェックする必要があるかもしれません。) 。
|
|
•
|
メッシュ点 [1,i] で計算される値のために使用されるスキームは、点 (t+k/2,x[i]) について PDE に差分を適用するという結果であり、後退差分法のための PDE に以下の離散化式を代入することで得ることができます。
|
u[xt] = (u[1,i]-u[1,i-1]-u[0,i]+u[0,i-1])/k/h,
u[t] = (u[1,i]-u[0,i])/k,
u[x] = (v[1,i]-v[1,i-1]+v[0,i]-v[0,i-1])/2/h,
u = (u[1,i]+u[0,i])/2,
t = t+k/2
|
|
前進差分法の場合は、単純に u[xt] と u[x] のための離散化式を以下と交換します。
|
u[xt] = (u[1,i+1]-u[1,i]-u[0,i+1]+u[0,i])/k/h,
u[x] = (v[1,i+1]-v[1,i]+v[0,i+1]-v[0,i])/2/h.
|
|
|
BackwardTime1Space[forward/backward]
|
|
|
•
|
時間微分項を後退差分で離散化し、空間微分項を前進差分、あるいは後退差分で離散化する方法は、導関数 u,u[x],u[t],u[xt] を含む PDE の解を求めるために使用できる、陰の一段階法となります。
|
|
•
|
この計算法の解度は O(h,k) となっています。
|
|
•
|
右進行波によって記述された PDEs は、後退差分法 ( BackwardTime1Space[backward] として指定されます。 ) を使用し、左境界としての境界条件を指定する必要があります。あるいは、左進行波が前進差分法 ( BackwardTime1Space[forward] として指定されます。 ) を必要とする場合は、右境界として境界条件を指定する必要があります。数値境界条件は必須ではありません。
|
|
•
|
これは陰のスキームであり、多くの問題について無条件に安定しています ( チェックする必要があるかもしれません。) 。
|
|
•
|
メッシュ点 [1,i] で計算される値のために使用されるスキームは、点 (t+k,x[i]) について PDE に差分を適用するという結果であり、後退差分法のための PDE に以下の離散化式を代入することで得ることができます。
|
u[xt] = (u[1,i]-u[1,i-1]-u[0,i]+u[0,i-1])/k/h,
u[t] = (u[1,i]-u[0,i])/k,
u[x] = (v[1,i]-v[1,i-1])/h,
u = u[1,i],
t = t+k
|
|
前進差分法の場合は、単純に u[xt] と u[x] のための離散化式を以下と交換します。
|
u[xt] = (u[1,i+1]-u[1,i]-u[0,i+1]+u[0,i])/k/h,
u[x] = (v[1,i+1]-v[1,i])/h.
|
|
|
ForwardTimeCenteredSpace or Euler
|
|
|
•
|
時間微分項を前進差分で離散化し、空間微分項を中心差分で離散化する方法は、時間に関して一階で、空間に関して任意の階数の、混合偏導関数を含まない PDEs の解を求めることができる、陽の一段階法となります。
|
|
•
|
この計算法の解度は O(h^2,k) となります。
|
|
•
|
空間に関して奇数オーダーの PDEs は、各変数の条件の総数が同じ指定された数値境界条件を必要とし、境界条件の総数は空間変数の微分オーダーよりも大きいものの一つとなります。時間に関して二階の PDEs は数値境界条件を必要としませんが、各境界で同じ条件数を必要とします。
|
|
•
|
これは陽のスキームなので、安定性は空間ステップに比例して制限を時間ステップに強制します。一般に、このスキームは xord が奇数である PDEs に対して不安定です。空間に関して偶数オーダーの PDEs には、問題によるいくつかの同じ a のための k < a*h^xord 形式の制限が必要となります。
|
|
•
|
メッシュ点 [1,i] で計算される値のために使用されるスキームは、点 (t,x[i]) に 1+2*ceil(xord/2) 点を使用する点に関する PDE に中心差分を適用する結果となります ( このメッシュは中心のどちらかの側の ceil(xord/2) 点を使用します。 ) 。
|
|
|
|
CenteredTimeCenteredSpace or CrankNicholson
|
|
|
•
|
時間微分項を中心差分で離散化し、空間微分項も中心差分で離散化する方法は、時間に関して一階で、空間に関して任意の階数の、混合偏導関数を含む PDEs の解を求めることができる、陰の一段階法となります。
|
|
•
|
この計算法の解度は O(h^2,k^2) となります。
|
|
•
|
空間に関して奇数オーダーの PDEs は、各変数の条件の総数が同じ指定された数値境界条件を必要とし、境界条件の総数は空間変数の微分オーダーよりも大きいものの一つとなります。時間に関して二階の PDEs は数値境界条件を必要としませんが、各境界で同じ条件数を必要とします。
|
|
•
|
これは陰のスキームであり、多くの問題について無条件に安定しています ( チェックする必要があるかもしれません。) 。
|
|
•
|
メッシュ点 [1,i] で計算される値のために使用されるスキームは、点 (t+k/2,x[i]) に 1+2*ceil(xord/2) 点を使用する点に関する PDE に中心差分を適用する結果となります ( このメッシュは中心のどちらかの側の ceil(xord/2) 点を使用します。 ) 。
|
|
|
|
BackwardTimeCenteredSpace or BackwardEuler
|
|
|
•
|
時間微分項を後退差分で離散化し、空間微分項を中心差分で離散化する方法は、時間に関して一階で、空間に関して任意の階数の、混合偏導関数を含む PDEs の解を求めることができる、陰の一段階法となります。
|
|
•
|
この計算法の解度は O(h^2,k) となります。
|
|
•
|
空間に関して奇数オーダーの PDEs は、各変数の条件の総数が同じ指定された数値境界条件を必要とし、境界条件の総数は空間変数の微分オーダーよりも大きいものの一つとなります。時間に関して二階の PDEs は数値境界条件を必要としませんが、各境界で同じ条件数を必要とします。
|
|
•
|
これは陰のスキームであり、多くの問題について無条件に安定しています ( チェックする必要があるかもしれません。) 。
|
|
•
|
メッシュ点 [1,i] で計算される値のために使用されるスキームは、点 (t+k,x[i]) に 1+2*ceil(xord/2) 点を使用する点に関する PDE に中心差分を適用する結果となります ( このメッシュは中心のどちらかの側の ceil(xord/2) 点を使用します。 ) 。混同偏導関数を含まない PDEs に対しては、方向時間ステップが u[0,i] 形式、他の全ての点は u[1,*] 形式での離散化の点のみになります。
|
|
|
|
Box
|
|
|
•
|
ボックス法は、時間に関して一階で、空間に関して任意の階数の、混合偏導関数を含む PDEs の解を求めることができる、陰の一段階法となります。
|
|
•
|
この計算法の解度は O(h^2,k^2) となります。
|
|
•
|
この計算法は偶数個の点を含む空間メッシュを使用し、一つの境界には一つの条件が他よりも多くなければなりません。計算法は対称性 ( 一つの境界が特別な条件と対比した他の条件を持つ場合は、異なります。 ) を持つので、この計算法は境界条件が適応されるように自動的に調整します。
|
|
•
|
空間に関して偶数オーダーの PDEs は一つの追加された数値境界条件を必要とし、空間に関して奇数オーダーの PDEs は数値境界条件を必要としません。
|
|
•
|
これは陰のスキームであり、多くの問題について無条件に安定しています ( チェックする必要があるかもしれません。) 。
|
|
•
|
メッシュ点 [1,i] で計算される値のために使用されるスキームは、特別な境界条件 ( あるいは偶数オーダー問題のための数値境界条件 ) が右か左の境界で与えられた場合、異なる依存となります。もし特別な条件が左でのものならば、スキームは、i の左には点 (t+k/2,x[i]-h/2) に floor(xord/2) 点を使用し、右には floor(xord/2) 点を使用する点に関する PDE に中心差分を用いるアプリケーションを通して得ることができます。もし特別な条件が右でのものならば、スキームは、i の左には点 (t+k/2,x[i]+h/2) に 1+floor(xord/2) 点を使用し、右には 1+floor(xord/2) 点を使用する点に関する PDE に中心差分を用いるアプリケーションを通して得ることができます。
|
|
|
|
LaxFriedrichs
|
|
|
•
|
Lax-Friedrichs 法は、時間に関して一階で、空間に関して奇数オーダーの、混合偏導関数を含まない PDEs の解を求めることができる、陽の一段階法となります。
|
|
•
|
この計算法の解度は O(h^2,k,h^(2*ceil(xord/2))/k) となります。注意 : 計算法は空間に関して偶数オーダーの PDEs に適用することができますが、安定した必要条件 ( 下記をご参照ください。 ) は、O(h^(2*ceil(xord/2))/k) 項のために誤差が O(1) であることを強制します。
|
|
•
|
空間に関して奇数オーダーの PDEs は、各変数の条件の総数が同じ指定された数値境界条件を必要とし、境界条件の総数は空間変数の微分オーダーよりも大きいものの一つとなります。時間に関して二階の PDEs は数値境界条件を必要としませんが、各境界で同じ条件数を必要とします。
|
|
•
|
これは陽のスキームなので、安定性は空間ステップに比例して制限を時間ステップに強制します。時間に関して奇数オーダーの PDEs に対し、安定性は問題によるいくつかの同じ a のための k < a*h^xord 形式の制限が必要となります。この条件が計算法の解度を含んだ結合であった場合、誤差は O(h^2,a'*h^xord,h/a') であるか、あるいは k = a'*h^xord, a'<a である場合は単純に O(h) であることがわかります。
|
|
•
|
メッシュ点 [1,i] で計算される値のために使用されるスキームは、値 u[0,i] が、j <> i としたときの Euler ステンシル u[0,j] の他の全ての点を使用した値 u[0,i] での補間によって交換されるような例外を含む Euler スキームと全く同一のものとなります。
|
|
|
|
LaxWendroff
|
|
|
•
|
Lax-Wendroff 法は、時間と空間に関して一階の線形 PDEs の解を求めるために使用することができます。
|
|
•
|
この計算法の解度は O(h^2,k*h,k^2) となっています。
|
|
•
|
Lax-Wendroff 法は一つの境界条件が各境界に指定されているような数値境界条件を常に必要とします。
|
|
•
|
これは陽のスキームなので、安定性は問題によった同じ a のための k < a*h 形式の制限が必要となります。
|
|
•
|
このスキームの微分は複素微分であり、入力した PDE の導関数の差分スキームの計算を必要としますが、( 線形性の制限での結果 ) 、多くの標準的な参考書で利用することができます。
|
|
|
|
Leapfrog
|
|
|
•
|
Leapfrog 法は、時間に関して一階で、空間に関して任意のオーダーの、混合偏導関数を含まない PDEs の解を求めることができる、陽の二段階法となります。
|
|
•
|
この計算法の解度は O(h^2,k^2) となっており、二段階法であり、始動の方法を使用しなければなりません。
|
|
•
|
空間に関して奇数オーダーの PDEs は、各変数の条件の総数が同じ指定された数値境界条件を必要とし、境界条件の総数は空間変数の微分オーダーよりも大きいものの一つとなります。時間に関して偶数オーダーの PDEs は数値境界条件を必要としませんが、各境界で同じ条件数を必要とします。
|
|
•
|
これは陽のスキームなので、安定性は空間ステップに比例して制限を時間ステップに強制します。一般に、このスキームは xord が偶数である PDEs に対して不安定です。空間に関して奇数オーダーの PDEs には、問題によるいくつかの同じ a のための k < a*h^xord 形式の制限が必要となります。
|
u[t] = (u[1,i]-u[-1,i])/2/k
|
|
を使用して、メッシュ点 [1,i] での値の計算するために使用することができ、全ての特別な導関数の計算に 1+2*ceil(xord/2) 点での中心差分を使用することができます。
|
|
|
|
DuFortFrankel
|
|
|
•
|
DuFort-Frankel 法は、時間に関して一階で、空間に関して偶数オーダーの、混合偏導関数を含まない線形 PDEs の解を求めることができる、陽の二段階法となります。同じ場合で、非線形 PDEs の解を求めるためにも使用することができます。
|
|
•
|
この計算法の解度は O(h^2,k^2,k^2/h^(2*ceil(xord/2))) となっており、二段階法であり、始動の方法を使用しなければなりません。注意 : 計算法は空間に関して奇数オーダーの PDEs に適用することができますが、安定した必要条件 ( 下記をご参照ください。 ) は、O(k^2/h^(2*ceil(xord/2))) 項のために誤差が O(1) であることを強制します。しかしながら、注意として、もし全ての PDE について偶数オーダーの導関数が現れなかった場合は、結果のスキームは同様の PDE に対する Leapfrog 法によって生成されるものと同様のスキームとなるので、この場合はそれらのコメントが適用されます。
|
|
•
|
この計算法は空間に関して偶数オーダーの PDEs に制限されているので、数値境界条件は必要としませんが、各境界での条件の数は同じでなければなりません。
|
|
•
|
これは陽のスキームなので、安定性は空間ステップに比例して制限を時間ステップに強制します。この制限は問題によったいくつかの a についての k < a*h^xord 形式のものとなります。この条件が計算法の解度を含んで計算された場合、誤差は O(h^2,a'^2*h^xord) 、あるいは k = a'*h^xord, a'<a であるときに単純に O(h^2) となります。
|
|
•
|
メッシュ点 [1,i] での値を計算するために使用するスキームは、値 u[0,i] がスキームの全ての出来事のために、値 u[1,i] と u[-1,i] の平均によって交換されるような例外を含んだ Leapfrog スキームと全く同じものとなります。
|
|
|
|
独自の解法の構成方法
|
|
|
•
|
pdsolve で使用する新しい計算法およびスキームの構成は、やや入り組んでおり、基本的な Maple のプログラミングの知識を必要とします。本質的には、入力する PDE と返される情報、および入力した表での計算法自身の情報に対して、任意の点 [1,xi] に対する望む計算法を導くルーチンを構成しなくてはなりません。
|
|
|
- PDE は時間に関して一階でなければなりません。
|
|
|
- スキームは二つ以上の追加段を使用することができません。
|
|
|
- 時間ステップと空間ステップの幅は定数となります。
|
|
|
- スキームは計算法のステンシルの外側の空間変数の離散的な値を使用することができません。
|
|
•
|
以下は CenteredTime1Space 法のためのコードであり、必要な事柄を説明するための例となります。
|
`pdsolve/numeric/findif/CenteredTime1Space` := proc(
PDE, # The PDE to be discretized
u, # The dependent variable, with dependencies (e.g. u(x,t))
x, # The space independent variable
t, # The time independent variable
v, # The name to be used for discretization, v[timeidx,spaceidx]
ix, # The value of the central space index
h, # The space step size (symbolic)
k, # The time step size (symbolic)
parms, # List of additional parameters for the method (the index
# for the 'method' argument.
# e.g. for "method=CenteredTime1Space[forward]" the list
# is ['forward'].
INFO) # Input/Output table of information
local sb;
# First validate that this method works for the input PDE
if INFO["timeord"]<>1 or INFO["spaceord"]<>1 then
error("The CenteredTime1Space scheme can only be used with first-order PDE");
elif parms=[] or not member(parms,{['forward'],['backward']}) then
error("The CenteredTime1Space scheme requires a parameter, 'forward' or 'backward'");
end if;
# Set output parameters
INFO["directional"] := true;
# Now compute the discretization of derivatives that may
# appear in the PDE
sb := diff(u,t)=(v[1,ix]-v[0,ix])/k,u=(v[1,ix]+v[0,ix])/2;
if parms=['forward'] then
sb := sb,diff(u,x)=(v[1,ix+1]-v[1,ix] + v[0,ix+1]-v[0,ix])/2/h,
diff(u,x,t)=(v[1,ix+1]-v[1,ix]-v[0,ix+1]+v[0,ix])/k/h;
else
sb := sb,diff(u,x)=(v[1,ix]-v[1,ix-1] + v[0,ix]-v[0,ix-1])/2/h,
diff(u,x,t)=(v[1,ix]-v[1,ix-1]-v[0,ix]+v[0,ix-1])/k/h;
end if;
# Now evaluate the PDE with respect to the discretization schemes and
# adjust explicit occurrences of the time/space variables
eval(PDE,{sb,x=x[ix],t=t+k/2});
end proc:
|
|
入力データは例のパラメータのリストで表記されたコメントによって公正に説明されており、数項目については議論する必要があります。
|
|
|
計算法に対する手続きの名前は、pdsolve による 'method' オプションの引数に基づいたルーチンの存在の確認として、`pdsolve/numeric/findif/<method name>` 形式でなければなりません。
|
|
|
入力する PDE は、pdsolve によって数値的に解かれた PDE と同等な表現 ( 方程式ではありません。 ) となります。もし pdsolve のための入力が方程式であった場合、表現に変換されます。これは diff 表記 ( 作用表記ではありません。 ) によって変換され、仮に D[1](u)(x,t) が入力した PDE で表示された場合は、このルーチンは diff(u(x,t),x) の代わりとみなします。
|
|
|
'v' の入力は、離散化のインデックス名前を構成するために使用されるべき基本の名前となり、メッシュ点 [1,ix] での値は v[1,ix] として表記されるべきです。
|
|
|
返された表現は、'ix' を入力するパラメータとして、v[1,ix] での離散的な解の計算に使用されるものです。
|
|
|
INFO 表は、導かれた計算法が入力した PDE と離散化に関するいくらかの情報の出力に対して有効的であった場合、より簡単な決定に使用できる入力をもちます。さらに、"directional" フラグは出力として設定することができます。
|
|
•
|
INFO 表は、このルーチンで使用できる以下の入力をもちます。
|
|
|
INFO["timeord"] 入力は、時間変数を含む PDE の微分オーダーを与える整数をもちます。
|
|
|
INFO["spaceord"] 入力は、空間変数を含む PDE の微分オーダーを与える整数をもちます。
|
|
|
INFO["mixed"] 入力は、PDE の混同偏導関数が表示されている場合、これを表すブール値を持ちます。
|
|
|
INFO["directional"] フラグは、指定された未知変数の解法が重要である全ての陰の計算法のために、ルーチンによって 'true' と設定すべきです。そうでなければ、数値は境界条件の指定区域に依存する計算法に切り替わるかもしれません。もし設定しなければ、false と仮定されます。
|
|
|
例えば、一階 PDE に適用された Box 法は、境界条件が左で与えられた場合には、境界条件が右で与えられたものと同じ形式のものになります。結果として、このスキームのために、INFO["directional"] は false ( デフォルト ) とすべきです。
|
|
|
注意 : 陽の計算法は、これらが未来の時間ステップに対するスキームについての一つの離散的な従変数値の表記のみであるため、自動的に true の INFO["directional"] の設定をもちます。
|
|
•
|
あとは非常に簡単な事柄です。ルーチンは PDE の差分スキームを生成すべきであり、スキームの未知関数表記のみはメッシュ点 v[timeidx,spaceidx] 、ステップ幅 h,k 、0 時間ステップでの時間値 t 、およびメッシュ x[spaceidx] の空間変数での従変数の値とならなくてはなりません。
|
|
|
出力スキームのアプリケーションで使用された 't' の値は 0 時間ステップでのもので、もし離散化が t=tv+k での離散化した解の値の計算に使用されたものであった場合、時間値 t はスキームが tv.となる計算で使用されます。これは、半分のメッシュ点での差分として導かれるスキーム、あるいは 'future' 解が出力スキームの t の値を適用しなければならない点での時間値に使用する後退差分スキームを意味しています。後者の場合、これは t から t+k までの設定によって行うことができます。
|
|
|
pdsolve コマンドは、計算法の手続きが分解を実行できた中で離散化された PDE について、 `expand` または `normal`を使用した設定の最適化によって PDE の直接的な離散化に対してうまく実行できます。独自の解法を書き込んだ場合の実験に対して、最適なものとなります。
|
|
|
|
