|
高次時間オーダー PDE の一階への変換
|
|
説明
|
|
|
•
|
pdsolve/numeric ページで一階以上の導関数を含む時間変数をもった連立 PDE の数値解法の前に言及されているように、方程式系は解についての一階の方程式系に変換されます。
|
|
|
この変換の標準的な計算方法は、時間について一階以上の導関数の環境をもつ新しい変数を導入し、一階の方程式でこれらの新しい変数を定義することによって実行されます。
|
|
>
|
diff(u(x,t),t,t)=diff(u(x,t),x,x)+diff(u(x,t),x,t);
|
|
| (1.1) |
これは、以下で定義される新しい変数 v(x,t) を導入することで一階に変換されます。
|
| (1.2) |
この新しい変数を入力した方程式系で使用すれば、以下の一階の方程式系を結果として得ます。:
|
>
|
diff(u(x,t),t)=v(x,t),
diff(v(x,t),t)=diff(u(x,t),x,x)+diff(v(x,t),x);
|
|
| (1.3) |
|
•
|
さらなる例として、従変数の四次時間導関数を含む PDE について、一階の方程式系に変換するためには三つの新しい変数を導入する必要があります。
|
|
•
|
このページでは、このアプローチが率直に自動化されているにも関わらず、これらが常に最も良い解を得るわけではないことの説明をしていることが最も重要となります。しかしながら、他の解法では入力した PDE および連立 PDE についての改質が必要とされ、初期条件・境界条件を含み、しばしば手動で行うことが最適であることがあります。
|
|
•
|
以下の PDE と初期条件・境界条件について考えます ( これは端点で固く固定された棒に関する線形重調和方程式となります。 ) 。:
|
|
>
|
PDE := diff(u(x,t),t,t)=-1/10*diff(u(x,t),x,x,x,x);
|
|
| (1.4) |
|
>
|
Init := u(x,0)=uv(x),D[2](u)(x,0)=0;
|
|
)(x, 0) = 0](/support/helpjp/helpview.aspx?si=2485/file01485/math65.png)
| (1.5) |
|
>
|
Bdry := u(0,t)=0,D[1](u)(0,t)=0,
u(1,t)=0,D[1](u)(1,t)=0;
|
|
)(0, t) = 0, u(1, t) = 0, (D[1](u))(1, t) = 0](/support/helpjp/helpview.aspx?si=2485/file01485/math73.png)
| (1.6) |
|
|
自動的なアプローチは入力を以下のように再公式化します。
|
|
>
|
PDE1 := diff(u(x,t),t)=v(x,t),
diff(v(x,t),t)=-1/10*diff(u(x,t),x,x,x,x);
|
|
| (1.7) |
|
>
|
Init1 := u(x,0)=uv(x),v(x,0)=0;
|
|
| (1.8) |
|
>
|
Bdry1 := u(0,t)=0,D[1](u)(0,t)=0,
u(1,t)=0,D[1](u)(1,t)=0;
|
|
)(0, t) = 0, u(1, t) = 0, (D[1](u))(1, t) = 0](/support/helpjp/helpview.aspx?si=2485/file01485/math98.png)
| (1.9) |
|
|
もう一つの方法として、以下のようにも再公式化することができます。:
|
|
>
|
PDE2 := diff(v(x,t),t)=3/10*diff(u(x,t),x),
diff(u(x,t),t)=-1/3*diff(v(x,t),x,x,x);
|
|
| (1.10) |
|
>
|
Init2 := u(x,0)=uv(x),v(x,0)=0;
|
|
| (1.11) |
|
>
|
Bdry2 := u(0,t)=0,v(0,t)=0,
u(1,t)=0,v(1,t)=0;
|
|
| (1.12) |
|
|
これは全ての 'point'='value' 形式の ( 離散化誤差を全く含まない ) 境界条件にとって良い選択であったかもしれません。
|
|
>
|
PDE3 := diff(v(x,t),t)=3/10*diff(u(x,t),x,x),
diff(u(x,t),t)=-1/3*diff(v(x,t),x,x);
|
|
| (1.13) |
|
>
|
Init3 := u(x,0)=uv(x),v(x,0)=0;
|
|
| (1.14) |
|
>
|
Bdry3 := u(0,t)=0,D[1](u)(0,t)=0,
u(1,t)=0,D[1](u)(1,t)=0;
|
|
)(0, t) = 0, u(1, t) = 0, (D[1](u))(1, t) = 0](/support/helpjp/helpview.aspx?si=2485/file01485/math154.png)
| (1.15) |
|
|
これは方程式系が x の最も低い最大微分オーダーをもつことから、良い選択であったかもしれません。
|
|
|
固定されたステップサイズのための区間 t=0..5 での方程式系から得られた解の誤差の計算を考えます。
|
|
|
注意 : 近似解を得るために uv(x) を使用します。
|
|
>
|
apxsol := .504413*(cosh(4.73*x)-cos(4.73*x))*cos(7.07505*t)
+.495587*(sin(4.73*x)-sinh(4.73*x))*cos(7.07505*t);
|
|
| (1.16) |
|
>
|
uv(x) := eval(apxsol,t=0);
|
|
| (1.17) |
|
>
|
pds1 := pdsolve({PDE1},{Init1,Bdry1}, numeric,
timestep=1/32, spacestep=1/64);
|
|
| (1.18) |
|
>
|
pds2 := pdsolve({PDE2},{Init2,Bdry2}, numeric,
timestep=1/32, spacestep=1/64);
|
|
| (1.19) |
|
>
|
pds3 := pdsolve({PDE3},{Init3,Bdry3}, numeric,
timestep=1/32, spacestep=1/64);
|
|
| (1.20) |
|
>
|
pl1 := pds1:-plot([u(x,t)-apxsol,color=red], t=5):
pl2 := pds2:-plot([u(x,t)-apxsol,color=blue], t=5):
pl3 := pds3:-plot([u(x,t)-apxsol,color=green], t=5):
plots[display]({pl1,pl2,pl3});
|
|
>
|
apxsol := .5000168*(cosh(10.99561*x)-cos(10.99561*x))*cos(38.23301*t)
+0.5*(sin(10.99561*x)-sinh(10.99561*x))*cos(38.23301*t);
|
|
| (1.21) |
|
>
|
uv(x) := eval(apxsol,t=0);
|
|
| (1.22) |
|
>
|
pds1 := pdsolve({PDE1},{Init1,Bdry1}, numeric,
timestep=1/64, spacestep=1/64);
|
|
| (1.23) |
|
>
|
pds2 := pdsolve({PDE2},{Init2,Bdry2}, numeric,
timestep=1/64, spacestep=1/64);
|
|
| (1.24) |
|
>
|
pds3 := pdsolve({PDE3},{Init3,Bdry3}, numeric,
timestep=1/64, spacestep=1/64);
|
|
| (1.25) |
|
>
|
pl1 := pds1:-plot([u(x,t)-apxsol,color=red], t=5):
pl2 := pds2:-plot([u(x,t)-apxsol,color=blue], t=5):
pl3 := pds3:-plot([u(x,t)-apxsol,color=green], t=5):
plots[display]({pl1,pl2,pl3});
|
|
|
これは x の最も低い最大微分オーダーの選択が最初のケースにとって最適であり、境界条件で導入した低い誤差が二番目のケースにとって最適な選択であったかのようにみえます。
|
|
|
この振るまいは、問題、PDE 、初期条件・境界条件の様子に基づいて様々に変化します。
|
|
|
さらにいくつかの場合、低い微分オーダーをとることは、微分項が混ざっていたり、難しい初期条件・境界条件を使用する場合は難しいかもしれません。
|
|
|
