|
ODE 数値解の新しいエラー制御
Maple 13 では、エラー制御のメカニズムが拡張され、添字を 1 とする変数に適用されるようになり、また、陽的 RK ペア法では通常は不可視な特定の一意性を持つクラスを検出できるようになりました。これは Maple 12 では実現していなかった機能で、結果的に、指定の許容誤差を大幅に超えてエラーが発生している場合でも計算が可能になっています。
エラー制御のメカニズムでは、低次元の多項式 (補間式) によって、解の各構成要素を許容誤差の範囲内にモデル化できるようにする 必要があります。
|
|
系の例
|
|
次に示す系の例は、当社の顧客から寄せられたものですが、dsolve,numeric を直接呼び出す計算方法はもう行いません。
|
>
|
dsys := {diff(o15(t),t) = 2*(70*o16(t)*o15(t)*sin(t15(t))+100*
cos(t15(t))*sin(t15(t))*o16(t)^2+50*cos(t15(t))*o15(t)^2*
sin(t15(t))-40*cos(t15(t))*T1(t)+170*sin(t15(t))*o16(t)^2+35*
o15(t)^2*sin(t15(t))+20*cos(t15(t))*T2(t)-68*T1(t)+100*cos(t15
(t))*o16(t)*o15(t)*sin(t15(t))+14*T2(t))/(100*cos(t15(t))^2-\
189), diff(o16(t),t) = -2*(70*o16(t)*o15(t)*sin(t15(t))+50*cos
(t15(t))*sin(t15(t))*o16(t)^2-20*cos(t15(t))*T1(t)+35*sin(
t15(t))*o16(t)^2+35*o15(t)^2*sin(t15(t))-14*T1(t)+14*T2(t))/(
100*cos(t15(t))^2-189), diff(t15(t),t) = o15(t), diff(t16(t),t
) = o16(t), diff(Lo15(t),t) = -(4.*(50.*cos(t15(t))*sin(
t15(t))+35.*sin(t15(t)))*o15(t)/(100.*cos(t15(t))^2-189.)+2.*(
70.*sin(t15(t))+100.*cos(t15(t))*sin(t15(t)))*o16(t)/(100.*cos
(t15(t))^2-189.))*Lo15(t)-(-140.*sin(t15(t))*o15(t)/(
100.*cos(t15(t))^2-189.)-140.*o16(t)*sin(t15(t))/(100.*cos(t15
(t))^2-189.))*Lo16(t)-Lt15(t), diff(Lo16(t),t
) = -(2.*(70.*sin(t15(t))+100.*cos(t15(t))*sin(t15(t)))*o15(t)
/(100.*cos(t15(t))^2-189.)+4.*(170.*sin(t15(t))+100.*cos(t15(t))*
sin(t15(t)))*o16(t)/(100.*cos(t15(t))^2-189.))*Lo15(t)-
(-140.*sin(t15(t))*o15(t)/(100.*cos(t15(t))^2-189.)-4.*(50.*cos(
t15(t))*sin(t15(t))+35.*sin(t15(t)))*o16(t)/(100.*cos(t15(
t))^2-189.))*Lo16(t)-Lt16(t), diff(Lt15(t),t)
= -(2.*(-50.*sin(t15(t))^2+50.*cos(t15(t))^2+35.*cos(t15(t)))*
o15(t)^2/(100.*cos(t15(t))^2-189.)+400.*(50.*cos(t15(t))*sin(
t15(t))+35.*sin(t15(t)))*o15(t)^2*cos(t15(t))*sin(t15(t))/
(100.*cos(t15(t))^2-189.)^2+2.*(-100.*sin(t15(t))^2+100.*cos(t15(t
))^2+70.*cos(t15(t)))*o16(t)*o15(t)/(100.*cos(t15(t))^2-189.)+
400.*(70.*sin(t15(t))+100.*cos(t15(t))*sin(t15(t)))*o16(t)*
o15(t)*cos(t15(t))*sin(t15(t))/(100.*cos(t15(t))^2-189.)^2+2.*
(170.*cos(t15(t))-100.*sin(t15(t))^2+100.*cos(t15(t))^2)*o16(t
)^2/(100.*cos(t15(t))^2-189.)+400.*(170.*sin(t15(t))+100.*cos(t15(
t))*sin(t15(t)))*o16(t)^2*cos(t15(t))*sin(t15(t))/(100.*cos(
t15(t))^2-189.)^2+80.*sin(t15(t))*T1(t)/(100.*cos(t15(t))^2-189.)+
400.*(-68.-40.*cos(t15(t)))*T1(t)*cos(t15(t))*sin(t15(t))/(100.*
cos(t15(t))^2-189.)^2-40.*sin(t15(t))*T2(t)/(100.*cos(t15(t))^2-\
189.)+400.*(14.+20.*cos(t15(t)))*T2(t)*cos(t15(t))*sin(t15(t))/(
100.*cos(t15(t))^2-189.)^2)*Lo15(t)-(-70.*cos(t15(t))*
o15(t)^2/(100.*cos(t15(t))^2-189.)-14000.*sin(t15(t))^2*o15(t)
^2*cos(t15(t))/(100.*cos(t15(t))^2-189.)^2-140.*o16(t)*cos(t15
(t))*o15(t)/(100.*cos(t15(t))^2-189.)-28000.*o16(t)*sin(t15(t)
)^2*o15(t)*cos(t15(t))/(100.*cos(t15(t))^2-189.)^2-2.*(-50.*sin(
t15(t))^2+50.*cos(t15(t))^2+35.*cos(t15(t)))*o16(t)^2/(100.*
cos(t15(t))^2-189.)-400.*(50.*cos(t15(t))*sin(t15(t))+35.*sin(
t15(t)))*o16(t)^2*cos(t15(t))*sin(t15(t))/(100.*cos(t15(t)
)^2-189.)^2-40.*sin(t15(t))*T1(t)/(100.*cos(t15(t))^2-189.)-400.*(-20.
*cos(t15(t))-14.)*T1(t)*cos(t15(t))*sin(t15(t))/(100.*cos(t15(
t))^2-189.)^2-5600.*T2(t)*cos(t15(t))*sin(t15(t))/(100.*cos(t15(t)
)^2-189.)^2)*Lo16(t), diff(Lt16(t),t) = 0, T1(t) =
piecewise(2.*Lo15(t)*(-68.-40.*cos(t15(t)))/(100.*cos(t15(t
))^2-189.)-2.*Lo16(t)*(-20.*cos(t15(t))-14.)/(100.*cos(t15(
t))^2-189.) <= 0,1,0 < 2.*Lo15(t)*(-68.-40.*cos(t15(t)))/(100.*
cos(t15(t))^2-189.)-2.*Lo16(t)*(-20.*cos(t15(t))-14.)/(100.
*cos(t15(t))^2-189.),-1), T2(t) = piecewise(2.*Lo15(t)*(14.+20.
*cos(t15(t)))/(100.*cos(t15(t))^2-189.)-28.*Lo16(t)/(100.*
cos(t15(t))^2-189.) <= 0,1,0 < 2.*Lo15(t)*(14.+20.*cos(t15(
t)))/(100.*cos(t15(t))^2-189.)-28.*Lo16(t)/(100.*cos(t15(t)
)^2-189.),-1)}:
ini := {o15(0) = 0, o16(0) = 0, t15(0) = 0, t16(0) = 0,
Lo15(0) = -2.39833305595243072, Lo16(0) = -6.71887507333190381, Lt15(0) = -1.29541249021244886, Lt16(0) = -3.50442072628983636}:
vars := [o15(t), o16(t), t15(t), t16(t), Lo15(t), Lo16(t), Lt15(t), Lt16(t), T1(t), T2(t)]:
|
この問題は、添字 1 の変数間で不連続な飛び越しが区分関数によって発生しているため、添字 1 の変数に対する多項式のエラー制御条件を満たすことができなかったために起きています。
|
>
|
dsn := dsolve(dsys union ini, vars, numeric);
|
|
| (1.1) |
|
Error, (in dsn) cannot evaluate the solution further right of 1.3631899, probably a singularity
| |
|
|
|
方法 1: 新しいエラー制御が無効
|
|
この場合は、Maple の前バージョンと同じ (不正な) 結果が返されます。
|
>
|
dsn := dsolve(dsys union ini, vars, numeric, interr=false);
|
|
| (2.1) |
|
![r1 := [t = 10., o15(t) = -2.65521291821903160, o16(t) = -1.01779342629280878, t15(t) = -4.09762284269594179, t16(t) = -3.12796915806187181, Lo15(t) = 7.91072208369697538, Lo16(t) = 25.1415481130806491, Lt15(t) = 36.5829765308151451, Lt16(t) = -3.50442072628984080, T1(t) = -1.00000000000000022, T2(t) = -1.00000000000000022]](/support/helpjp/helpview.aspx?si=2042/file01042/math124.png)
| (2.2) |
エラーを確認するために、最初の転換点の領域を調べます。
|
>
|
plots[odeplot](dsn,[t,T1(t)],t=1.363..1.3635);
|
エラーは O(1) で、問題の関数は 1.75 または -1.25 ではなく、本来は 1 または -1 になるはずです。
|
|
|
方法 2: 不連続な変数を「離散」変数として扱う
|
|
この方法では、不連続な変数を 離散 変数として扱い、イベントを使用して転換点に対処します。
短所: 作業量が増えることと、イベントに関する一定の知識が必要になります。
長所: 得られる解が正確であることです (3 種類の方法のうち最も精度が高い)。
系の残りの部分に T1(t) と T2(t) がないと仮定すると、エラーはこれらの変数にのみ隔離されることになりますが、今回はこれに当てはまる事例ではないことから、転換点が系の残りの部分にも波及するエラーを引き起こすでしょう。
イベント手法で正しく対処できる理由は、イベントベースの転換点が連続する 1 つの変数としてではなく離散的に処理されるため連続的な積分を飛び越し部分で行わずに、この地点が「停止後に再開する」箇所として処理されるためです。
したがって、最初の手順では、離散変数を使用して区分関数の条件文を書き換えることになります。
|
>
|
pw := indets(dsys,specfunc(anything,piecewise)):
nops(pw);
|
|
| (3.1) |
|
| (3.2) |
|
>
|
cond1 := lhs(op(1,pw[1]));
|
|
| (3.3) |
|
>
|
eval(eval(cond1,t=0),ini);
|
|
| (3.4) |
新しい離散変数 disc1(t) は最初 -1 ですが、cond1 が増加して零を超えると 1 に変化します。やがて cond1 が減少して零を超えた時点で -1 に変化します。
|
>
|
dvars:= disc1(t):
dini := disc1(0)=-1:
evts := [increasing(cond1), disc1(t)=1],
[decreasing(cond1), disc1(t)=-1]:
|
第 2 区分は次のようになります。
|
| (3.5) |
|
>
|
cond2 := lhs(op(1,pw[2]));
|
|
| (3.6) |
|
>
|
eval(eval(cond2,t=0),ini);
|
|
| (3.7) |
新しい離散変数 disc2(t) は最初は -1 ですが、cond2 が増加して零を超えると 1 に変化します。やがて cond2 が減少して零を超えた時点で -1 に変化します。
|
>
|
dvars:= dvars,disc2(t):
dini := dini, disc2(0)=-1:
evts := evts, [increasing(cond2), disc2(t)=1],
[decreasing(cond2), disc2(t)=-1]:
|
ここで、区分の条件文を系の離散変数で書き換えて解を求めます。
|
>
|
ndsys := subs(cond1=disc1(t),cond2=disc2(t),dsys):
|
|
>
|
nini := ini union {dini}:
|
|
>
|
dsn := dsolve(ndsys union nini, vars, numeric, discrete_variables=[dvars], events=[evts]);
|
|
| (3.8) |
|
![r2 := [t = 10., o15(t) = -2.65569380903979946, o16(t) = -1.01715167157836439, t15(t) = -4.09798438177018554, t16(t) = -3.12822535149047765, Lo15(t) = 7.91609608943657506, Lo16(t) = 25.1466535536435565, Lt15(t) = 36.5804794483301805, Lt16(t) = -3.50442072628983903, T1(t) = -.99999999999999989, T2(t) = -.99999999999999989, disc1(t) = 1., disc2(t) = 1.]](/support/helpjp/helpview.aspx?si=2042/file01042/math265.png)
| (3.9) |
最初の転換点をプロットします。
|
>
|
plots[odeplot](dsn,[t,T1(t)],t=1.363..1.3635);
|
0..10 の範囲の離散変数をプロットします。
|
>
|
plots[odeplot](dsn,[[t,T1(t)],[t,T2(t)+0.5]],t=0..10, numpoints=1000);
|
相違が許容誤差よりはるかに大きいことに注意してください。
|
>
|
map(rhs,r1)-map(rhs,r2[1..-3]);
|
|
![[0., 0.480891e-3, -0.641754e-3, 0.361539e-3, 0.256193e-3, -0.5374005e-2, -0.510544e-2, 0.249708e-2, 0., 0., 0.]](/support/helpjp/helpview.aspx?si=2042/file01042/math296.png)
| (3.10) |
|
|
|
方法 3: 問題の回避
|
|
この手法では、系から T1 と T2 を除去します。すると、O(1) の一意性が O(t) またはより優れた一意性 (これで飛び越しの不連続性は微分部に置かれます) に変化します。
これは方法 1 よりは精度が高く、方法 2 (しかし、方法 2 よりは簡単です)。また、関連性のある区分関数を解の関数としてプロットできます。
|
>
|
T12, n2dsys := selectremove(a->member(lhs(a),{T1(t),T2(t)}),dsys):
|
|
>
|
n2dsys := subs(T12,n2dsys):
|
|
>
|
dsn := dsolve(n2dsys union ini, remove(member,vars,{T1(t),T2(t)}),numeric);
|
|
| (4.1) |
|
![r3 := [t = 10., o15(t) = -2.65568401531980136, o16(t) = -1.01716098131132715, t15(t) = -4.09797928769869557, t16(t) = -3.12821680403911229, Lo15(t) = 7.91598878067716782, Lo16(t) = 25.1465862728621836, Lt15(t) = 36.5805308585457638, Lt16(t) = -3.50442072628983814]](/support/helpjp/helpview.aspx?si=2042/file01042/math331.png)
| (4.2) |
|
>
|
plots[odeplot](dsn,[[t,subs(T12,T1(t))],[t,subs(T12,T2(t)+0.5)]], 0..10, numpoints=1000);
|
|
>
|
map(rhs,r1[1..-3])-map(rhs,r3);
|
|
![[0., 0.471097e-3, -0.632445e-3, 0.356445e-3, 0.247646e-3, -0.5266697e-2, -0.503816e-2, 0.244567e-2, 0.]](/support/helpjp/helpview.aspx?si=2042/file01042/math347.png)
| (4.3) |
|
>
|
map(rhs,r2[1..-5])-map(rhs,r3);
|
|
![[0., -0.9794e-5, 0.9309e-5, -0.5094e-5, -0.8547e-5, 0.107308e-3, 0.6728e-4, -0.5141e-4, 0.]](/support/helpjp/helpview.aspx?si=2042/file01042/math354.png)
| (4.4) |
第 2 の手法と第 3 の手法の結果の相違は、デフォルトの許容誤差の範囲内にあり、最初の手法と第 3 の手法との間の相違に比べて桁違いに小さいことがわかります。
|
|
