updates/Maple11/de

はじめる前に
新機能一覧
Maple ワークシートの作成
Mapleワークシートを共有
Maple ウィンドウのカスタマイズ
Maple システムのカスタマイズ
コネクティビティ
Mathematics
数学
物理
プログラミング
グラフィックス
Student パッケージ
Science and Engineering
科学・エンジニアリング
アプリケーションと例題
リファレンス
システム
MapleSim
MapleSim ツールボックス
MapleSim ヘルプ
Tasks
Toolboxes
エラーメッセージガイド
マニュアル
数学アプリ

Maple 11 における微分方程式ソルバの追加・改良点

サマリ

Maple 11 での微分方程式(DE)のテーマは、微分方程式の持つ幾何学的構造の探検を可能にすることです。例えば、2 階非線形常微分方程式に対しては、楕円関数を用いてすべてのクラスに関する解法を与えるアルゴリズムが開発されました。1 階の微分方程式(ODE)については、共形対称性を加味するのと同様に、非線形の Abel AIR クラスは超幾何関数および楕円関数に関して完全に求解可能です。偏微分方程式(PDE)については、このバージョンの目的でもあることから、幾何学的構造が関連する不変解の群同様にその対称性に反映されています。これは PDE 系の厳密解に対して、Maple ライブラリにとって過去 5 年間の中でもっとも先端的かつ重要な成果を意味しています。

厳密解

常微分方程式 (ODE)

2 階非線形方程式に対する楕円関数による新しい解の構成

新しくなった非線形 2 解微分方程式の求解アルゴリズムに、 WeierstrassP および 12 種の JacobiPQ 楕円関数に関する解表現が扱えるようになりました。これらの方程式のクラスは、1 点での対称性と、Maple 11 の新アルゴリズムの威力を発揮する問題にも関連する 1 階非線形 Abel 型 ODE への次数低減へとつながってきます。

例題

>	PDEtools:-declare(y(x), prime=x);

(2.1.1.1)

>	diff(y(x), x, x) = 4/xdiff(y(x), x) + 6/x^2y(x)^2 - 6/x^2y(x) - 1/2x^2*a;

`y''` = 4*`y'`/x+6*y^2/x^2-6*y/x^2-(1/2)*x^2*a

(2.1.1.2)

>	dsolve(%);

(2.1.1.3)

>	diff(y(x), x, x) = 2/xa/(x+a)diff(y(x), x) + 2k^2/x^2/(x+a)^4y(x)^3 + (-2x^3a+(-k^2-1-6a^2)x^2 - 6a^3x - 2a^4)/x^2/(x+a)^4y(x);

`y''` = 2*a*`y'`/(x*(x+a))+2*k^2*y^3/(x^2*(x+a)^4)+(-2*x^3*a+(-k^2-1-6*a^2)*x^2-6*a^3*x-2*a^4)*y/(x^2*(x+a)^4)

(2.1.1.4)

>	dsolve(%);

y = _C2*((k^2+1)/(1+k^2*_C2^2))^(1/2)*JacobiSN((-1/(x+a)+_C1)*((k^2+1)/(1+k^2*_C2^2))^(1/2), _C2*k)*x

(2.1.1.5)

上記の出力結果を得るためのアルゴリズムでは、その解が WeierstrassP または 12 種の JacobiPQ 関数において線形であっても、積分式を含まない解を求めることが可能です。

新オプション 'singsol = <none, essential, all>' と高次 1 階非線形常微分方程式の新しい楕円関数解

高次の方程式は高階微分で非線形なものとなり、それらの一般解に加えて、そのような方程式系は dsolve コマンドによって自動的に計算される特異解を与えます。

Maple 11 の目新しい機能は、2 つの要素から成ります: ユーザは singsol = ...; の形式の追加の引数を与えることで、全部(all), なし(none), または本質的な(essential)特異解を計算するか指定することができます。また、未計算積分や楕円関数の RootOf 式を持った形式的な分数ベキを計算する標準的な一般解の代わりに、dsolve は未計算積分や RootOf を含まない一般解を計算しますが、可能な場合においては JacobiPQ, WeierstrassP および WeierstrassPPrime などの楕円関数に関して線形かつ陽的な解を直接計算します。

例題

について 2 次の例題です: singsol = none のオプションのみを与えることで一般解を得ます(これにより dsolve が問題を場合分けしたり特異解を計算するための時間を節約します)

>	(-2/x^3y(x)+1/x^2diff(y(x),x))^2-1-b^2/x^8y(x)^4-(-b^2-1)/x^4y(x)^2 = 0;

(-2*y/x^3+`y'`/x^2)^2-1-b^2*y^4/x^8-(-b^2-1)*y^2/x^4 = 0

(2.1.2.1)

>	dsolve(%, singsol = none);

(2.1.2.2)

高階な ODE に対する特異解は、関連する高次の多項式の根となります; 以下はオプション引数 implicit を与えて、実際に根を計算する代わりに陰的な特異解を表現している多項式です。

>	(-2/x^3y(x)+1/x^2diff(y(x),x))^2+1/x^2y(x)a-4/x^6*y(x)^3+b = 0;

(-2*y/x^3+`y'`/x^2)^2+y*a/x^2-4*y^3/x^6+b = 0

(2.1.2.3)

>	dsolve(%, implicit); # implicit avoids computing the roots of the cube

(2.1.2.4)

について、その解が WeierstrassPPrime を含むような 3 次の例題です。

>	-(x+a)^6(y(x) + xdiff(y(x),x))^3 + 1/2a^3/b^3 - 27/2x^4y(x)^4 - 27b/cx^2y(x)^2 - 27/2b^2/c^2 - 3/2a/b(x+a)^4(y(x)+x*diff(y(x),x))^2 = 0;

-(x+a)^6*(y+`y'`*x)^3+(1/2)*a^3/b^3-(27/2)*x^4*y^4-27*b*x^2*y^2/c-(27/2)*b^2/c^2-(3/2)*a*(x+a)^4*(y+`y'`*x)^2/b = 0

(2.1.2.5)

>	dsolve(%, implicit);

y^4-(1/27)*(-54*b^4*x^2*y^2*c+a^3*c^2-27*b^5)/(x^4*b^3*c^2) = 0, y = WeierstrassPPrime(1/(x+a)+_C1, a/b, b/c)/x

(2.1.2.6)

>	map(odetest, [%], %%); # verify these results for correctness

(2.1.2.7)

AIR クラスの 1 階 Abel 型方程式に対する超幾何関数に関する新しい解

3 パラメータ AIR クラスの 1 階 Abel 型 ODE 系の新しい解法アルゴリズムにより、 F, P, Q, R, S の有理関数写像 x -> F(x), y(x) -> (P(x) y(x) + Q(x))/(R(x) y(x) + S(x)) による変数変換の元である種の代表元との同値関係を解くことで、有理係数を持つ方程式のすべてのクラスを解くことが可能です。加えて、制限されたクラスについては、実装されたアルゴリズムで対応する ODE の係数が非有理数に写像するような場合でも計算が可能です。Abel AIR クラスは、文字通り Maple 11 により現在すべてのクラスが可解なものとして表現されます。

例題

非有理係数を含む 2 つの例題です。

>	diff(y(x),x) = y(x)/x/(ln(x)y(x)+(-ln(x)+b)(-ln(x)+c));

(2.1.3.1)

>	dsolve(%);

(2.1.3.2)

陰的な解が HeunC 関数により解ける例です。

>	diff(y(x),x) = y(x)(y(x)-1)exp(x)/(y(x)a+exp(2x)-exp(x)*c);

`y'` = y*(y-1)*exp(x)/(y*a+exp(2*x)-exp(x)*c)

(2.1.3.3)

>	sol := dsolve(%);

sol := _C1+(HeunCPrime(0, -c, -1, 0, a+1/2, y/(y-1))*y-exp(2*x)*(y-1)*HeunC(0, -c, -1, 0, a+1/2, y/(y-1))/exp(x))*(y-1)^c/(((y-1)*(-exp(2*x)/exp(x)+c)*HeunC(0, c, -1, 0, a+1/2, y/(y-1))+y*HeunCPrime(0, c, -1, 0, a+1/2, y/(y-1)))*y^c) = 0

(2.1.3.4)

この解を超幾何関数として見るには convert(sol, hypergeom); を実行してください; (convert/to_special_function を参照) 以下は、超幾何関数を用いて直接的に表現される解を持つクラスの例です。

>	diff(y(x),x) = ((3+3x^4+3x^2-6x^3+6x)y(x)^2 + (3x^5+30x^4-26x^6-3x^3-15x^7+2x^8)y(x) - x^7(29x^3-31-27x+x^5+21x^2+5x^4))/((x^4+x^2+1)y(x) + x^3(x^3+x-3x^2+x^5-1-x^4))/x;

`y'` = ((3+3*x^4+3*x^2-6*x^3+6*x)*y^2+(3*x^5+30*x^4-26*x^6-3*x^3-15*x^7+2*x^8)*y-x^7*(29*x^3-31-27*x+x^5+21*x^2+5*x^4))/(((x^4+x^2+1)*y+x^3*(x^3+x-3*x^2+x^5-1-x^4))*x)

(2.1.3.5)

>	dsolve(%);

(2.1.3.6)

共形対称性を許容した 1 階の方程式に対する新しい解

求積法で計算される解を持つことのできる共形対称な非線形 1 階 ODE についての最小の考察は、Sophus Lie の独自の研究でした。ある共形対称の無限小生成作用素 [(x, y), (x, y)] は、条件 [x] = [y] かつ [y] = -[x] を満たします。ここで、 [x] はに関する微分を意味します。Maple ライブラリのフレームワークにおける主な考察は、この種の対称性を持つ ODE などは、体系的に積へ分離可能な無限小を許容した方程式系に変換され得るという問題 - これは Maple 6 で実装された解法アルゴリズム - です。

例題

Kamke's の例題 1.83

>	diff(y(x),x) = tan(x*y(x));

(2.1.4.1)

>	dsolve(%, implicit);

(2.1.4.2)

やや複雑な例

>	diff(y(x),x) = -(-2-(1/((1/2y(x) + 1/2x)^2) + c)y(x) + (1/((1/2y(x) + 1/2x)^2)+c)x+2I((1/2y(x) + 1/2x)^2 + b))/(2 - (1/((1/2y(x) + 1/2x)^2) + c)y(x) + (1/((1/2y(x) + 1/2x)^2) + c)x + 2I((1/2y(x) + 1/2x)^2+b));

`y'` = -(-2-(1/((1/2)*y+(1/2)*x)^2+c)*y+(1/((1/2)*y+(1/2)*x)^2+c)*x+(2*I)*(((1/2)*y+(1/2)*x)^2+b))/(2-(1/((1/2)*y+(1/2)*x)^2+c)*y+(1/((1/2)*y+(1/2)*x)^2+c)*x+(2*I)*(((1/2)*y+(1/2)*x)^2+b))

(2.1.4.3)

>	dsolve(%, implicit);

((1/2)*I)*exp(-(((1/2)*y+(1/2)*x)^2*c-1)/((1/2)*y+(1/2)*x))*(y-x)-(Int(exp(-(_a^2*c-1)/_a)*(_a^2+b), _a = `` .. (1/2)*y+(1/2)*x))-_C1 = 0

(2.1.4.4)

>	odetest(%, %%); # verifies this result, see ?odetest

(2.1.4.5)

3 次以上の線形 ODE に対する超幾何関数解の拡張された計算法

超幾何関数で表現可能な解を計算する 3 次以上の線形 ODE に対する解法アルゴリズムに改良が施されました。

>	diff(y(x),x,x,x,x) = -15/2x/(x^4+3)diff(y(x),x)-45/2x^2/(x^4+3)diff(y(x),x,x)-10x^3/(x^4+3)diff(y(x),x,x,x)+15/16*y(x)/(x^4+3);

`y''''` = -(15/2)*x*`y'`/(x^4+3)-(45/2)*x^2*`y''`/(x^4+3)-10*x^3*`y'''`/(x^4+3)+(15/16)*y/(x^4+3)

(2.1.5.1)

>	dsolve(%);

y = _C1*hypergeom([1/8, 3/8, 5/8, 7/8], [1/2, 3/4, 5/4], -3/x^4)/x^(1/2)+_C2*hypergeom([3/8, 5/8, 7/8, 9/8], [3/4, 5/4, 3/2], -3/x^4)/x^(3/2)+_C3*x^(1/2)*hypergeom([-1/8, 1/8, 3/8, 5/8], [1/4, 1/2, 3/4], -3/x^4)+_C4*hypergeom([5/8, 7/8, 9/8, 11/8], [5/4, 3/2, 7/4], -3/x^4)/x^(5/2)

(2.1.5.2)

偏微分方程式 (PDE): 可分性、対称解析、および不変解

自動的に計算される極小対称作用素や関連する厳密な不変式による解の構成などと同様に、PDE 系の古典的な対称解析を行うための 19 種の新コマンドが PDEtools に追加されました。新コマンドには、1 つの PDE 系を別の系に変換(簡約)するためのものも含まれています。separability コマンドは PDE 系の和または積へと可分であるかを決定するよう拡張されています。

PDE 系の可分性

例題

3 つの独立変数を持つ 2 つの未知関数による非線形 PDE 系 - まず出力表示の冗長表現を避けるために declare 機能を用います(これにより微分は添え字で表示)。

>	with(PDEtools);

(2.2.1.1)

>	declare((u,v)(x,y,z));

(2.2.1.2)

sys := [diff(u(x,y,z),x,y,z)*v(x,y,z) + diff(u(x,y,z),x,y)*diff(v(x,y,z),z) + diff(u(x,y,z),x,z)*diff(v(x,y,z),y) + diff(u(x,y,z),x)*diff(v(x,y,z),y,z) + diff(u(x,y,z),y,z)*diff(v(x,y,z),x) + diff(u(x,y,z),y)*diff(v(x,y,z),x,z) + diff(u(x,y,z),z)*diff(v(x,y,z),x,y) + u(x,y,z)*diff(v(x,y,z),x,y,z) = 0, diff(u(x,y,z),x,y,z) + diff(v(x,y,z),x,y,z) = 0];

sys := [u[x, y, z]*v+u[x, y]*v[z]+u[x, z]*v[y]+u[x]*v[y, z]+u[y, z]*v[x]+u[y]*v[x, z]+u[z]*v[x, y]+u*v[x, y, z] = 0, u[x, y, z]+v[x, y, z] = 0]

(2.2.1.3)

1 つの PDE の場合のように、連立 PDE 系についても特に引数を指定しない場合は和に関する可分性を確認します。

>	separability(sys);

(2.2.1.4)

これらの可分性の結果は、関連する与式の可積分条件だけで判定され、変数の分離に対するものではありません。以下は積による分離です;

>	remain := separability(sys, `*`);

remain := (u[z]^2*v^2+v[z]^2*u^2)*(u[x]^2*v^2+v[x]^2*u^2)*(u[z]*u[x]*v^2*u[y]+v[z]*v[y]*u^2*v[x])*v[y]*u[y], -(u[z]^2*v^2+v[z]^2*u^2)*(v[y]^2*u^2+u[y]^2*v^2)*(u[z]*u[x]*v^2*u[y]+v[z]*v[y]*u^2*v[x])*u[x]*v[x], -(u[x]^2*v^2+v[x]^2*u^2)*(v[y]^2*u^2+u[y]^2*v^2)*(u[z]*u[x]*v^2*u[y]+v[z]*v[y]*u^2*v[x])*u[z]*v[z]

(2.2.1.5)

上記の結果は、この 3 つの方程式が 0 でない限り分離可能な解が存在しないことを意味しています。すなわち、各変数それぞれに依存する 3 つの関数の積を含むような解は存在しないことを意味しています。これは、未知関数 {u(x, y, z) v(x, y, z)} のどちらかまたは両方が変数のいくつかに依存している(つまり、1 つ以上の変数に関して定数となる)ということを意味しているものではありません。例えば、次のような制限を与えて考えてみます;

>	restriction := [u(x,y,z) = u(x,z), v(x,y,z) = v(y,z)];

(2.2.1.6)

この制限の元でを評価してみます。

>	eval([remain], restriction);

(2.2.1.7)

つまり、この制限により可分性の条件が満たされていることがわかります。どんな場合でも、+ や * で可分な解は HINT オプションを使って pdsolve コマンドで計算できます。+ で可分な解を持つような特別な系の解は、3 つの各変数に依存した関数ととなります。

>	pdsolve(sys, HINT = `+`);

${v = _F1(x)+_F2(y)+_F3(z), u = _F4(x)+_F5(y)+_F6(z)}$

(2.2.1.8)

>	pdetest(%, sys); # verify results - see ?pdetest

(2.2.1.9)

>	pdsolve(sys, HINT = `*`);

(2.2.1.10)

>	map(pdetest, [%], sys);

(2.2.1.11)

対称解析と PDE 系の不変式

例題

PDEtools パッケージの新しい対称性のためのコマンドを用いて、対称性に関する“解析と求解”のサイクルを自動的に、もしくはステップ毎に実現できます。以下の例題を考えてみます。

>	declare(u(x,t));

(2.2.2.1)

>	PDE := diff(u(x,t),x,x) - diff(u(x,t), t) = 0;

(2.2.2.2)

determining PDE コマンドは解が極小作用素となる系を返します。

>	DetSys := DeterminingPDE(PDE);

$DetSys := {_xi[1][x] = (1/2)*_xi[2][t], _xi[1][t] = -2*_eta[1][u, x], _xi[1][u] = 0, _eta[1][x, x] = _eta[1][t], _eta[1][t, u] = -(1/4)*_xi[2][t, t], _eta[1][u, u] = 0, _xi[2][t, t, t] = 0, _xi[2][x] = 0, _xi[2][u] = 0}$

(2.2.2.3)

この系を pdsolve コマンドで解くことで、ユーザは PDE の点対称変換の極小生成作用素を計算することが可能です。

>	pdsolve(DetSys);

${_xi[2](x, t, u) = (1/2)*_C1*t^2+_C2*t+_C3, _xi[1](x, t, u) = (1/2)*(_C1*t+_C2)*x-2*t*_C4+_C9, _eta[1](x, t, u) = ((1/8)*(-2*t-x^2)*_C1+_C4*x+_C5)*u+exp(_c[1]^(1/2)*x)*_C8*exp(_c[1]*t)*_C6+_C8*exp(_c[1]*t)*_C7/exp(_c[1]^(1/2)*x)}$

(2.2.2.4)

もしくは、PDE を引数にして Infinitesimals コマンドを呼び出す方法もあります。

>	G := Infinitesimals(PDE);

(2.2.2.5)

極小作用は PDE 不変式の変換を計算できます。例えば上の結果の中から最初の極小作用を取り出し、新しい変数を v(r, s) とします。

>	SymmetryTransformation(G[1], u(x,t), v(r,s));

${v(r, s) = u, s = t+_epsilon, r = x}$

(2.2.2.6)

この変換を {x, t, u} について解きます。

>	TR := solve(%, {x,t, u(x,t)});

$TR := {t = s-_epsilon, u = v(r, s), x = r}$

(2.2.2.7)

この変換を PDE に適用することで、再度 PDE を得ます( dchange を参照)。

PDE;

(2.2.2.8)

>	dchange(TR, %, [r,s, v(r,s)]);

(2.2.2.9)

これらの極小作用は、相似変換( similarity transformation )とも呼ばれ、独立変数の個数を減らす変換の計算にも用いられます。これはより一般の不変式変換の特別な場合であります。例えば、先ほどの G の最初の要素を取り出して相似変換を計算してみます。

>	SimilarityTransformation(G[1], u(x,t), v(r,s), jetnotation = false);

(2.2.2.10)

上式で、新しい変数は v(s) です。この変換を {x, t, u(x, t)} について解きます。

>	solve(%, {x, t, u(x,t)});

(2.2.2.11)

以下は独立変数の個数を減らすための PDE における変数変換の結果です;

>	PDEtools:-dchange(%, PDE);

(2.2.2.12)

この代わりに、独立変数の個数を減らす不変変換(およびその逆変換)は、一度にすべての対称性を計算する InvariantTransformation コマンドを呼び出すことで、上で計算したのすべての対称性が得られます。

>	InvariantTransformation([G], u(x,t), v(r,s));

$[{r = t, v(r, s) = u}, {t = r, u = v(r, s)}], [{v(r, s) = u, r = x}, {u = v(r, s), x = r}], [{_p1 = t, r = t/x^2, v(r, s) = u}, {t = _p1, u = v(r, s), x = (_p1/r)^(1/2)}], [{_p1 = x, r = t, v(r, s) = u*exp((1/4)*x^2/t)}, {t = r, x = _p1, u = v(r, s)/exp((1/4)*_p1^2/r)}], [{_p1 = t, r = t/x, v(r, s) = u*x^(1/2)*exp((1/4)*x^2/t)}, {t = _p1, x = _p1/r, u = v(r, s)/((_p1/r)^(1/2)*exp((1/4)*_p1/r^2))}]$

(2.2.2.13)

例えば、上式の 2 番目のリストに変数変換を適用してみます。以下はその変換、新しい変数、縮約された PDE です。

>	%[2][2], map(rhs, [op(%[2][2])]);

(2.2.2.14)

>	dchange(%[1], PDE, %[2]);

(2.2.2.15)

この変換された PDE を解き(より簡単な ODE の問題に帰着)、元の PDE の解に変数を戻してみます。もしくは、ユーザは PDE の不変式を直接計算することで、これまでの対称性に基づくサイクルをすべて自動的に適用することが可能です。

>	InvariantSolutions(PDE);

u = _C1, u = x*_C1+_C2, u = _C1/(t^(1/2)*exp((1/4)*x^2/t)), u = _C1+erf((1/2)/(t/x^2)^(1/2))*_C2, u = (x*_C1+_C2*t)/(t*(t/x)^(1/2)*x^(1/2)*exp((1/4)*x^2/t))

(2.2.2.16)

Maple のこの新しい InvariantSolutions コマンドは、世の中にある記号計算システムの中でも非常にユニークなものです。PDE 系の独立変数のすべてについて縮約を自動的に行い、一回の計算で 1 つの対称性を計算するだけでなく、同時にすべての対称性を生成することで、解を導き出します。また、ユーザのより細かい要求に応えるのと同様に、解の(個々の変数の)依存性を InvariantSolutions コマンドに指示することも可能です。

数値解とプロット関連

•	ODE 系や微分代数方程式(DAE)の初期値問題の数値解をより効率的に計算するための改良が実装されました。これまでは難解であったり(Maple では)大規模であった問題にフォーカスされています。