1. 過剰決定系
最初の例として、1 つの従属変数 y(x) を持つ 2 種類の不等式の系を過剰決定します。
|
>
|
![sys1 := [x*(diff(y(x), x))^2*(diff(y(x), x, x))^2-2*x*(diff(y(x), x))*(diff(y(x), x, x))*y(x)*(diff(y(x), x, x, x))+x*y(x)^2*(diff(y(x), x, x, x))^2-y(x)*(diff(y(x), x, x))+(diff(y(x), x))^2 = 0, -(diff(y(x), x))*(diff(y(x), x, x))+y(x)*(diff(y(x), x, x, x))+2*y(x)^2*(diff(y(x), x, x))^2-4*y(x)*(diff(y(x), x, x))*(diff(y(x), x))^2+2*(diff(y(x), x))^4 = 0]](/support/helpjp/helpview.aspx?si=1448/file00448/math112.png)
|
|
![sys1 := [x*(diff(y(x), x))^2*(diff(y(x), `$`(x, 2)))^2-2*x*(diff(y(x), x))*(diff(y(x), `$`(x, 2)))*y(x)*(diff(y(x), `$`(x, 3)))+x*y(x)^2*(diff(y(x), `$`(x, 3)))^2-y(x)*(diff(y(x), `$`(x, 2)))+(diff(y(x), x))^2 = 0, -(diff(y(x), x))*(diff(y(x), `$`(x, 2)))+y(x)*(diff(y(x), `$`(x, 3)))+2*y(x)^2*(diff(y(x), `$`(x, 2)))^2-4*y(x)*(diff(y(x), `$`(x, 2)))*(diff(y(x), x))^2+2*(diff(y(x), x))^4 = 0]](/support/helpjp/helpview.aspx?si=1448/file00448/math115.png)
| (4.1) |
ここで次のように rifsimp を呼び出します。
|
>
|
|
|
>
|
|
|
![ans1 := TABLE([Case = [[y(x) <> 0, diff(y(x), `$`(x, 3))]], Solved = [diff(y(x), `$`(x, 2)) = (diff(y(x), x))^2/y(x)], Pivots = [y(x) <> 0]])](/support/helpjp/helpview.aspx?si=1448/file00448/math130.png)
| (4.2) |
ここで系がより単純な方程式に簡素化され (y(x) の場合は完全にゼロになるわけではありません)、Maple の dsolve で処理できるようになります。
|
>
|
![dsolve(ans1[Solved], {y(x)}, explicit)](/support/helpjp/helpview.aspx?si=1448/file00448/math140.png)
|
|
| (4.3) |
厳密解を得るために dsolve を支援するだけでなく、この簡素化された形式は、initialdata 関数と組み合わせて使用することで、式を解くのに必要な初期データを得ることができます。
|
>
|
![id1 := initialdata(ans1[Solved])](/support/helpjp/helpview.aspx?si=1448/file00448/math154.png)
|
|
![id1 := TABLE([Finite = [y(x[0]) = _C1, (D(y))(x[0]) = _C2], Infinite = []])](/support/helpjp/helpview.aspx?si=1448/file00448/math157.png)
| (4.4) |
この場合は、2 階 ODE で想定された結果が得られましたが、さらに複雑な PDE 系のための必須初期データを計算することもできます (詳細は initialdata を参照してください)。これで、initialdata で計算した型の初期条件を持つ縮小された系に対し、数値計算方法を正しく適用できるようになりました。
|
>
|
![dsolve(`union`(convert(ans1[Solved], set), {y(0) = 1, (D(y))(0) = 2}), {y(x)}, numeric, output = Array([0, .25, .5, .75, 1.0]))](/support/helpjp/helpview.aspx?si=1448/file00448/math168.png)
|
|
| (4.5) |
最終的には、この問題の局所解の形式べき級数 (Taylor 級数) を得ることもできます (詳細については rtaylor を参照してください)。
|
>
|
![tay_ser := rtaylor(ans1[Solved], order = 4)](/support/helpjp/helpview.aspx?si=1448/file00448/math180.png)
|
|
![tay_ser := [y(x) = y(x[0])+(D(y))(x[0])*(x-x[0])+(1/2)*(D(y))(x[0])^2*(x-x[0])^2/y(x[0])+(1/6)*(D(y))(x[0])^3*(x-x[0])^3/y(x[0])^2+(1/24)*(D(y))(x[0])^4*(x-x[0])^4/y(x[0])^3]](/support/helpjp/helpview.aspx?si=1448/file00448/math183.png)
| (4.6) |
前述の初期データを与えると、上記のように級数が完全に決定されます：
|
>
|
![tay_ser := rtaylor(ans1[Solved], id1, order = 4)](/support/helpjp/helpview.aspx?si=1448/file00448/math189.png)
|
|
![tay_ser := [y(x) = _C1+_C2*(x-x[0])+(1/2)*_C2^2*(x-x[0])^2/_C1+(1/6)*_C2^3*(x-x[0])^3/_C1^2+(1/24)*_C2^4*(x-x[0])^4/_C1^3]](/support/helpjp/helpview.aspx?si=1448/file00448/math192.png)
| (4.7) |
2. 矛盾する系
rifsimp は、系が矛盾していないか判別するのによく使用されます。
|
>
|
![sys2 := [(diff(u(x, y), x))^2+diff(u(x, y), y, y) = 0, u(x, y)^3-(diff(u(x, y), x))+x = 0]](/support/helpjp/helpview.aspx?si=1448/file00448/math201.png)
|
|
>
|
|
|
![TABLE([status = "system is inconsistent"])](/support/helpjp/helpview.aspx?si=1448/file00448/math208.png)
| (4.8) |
3. 拘束機械系
この例では、rifsimp を拘束機械系 (代数方程式、つまり DAE 系) のプリプロセッサとして使用します。Lagrange 法の公式を使用すると、以下の重力下における摩擦なし形状 Phi(x,y) のワイヤのビーズ運動の非ゼロ質量 m を求めることができます：
|
>
|
)(x(t), y(t)), m*(diff(y(t), t, t)) = -lambda(t)*(D[2](Phi))(x(t), y(t))-m*g]](/support/helpjp/helpview.aspx?si=1448/file00448/math224.png)
|
|
)(x(t), y(t)), m*(diff(y(t), `$`(t, 2))) = -lambda(t)*(D[2](Phi))(x(t), y(t))-m*g]](/support/helpjp/helpview.aspx?si=1448/file00448/math227.png)
| (4.9) |
空気抵抗を考慮せずに、重力下で降下する質量は次のように表されます：
|
>
|
|
|
| (4.10) |
ここでは振子を例に挙げ、ビーズは円を描くようにワイヤ上を移動します。
|
>
|
|
|
| (4.11) |
|
>
|
|
|
![[m <> 0, x(t)^2+y(t)^2-1 = 0, m*(diff(x(t), `$`(t, 2))) = -2*lambda(t)*x(t), m*(diff(y(t), `$`(t, 2))) = -2*lambda(t)*y(t)-m*g]](/support/helpjp/helpview.aspx?si=1448/file00448/math252.png)
| (4.12) |
このような拘束系を数値ソルバで表すことは非常に困難です。もちろん、極座標を使用した拘束を排除することはできますが、通常、複雑な拘束を受ける機械系にとってこれは不可能です。たとえば、拘束を別の関数に置き換えることでワイヤの形状を変えてみます。
|
>
|
|
|
| (4.13) |
(ページ Rif の参考文献およびパッケージ情報 (Reidら著: 1996 年) を参照してください)。振子はこのような系における代表的な例で、アプリケーションでの重要性が高いことから近年の中心的な研究対象です。
|
>
|
|
|
![SimpSys := TABLE([Case = [[x(t) <> 0, diff(x(t), t)], [(y(t)-1)*(y(t)+1) <> 0, diff(lambda(t), t)]], Solved = [diff(y(t), `$`(t, 2)) = -(2*lambda(t)*y(t)+m*g)/m, diff(lambda(t), t) = -(3/2)*(diff(y(t), t))*m*g, diff(x(t), t) = -y(t)*(diff(y(t), t))/x(t)], Pivots = [m <> 0, x(t) <> 0, y(t)-1 <> 0, y(t)+1 <> 0], Constraint = [m*(diff(y(t), t))^2+2*y(t)^2*lambda(t)-y(t)*m*g+y(t)^3*m*g-2*lambda(t) = 0, x(t)^2+y(t)^2-1 = 0]])](/support/helpjp/helpview.aspx?si=1448/file00448/math276.png)
| (4.14) |
|
>
|
|
|
![TABLE([Finite = [lambda(t[0]) = _C1, x(t[0]) = _C2, y(t[0]) = _C3, (D(y))(t[0]) = _C4, g = _C5, m = _C6], Pivots = [_C6 <> 0, _C2 <> 0, _C3-1 <> 0, _C3+1 <> 0], Infinite = [], Constraint = [_C6*_C4^2+2*_C3^2*_C1-_C3*_C6*_C5+_C3^3*_C6*_C5-2*_C1 = 0, _C2^2+_C3^2-1 = 0]])](/support/helpjp/helpview.aspx?si=1448/file00448/math283.png)
| (4.15) |
前述した例ほどの簡潔さはないものの、rifsimp によって初期条件で満たす必要があるさらに別の拘束 (出力 initialdata における Constraint および Pivots) が見つかり、時間導関数の方程式 lambda(t) を得ることができます。ここで Maple の dsolve[numeric] を使用すると前述した例のような数値解が得られます。
4. Lie 対称性決定
この例では、rifsimp で ODE の Lie 点の対称性を決定する方法を示します：
|
>
|
|
|
| (4.16) |
以下のようにして、Lie 点の対称性を求める PDE 系を得ます：
|
>
|
(x, y), y(x))), _y1)}](/support/helpjp/helpview.aspx?si=1448/file00448/math319.png)
|
|
| (4.17) |
rifsimp を使用してこの系を次のように大幅に簡約します：
|
>
|
|
|
![rifsys := TABLE([Solved = [diff(eta(x, y), `$`(x, 2)) = (-eta(x, y)*x+xi(x, y)*y)/x^3, diff(xi(x, y), x) = eta(x, y)/y, diff(eta(x, y), y) = eta(x, y)/y, diff(xi(x, y), y) = 0]])](/support/helpjp/helpview.aspx?si=1448/file00448/math332.png)
| (4.18) |
この時点で pdsolve を使用すると以下のように Lie 点の対称性を得ることができます：
|
>
|
![pdsolve(rifsys['Solved'])](/support/helpjp/helpview.aspx?si=1448/file00448/math340.png)
|
|
| (4.19) |
