pdsolve/numeric/method - 指定された解法で PDE の数値解を求める
使い方
pdsolve(PDE, conditions, numeric, meth, options)
パラメータ
PDE - 単独の二変数偏微分方程式
conditions - 初期条件および境界条件の集合またはリスト
numeric - 数値解を得るためにキーワード
meth - method=meth 形式の名前 ( 記号およびインデックス ) ; 解の計算に使用する方法の選択
options - ( オプション ) keyword = value 形式の方程式
|
説明
|
|
•
|
デフォルトの thete 法を使用する pdsolve/numeric のインターフェースに加え、教育目的で数値 pdsolve を扱うために追加機能を使用することができます。
|
•
|
11 の標準的な方法の使用を含むこの追加機能は、数値的な境界条件の指定、二段階の計算法の最初のスキームの指定、そして使用者が設定した算法への拡張があります。
|
•
|
これらの計算法を含む pdsolve/numeric の使用は、単独の時間上の一階 parabolic/hyperbolic PDE に制限されています。
|
•
|
デフォルトの計算法を含み、pdsolve/numeric/method は解のプロットのために使用することができるモジュール、およびそれらの値を設定する手続きを得るモジュールを返します。plot, animate, plot3d, value, settings モジュールの詳細は pdsolve/numeric をご参照ください。
|
•
|
method 引数は問題に使用する計算法を指定します。利用可能な計算法は以下のものとなります。:
|
ForwardTime1Space[backward/forward]
CenteredTime1Space[backward/forward]
BackwardTime1Space[backward/forward]
Euler/ForwardTimeCenteredSpace
CrankNicholson/CenteredTimeCenteredSpace
BackwardEuler/BackwardTimeCenteredSpace
Box[backward/forward]
LaxFriedrichs
LaxWendroff
Leapfrog
DuFortFrankel
|
|
オプション
|
|
•
|
pdsolve/numeric/method に使用可能な追加された計算法のためにオプションは以下となります。
|
'time' = name
'range' = l..r
'spacestep' = numeric
'timestep' = numeric
'bcpts' = integer
'optimize' = boolean/symbol
'numericalbcs' = set, list or expression
'startup' = symbol
'startupsteps' = integer
'startupnumbc' = set, list or expression
•
|
time=name, range=l..r, spacestep=numeric, timestep=numeric, bcpts=integer, optimize=boolean/symbol オプションはデフォルトの計算法と同様となります。
|
|
- time と range オプションは、一つの端点のみを含む境界条件 ( そうでなければ、時間変数 / 空間変数および区間が自動的に発覚するような条件 ) の時間変数と空間変数からなる一階の問題のために指定しなければなりません。
|
|
- spacestep オプションの値は、離散的な問題に使用する空間ステップの幅となります。デフォルトでは空間の区間の 1/20 となっています。
|
|
- timestep オプションの値は、離散的な問題に使用する時間のステップ幅となります。デフォルトでは spacestep の値となっています。
|
|
- bcpts オプションの値は、境界条件での導関数の近似で使用する点の数となります。
|
|
- optimize オプションは、解のモジュールの生成に使用される最適化のレベルを制御します。詳細は pdsolve/numeric をご参照ください。
|
•
|
numericalbcs=set, list or expression 引数は数値境界条件を表現します。これはいくつかの PDEs を含む計算方法の使用に対して要求されます。
|
|
例えば Euler 法は、空間の領域で 3 連続の点の解の値と関連します。この意味は、境界条件が空間の領域の 2 点で解の値を指定することになります。もし領域で PDE が一階であるならば、一つだけの自然な境界条件を指定することができ、スキームは直接使用することができません。
|
|
数値境界条件は不連続方程式の追加を許可し、空間のメッシュ上で、許可された計算法は PDE のために使用することができます。これらは、しばしば補外法および要求で追加された境界条件の付近での PDE の打ち切りから導かれることがあります。
|
|
pdsolve コマンドでは、数値境界条件の二つの形式が使用できます。
|
|
数値境界条件から指定される一つ目の形式は、入力した PDE の境界付近で適用された指定の ( おそらくより小型な ) 計算法の一つから得られます。
|
|
二つ目の形式は、極めて詳細な制御を認めますが、数値境界条件からの要求を、境界での従変数の離散的な値の間の関係を表す不連続方程式として設定します。この数値境界条件の形式は、より良い時間と空間のステップ幅で、問題の時間変数と空間変数に依存することができ、離散的な従変数の値で線形でなければなりません。
|
|
numericalbcs 引数は、以下の形式を用いて指定することができます。:
|
`numericalbcs` = [method, point]
`numericalbcs` = expression
`numericalbcs` = {expression, 'timestep'=symbol, 'spacestep'=symbol}
`numericalbcs` = {[method[1], point[1]], ..., [method[n], point[n]],
expression[1], ..., expression[m],
'timestep'=symbol, 'spacestep'=symbol}
|
pdsolve ( 例えば Euler 法 ) および特別に設定された計算法で許された方法を表現している method では、method, timestep,spacestep で生成されるスキームに使用する境界と補正を表している point が数値境界条件で時間と空間のステップ幅に使用する名前を指定し、expression 引数は数値境界条件の形式を表しています。
|
•
|
数値境界条件で指定された計算法のために、点の値は計算法で得られる離散化公式で使用される最も左か最も右の点となります。例えば、もし点が 'n'-1 ( 'n' はいくつかの記号 ) と 2 の幅をもつ計算法から指定された場合、境界条件は右境界で適用され、'n'-1 と 'n'-2 の不連続点での解の値は離散化公式で使用されます。他の例として、もし点が同じ計算法で 1 として指定された場合、境界条件は左境界で適用され、1,2 の不連続点での解の値が離散化公式で使用されます。
|
|
より詳細な例として、LaxFriedrichs 法を使用した u[t]+0.1*u[x]=0 なる PDE の解を考えます。この問題は空間変数について一階であることから、一つの自然な ( 左境界での ) 境界条件だけが適用できます。しかしながら、計算法は 3 の幅をもっていて、二つの条件が必要となるので、数値境界条件の指定が必要となります。
|
|
追加した引数 numericalbcs=[Box,n] は、必要とされる右境界へ追加する関係を得るために Box 法の使用を指定し、実際には数値境界条件で指定される表現を使用するための対応として、より複雑な指定 numericalbcs={(v[1,n]-v[0,n] + v[1,n-1]-v[0,n-1])/2/k + 0.1*(v[1,n]-v[1,n-1] + v[0,n]-v[0,n-1])/2/h =0, timestep=k, spacestep=h} を使用することになります。
|
•
|
数値境界条件によって指定された表現は非常に特殊な形式で設定しなければなりません。離散的な従変数の値は、空間インデックス (u[timeidx,spaceidx]) に従った時間インデックスの最初を含む、インデックス変数として与えられなければなりません。時間インデックスの値は、1 の値を新しい時間ステップでの値で表されたものとして、-1 から 1 の区間で与えることができます ( 注意 : 区間は一段法の 0 から 1 区間からなります。 ) 。空間インデックスの値は、'n' を数値境界条件上の任意定数の記号とすると、1 ( 左境界 ) から 'n' ( 右境界 ) の区間で与えることができます。もし空間インデックスで 'n' が与えられた場合、i=0,1,2,... で n-i 形式とならなければなりません ( 例えば、u[1,n/2] と u[1,2*n] は認められません。 ) 。
|
|
例えば、最初の内部点での解の値と同様な右境界上の解の値を強制するような数値境界条件は、numericalbcs = u[1,n]-u[1,n-1] として指定されます。左境界のための類似の条件は numericalbcs = u[1,1]-u[1,2] として指定されます。
|
|
独立変数の値は常にいくつかの記述を必要とします。空間変数は、同様の様式での空間変数のインデックスと同じような従変数の値のための離散的な空間変数として与えられなければなりません ( 例えば x[n-1] あるいは x[1] ) 。時間の値は 0 時間インデックスと相対的なもので与えられなければなりません。例えば、従変数の値 u[1,0] の時間のために、'timevar' の代わりに時間の値 'timevar' + 'timestep' を使用します。
|
|
手軽な例として、右境界の数値境界条件として PDE u[t]+u[x]-sin(x*t) の後退時間と後退空間の打ち切りの使用を考えます。使用者は、u[t]=(u[1,n]-u[0,n])/k と u[x]=(u[1,n]-u[1,n-1])/h として導関数の打ち切りを得ることができます。これは、t+k,x[n] で打ち切りの拡張点が、それらの値を含む離散化公式での時間変数と空間変数と交換できることによります。それゆえ、この数値境界条件の指定は numericalbcs = {(u[1,n]-u[0,n])/k+u[x]=(u[1,n]-u[1,n-1])/h-sin((t+k)*x[n]), timestep=k, spacestep=h} となります。
|
|
注意 : この条件は numericalbcs = [BackwardTime1Space[backward],n] として使用する計算法の形式を指定することができます。
|
|
数値境界条件はもし必要ならば設定しなければなりません。そうでなければ、Maple はエラーを返すことになります。
|
•
|
startup=symbol, startupsteps=integer, startupnumbc=set, list or expression 引数は、未来の時間ステップでの解の計算のための現在および過去の時間ステップにおける解の値を必要とするような方法である Leapfrog や DuFortFrankel のような二段解法で必要とされる追加段階の計算方法を指定します。これは、一つの時間ステップが既知 ( これは問題の初期条件となります。 ) であるような PDE の時間積分の最初のステップで問題となります。
|
|
startup 引数は二つの形式をとることができ、計算法が指定された形式、あるいは関数が指定された形式となります。計算法が指定された形式には、method を数値 pdsolve を含んで使用することが可能な任意の一段法とすると、startup = method として引数を指定しなければなりません。関数が指定された形式には、expressionを時間変数と空間変数について完全に指定された任意の関数とすると、startup = expression として引数を指定しなければなりません。もし startup 引数が設定されなかった場合は、計算法の形式はデフォルトの Theta スキームを含んで使用されます。
|
|
startupsteps 引数 ( デフォルトの値は 1 となっています。 ) は、startup が二段法よりも低い解度であった場合、解の解度をより良いものとなします。これは、主な計算法によって必要とされる追加段の計算に使用する二重のステップの数よりも大きいものの一つを指定します。最適なものは、例を通して説明します。
|
|
Leapfrog 法が主要な方法として使用され、Euler 法が標準的な方法として使われるような PDE を考えます。
|
|
* startupsteps=1 でステップ幅が k の場合、Euler 法は必要な追加段として用いる t[0]+k での解の値の計算に使用されます。
|
|
* startupsteps=2 の場合、Euler 法は t[0]+k/2 での解の値の計算に使用され、Leapfrog 法は必要な追加段として用いる t[0]+k での解の値の計算に使用されます。
|
|
* startupsteps=3 の場合、Euler 法は t[0]+k/4 での解の値の計算に使用され、Leapfrog 法は t[0]+k/2 および t[0]+k での解の値の計算に使用されます。
|
|
* 一般に startupsteps=n では、Euler 法は t[0]+k/2^(n-1) での解の値の計算に使用され、Leapfrog スキームは t[0]+k/2^(n-2), t[0]+k/2^(n-3), ..., t[0]+k/2, t[0]+k での n-1 回の解の値の計算に使用されます。t[0]+k での値は、必要な追加段として使用されます。
|
|
このアプローチは、解度の無駄が小さい、あるいは全く含まない起動時の方法として使用されるため、低い解度の計算法を許可します。
|
|
startupnumbc 引数は、起動時のスキームのために使用する数値境界条件を指定します。デフォルトでは、pdsolve の最初の試みは、起動時の解法と主要な解法と同様の数値境界条件を使用することで、起動時の計算法を得ようと試みることになります。もしこれがうまくいかない場合は、pdsolve は数値境界条件を含まない起動時の計算法を使用した起動時の解法を得ようと試みます。このオプションは、これらのアプローチのどちらも十分でない場合を補助するために設定します。例えば、もし主要な計算法が二つの数値境界条件を必要とし、起動時の計算法を一つだけ必要とした場合、このアプローチはうまくいかず、startupnumbc を使用しなければなりません。注意として、不正確な条件数を含む startupnumbc の指定、あるいは必要とされないオプションの使用に関しては、エラーによって結果を返します。
|
|
|
例
|
|
ForwardTime1Space[backward] 法を使用した左の境界条件を含む基本的な一成分の波動方程式
>
|
PDE1 := {diff(u(x,t),t)=-1/2*diff(u(x,t),x)};
|
| (3.1) |
>
|
IBC1 := {u(x,0)=sin(Pi*x), u(0,t)=-sin(Pi*t/2)};
|
| (3.2) |
>
|
smod1 := pdsolve(PDE1, IBC1, numeric, time=t, range=0..1,
method=ForwardTime1Space[backward]);
|
| (3.3) |
この問題の完全な解が得られます。t=1.25 での誤差関数を得ます。
>
|
esol1 := sin(Pi*(x-t/2));
|
| (3.4) |
>
|
errfunc := smod1:-value(u(x,t)-esol1,t=1.25,output=listprocedure);
|
| (3.5) |
>
|
errfunc := rhs(errfunc[3]);
|
| (3.6) |
区間での最大ノルム誤差を計算します。
>
|
maxerr := max(seq(errfunc(i/100),i=0..100));
|
| (3.7) |
ステップ幅を半分にすると誤差も半分になることが確認できます。
>
|
smod1:-settings(timestep=1/40, spacestep=1/40);
errfunc2 := smod1:-value(u(x,t)-esol1,t=1.25,output=listprocedure);
|
| (3.8) |
>
|
errfunc2 := rhs(errfunc2[3]);
|
| (3.9) |
>
|
maxerr2 := max(seq(errfunc2(i/100),i=0..100));
|
| (3.10) |
| (3.11) |
一成分の非同次波動方程式に CrankNicholson 法を使用します。これは数値境界条件が必要となります。
>
|
PDE2 := diff(u(x,t),t)+1/10*exp(-t/10)*diff(u(x,t),x)
= cos(Pi*(x-exp(-1/10*t)))*Pi*exp(-1/10*t)/5;
|
| (3.12) |
>
|
IBC2 := {u(x,0)=-sin(Pi*x), u(0,t)=u(2,t)};
|
| (3.13) |
数値境界条件として、最初の内部点を用いた右境界での PDE の Box 打ち切りを選択します。
>
|
discs := diff(u(x,t),t) = (u[1,n]-u[0,n] + u[1,n-1]-u[0,n-1])/2/k,
diff(u(x,t),x) = (u[1,n]-u[1,n-1] + u[0,n]-u[0,n-1])/2/h;
|
| (3.14) |
>
|
nbc := {eval(PDE2,{discs,t=t+k/2, x=x[n]-h/2}),
timestep=k, spacestep=h};
|
| (3.15) |
注意 : time と range の指定は必要としません。これらは暗黙のうちに周期的な境界条件で定義されます。
>
|
smod2 := pdsolve(PDE2, IBC2, type=numeric, numericalbcs=nbc,
method=CrankNicholson);
|
| (3.16) |
再び、問題の完全解を得ます。誤差の検討をします。
>
|
esol2 := sin(Pi*(x-exp(-1/10*t)));
|
| (3.17) |
>
|
errfunc := rhs(smod2:-value(u(x,t)-esol2,t=1,
output=listprocedure)[3]);
|
| (3.18) |
>
|
maxerr := max(seq(errfunc(i/50),i=0..100));
|
| (3.19) |
さらに、ステップ幅を半分にします。 :
>
|
smod2:-settings(timestep=1/40, spacestep=1/40);
errfunc2 := rhs(smod2:-value(u(x,t)-esol2,t=1,
output=listprocedure)[3]);
|
| (3.20) |
>
|
maxerr2 := max(seq(errfunc2(i/50),i=0..100));
|
| (3.21) |
| (3.22) |
予期したものよりも良く改良されています。
左境界の Box スキーム ( 明確に導かれたスキームは除きます。 ) は、numericalbcs=[Box,n] を用いて指定することができます。
>
|
smod2b := pdsolve(PDE2, IBC2, numeric, numericalbcs=[Box,n],
method=CrankNicholson);
|
| (3.23) |
>
|
errfunc := rhs(smod2b:-value(u(x,t)-esol2,t=1,
output=listprocedure)[3]);
|
| (3.24) |
>
|
maxerr := max(seq(errfunc(i/50),i=0..100));
|
| (3.25) |
二段スキームの例として、孤立点を含む convection/diffusion 方程式にDuFortFrankel 法を適用させます。
>
|
PDE3 := diff(u(x,t),t)=1/20*diff(u(x,t),x,x)+1/2*diff(u(x,t),x);
|
| (3.26) |
>
|
IBC3 := {u(x,0)=cos(Pi*x)^2, D[1](u)(-1/2,t)=0,D[1](u)(1/2,t)=0};
|
| (3.27) |
>
|
smod3 := pdsolve(PDE3, IBC3, type=numeric, method=DuFortFrankel,
startup=Euler, timestep=1/50);
|
| (3.28) |
異なる時間でのプロットの級数解から、いくつかのプロットで生成し表示することができます。
>
|
p1 := smod3:-plot(t=0, color=red):
p2 := smod3:-plot(t=1/6,color=maroon):
p3 := smod3:-plot(t=1/3,color=blue):
p4 := smod3:-plot(t=2/3,color=green):
p5 := smod3:-plot(t=1, color=yellow):
plots[display]({p1,p2,p3,p4,p5});
|
|
|