12. 数値計算

12.1 数値計算について
12.2 DCADRE
12.3 ELLIPT
12.4 FOURIER
12.5 NDIFFQ
12.6 数値計算に関する諸定義

12.1 数値計算について

12.2 DCADRE

以下は旧弊である。FORTRANライブラリと現在のMAXIMAとのインターフェイスを構築 する事は、"maxima/src/fortdef.lsp"の例題を参照されたい。

- Romberg積分のIMSL版は漸くMacsymaで利用可能である。文献は PRINTFILE(DCADRE,USAGE,IMSL1);、デモは、batch("dcadre.mc");を実行されたい。 これは数値積分パッケージで、Romberg外挿を用いたものである。 DCADREパッケージはIMSL FORTRANライブラリのルーチンDCADREを呼出す為に書かれて いる。この文書はその様なプログラム向けである。バグやコメントはKMPに送られたい。

このパッケージをロードする為には、次の行を実行せよ。

LOADFILE("imsl")$

このパッケージのデモは、次の行を実行せよ。

batch("dcadre.mc");

動作する関数は次の構文を持つ；

• IMSL_ROMBERG(fn,low,hi) ここでfnは1引数の関数；lowとhiは積分の下限と上限でなければならない。 fnは浮動点値を返す筈である。
• IMSL_ROMBERG(exp,var,low,hi) ここでexpはvar=lowからhiまでの領域で積分されるべき関数である。 そのexpを評価した結果は浮動点数値である。
• FAST_IMSL_ROMBERG(fn,low,hi) この関数は、エラーチェックを一切行わないが、IMSL_ROMBERG関数以上の 処理速度を達成するかもしれない。fnはLisp関数(又は変換されたMacsyma関数)で、 浮動点引数を許容し、常に浮動点値を返す事を仮定している。

• [SUCESSS,answer,error]、ここでanswerは積分の結果、errorは出力の絶対誤差を 評価したものでDCADREは以下のPURPOSEで記述されたもの、 又は
• [WARNING,n,answer,error]、ここでnは警告コード、answerは答、そしてerrorは 出力の絶対誤差を評価したものでDCAREは以下のPUPPOSEで記述されたものである。 次の警告コードが生じる；
  • 65 = 一つ又はそれ以上の特異点が続いて発生。
  • 66 = 或る部分区間(複数の場合もある)で、何時もの振舞いが認められない にも関わらず、単に評価された誤差が小さい為に、積分の評価として 採用する。
  又は、
• [ERROR,errorcode]、ここで、errorcodeはIMSLが生成したエラーコードである。 次のエラーコードが生じ得る；
  • 131 = 内部作業領域が不十分な為に失敗
  • 132 = 失敗。これは関数にあまりにも多くのノイズが(与えられた要求誤差に 対し)含まれているか、被積分関数の病的な挙動による。
  • 133 = RERRが0.1よりも大きいか、0.0よりも小さいか、又は計算機の機械精度 と比べてあまりにも小さい場合。

• ROMBERG_AERR [デフォルト値 1.0E-5] -- 答で望ましいの絶対誤差。
• ROMBERG_RERR [デフォルト値 0.0] -- 答で望ましい相対誤差。

注意： IMSLがエラーを出す場合、メッセージは利用者のコンソール上に表示され、 エラーの内容が述べられる(このエラーメッセージはIMSVERBOSEをFALSEに 設定する事で抑制しても良い)。

注意： ここで変換されたFortranルーチンを使う為、再帰的な呼出しをしてはなら ない。それ自身を呼出さないが、利用者はIMSL_ROMBERGの計算の最中に^Aを 押してはならず、同じパッケージを用いて別の計算を始めて結果を得る事を 期待してはならない事に注意しなければならない -- IMSL_ROMBERGは、それを投入した時に一つのプロジェクトが既に実行され ているかどうか文句を言う。これは最少問題の原因となる。

Purpose (IMSL文書を編集した版) ----------------------------------------------------

DCADREは次の問題を解く事を試みる。一引数の実数値関数F、二つの実数値AとBが 与えられて次を満す数値DCADREを見付ける。

| / B | [ | / B | ] | [ | [ | [ | ] | I F(x)dx - DCADRE | <= max [ ROMBERG_AERR, ROMBERG_RERR * | I F(x)dx | ] | ] | [ | ] | ] | / A | [ | / A | ]

計算手順(IMSLのドキュメントを編集した版)

このルーチンはある手法を用いており、その手法では、DCADREは与えられた積分区間 の適当に選ばれた部分区間上でF(X)の積分に対する評価の和として計算されたもので ある。積分区間それ自身で開始する事、最初にその様な部分区間として、与えられた 部分区間上でRomberg外挿が亨受可能な評価を見付けるまで用いられる。 もし、この試みが失敗すると部分区間は同じ長さの二つの部分区間に分割され、各々 に分けて考察される。プログラミングノート(IMSL文書の修正版)

• 1. DCADRE (IMSL_ROMBERG向けの変換されたFORTRANを元にしたもの)は多くの場合で、 不連続点と代数的不連続点を飛す事が出来る。詳しくは参考文献を見よ。
• 2. 相対誤差助変数ROMBERG_RERRは区間[0.0,0.1]の中になければならない。例えば、 ROMBERG_RERR=0.1は積分の誤差が小数点第一位の精度までとする事を意味する。 ここで、ROMBERG_RERR=1.0E-4とすれば、精度は小数点第四位となる。DCADREは 要求される相対的精度が満足出来ないと決定すれば、IERが133(ROMBERG_RERRは 十分大きくなければならないが、100.0に加えられた時に結果が100.0よりも大 きい(これは本当に小さな浮動小数に対しては機械演算の性質からTRUEであるとは 限らない))に設定する。
• 3. 絶対誤差助変数のROMBERG_AERRは非負でなければならない。ROMBERG_AERRに合理的 な値を与える為に、利用者は計算されている積分のおおよその大きさが判っていな ければならない。多くの場合で、AERR=0.0を用いる事で十分である。この場合、 要求される相対誤差のみが計算で満足されている。
• 4. 参考文献からの引用で"とても用心深い人がDCADREを認可するのは、IER [警告やエラーコード]が0か65の時に限る。単に合理的なだけの人はIERが66であれ ば従順に受け取る。冒険者はIERが131か132であってもDCADREを承認するが、大抵 は正しい。"IERが0でなかった時でも、DCADREは計算された最適の評価を返す。

この手法に関する文献は、 de Boor, Calr, "CADRE: An Algorithm for Numerical Quadrature," Mathematical Software (John R. Rice, Ed.), New York, Academic Press, 1971, Chapter 7. を見よ。

12.3 ELLIPT

AM(U,M) - 母数Mの振幅関数 AM1(U,M1) - 補母数M1の振幅関数 AM(U,M):=AM1(U,1-M); so use AM1 if M ~ 1 SN(U,M):=SIN(AM(U,M)); CN(U,M):=COS(AM(U,M)); DN(U,M):=SQRT(1-M*SN(U,M)^2); (これらの関数はこの様に定義されている。他のCD,NS等も同様に定義される。) 完全楕円積分 ELLIPTK(M) - 第一種完全楕円積分 ELLIPTK1(M1) - 第一種完全楕円積分、補母数付き ELLIPTK(M):=ELLIPTK1(1-M); so use if M ~ 1 ELLIPTE(M) - 第二種完全楕円積分 ELLIPTE1(M1) - 第二種完全楕円積分、補母数付き ELLIPTE(M):=ELLIPTE1(1-M); so use if M ~ 1

12.4 FOURIER

12.5 NDIFFQ

Define_Variable(N,0.3,FLOAT);
Define_Variable(H,0.175,FLOAT);
F(X,E):=(Mode_Declare([X,E],FLOAT),N*EXP(X)/(E+X^(2*H)*EXP(H*X)));
Compile(F);
Array([X,E],FLOAT,35);
Init_Float_Array(X,1.0E-3,6.85); /* Xに区間を設定 */
E[0]:5.0;                        /* 初期条件*/
Runge_Kutta(F,X,E);              /* 解を計算 */
Graph2(X,E);                     /* 解の表示  */

12.6 数値計算に関する諸定義

関数: FFT (real-array, imag-array)

Fast Fourier Transform. このパッケージはLOAD(FFT);で読み込まれる。IFT命令もあり、これは逆Fourier変換 向けである。これらあの関数は1又は2次元のFLOATING-POINTの配列での (複素)Fast Fourier変換を実行する。

これらの配列は次の形式である：

ARRAY(<ary>,FLOAT,<dim1>); 又は ARRAY(<ary>,FLOAT,<dim1>,<dim2>);

1次元配列に対しては

<dim1> = 2^n-1

2次元配列に対しては

<dim1>=<dim2>=2^n-1

(つまり、配列は平方である)。 (MACSYMAの配列の添数は0から始まる事から、上の二つ場合は2^nと(2^2)^2の配列の 要素がある。)このパッケージは又二つの他の関数、POLARTORECTとRECTTOPOLARを含む。 詳細はDESCRIBE(命令)を実行せよ。実装の詳細に関しては、 PRINTFILE(FFT,USAGE,SHARE); を実行せよ。

変数: FORTINDENT

デフォルト値:[0] - FORTRAN命令によって表示される式の左側のマージンを制御する。 0は通常の表示(つまり、6個のスペース)を与え、正の値は式が右側よりに表示される 事となる。

関数: FORTMX (name,matrix)

MACSYMA行列をFORTRANの割り当て命令の形式でname(i,j)=<関連する行列の元>の列に 変換する。この命令は今ではお古である。FORTMX(name,matrix);は今では FORTRAN(name=matrix);として実行して良い(もし、"name"に値が設定されていれば、 FORTRAN('name=matrix);が必要である)。FORTMX命令が除去されても良い様に、それを 用いたコードを変換して欲しい。

関数: FORTRAN (exp)

expをFORTRANの正規の式に変換し、各行と継続行に6個のスペースが挿入されて、 (Maximaの羃乗である)^は**となる。オプションのFORTSPACES[FALSE]がTRUEであれば、 FORTRAN命令はスペースで80行を満たす。もし、FORTRANが値の設定された原子に対し て呼出されると、例えば、ここでX:A*B$が実行されていた状態でFORTRAN(X);を実行 すればX={Xの値}、ここではX=A*Bが生成される。更に、例えば、M:MATRIX(...); が実行されていれば、FORTRAN(M);は適切な配列の割当て文の書式で name(i,j)=<関連する行列の元>を生成する。FORTINDENT[0]は表示される式の 左マージンを制御し、0は通常のマージン(つまり、6スペースのインデント)となり、 これを増やすと表示式はより右寄りになる。

変数: FORTSPACES

初期値:[FALSE] - TRUEであれば、FORTRAN命令でスペースを用いて80カラムを満す。

関数: HORNER (exp, var)

expをHornerの規則に沿って並び換えられた表現(Horner形式)に変換し、varが指定 されていれば、それを主変数として用いる。expのCRE形式の主変数が用いられて いれば、varを省略しても構わない。式expを数値的に評価する場合、HORNERを用いる 事で安定性が改善される事もある。FORTRANで用いるプログラムをMACSYMAで生成する 場合にも便利である。

(C1) 1.0E-20*X^2-5.5*X+5.2E20; 2 (D1) 1.0E-20 X - 5.5 X + 5.2E+20 (C2) HORNER(%,X),KEEPFLOAT:TRUE; (D2) X (1.0E-20 X - 5.5) + 5.2E+20 (C3) D1,X=1.0E20; ARITHMETIC OVERFLOW (C4) D2,X=1.0E20; (D4) 6.9999999E+19

関数: IFT (real-array, imag-array)

逆Fourier変換。このパッケージの読み込みはLOAD(FFT);を実行せよ。これらの関数 (FFTとIFT)は1次元、又は2次元のどちらかの浮動点配列で高速Fourier変換を実行し、 配列はARRAY(<ary>,FLOAT,<dim1>);やARRAY(<ary>,FLOAT,<dim1>,<dim2>);の形式 である。一次元配列では<dim1>は2^n-1に等しくなければならず、二次元配列では <dim1>=<dim2>=2^n-1でなければならない(つまり、配列は平方)。(MAXIMAの配列の 添字は0から開始し、上の二つの場合では、2^n個と(2^n)^2個の配列の要素となる)。 実装に関するより詳細はPRINTFILE(FFT,USAGE,SHARE);を実行せよ。

関数: INTERPOLATE (func,x,a,b)

Xをかえてfuncの零点を見付ける。最後の二つの引数は捜索の領域を与える。この関数 は各末端で符号が異っていなければならない。もし、この条件が満されていなければ、 この関数の動作はINTPOLERROR[TRUE]で支配される。もし、INTPOLERRORがTRUEであれ ば、エラーが発生し、それ以外はINTPOLERRORの値が返される(それ故、グラフ表示で は、INTPOLERRORは0.0に設定されるべきである)。それ以外(MACSYMAが第一引数を指定 された領域で評価可能で、しかも、それは連続として与えられている)では、 INTERPOLATEが零点を返す事が保証されている(もし、零点が一つよりも多ければ、 それらの内の一つを返す)。INTERPOLATEの精度はINTPOLABS[0.0]とINTPOLREL[0.0]で 制御され、これらは非負の浮動小数でなければならない。INTERPOLATEは最初の引数が INTPOLABS以下として評価された時や、根に対する次の近似解が INTPOLREL*<近似の一つ>と最早、違わなくなった場合に停止する。INTPOLABSと INTPOLRELのデフォルト値は0.0なので、INTERPOLATEが単精度演算で可能な限り 良い答が得られる。最初のargは方程式で良い。

最後の二つの引数の順序は無関係である。その為、

INTERPOLATE(SIN(X)=X/2,X,%PI,.1); は次と同値である。 INTERPOLATE(SIN(X)=X/2,X,.1,%PI);

用いられた手法は最後の二つの引数によって指定された領域での二分検索である。 関数が十分線形に近いと考えた時、線形補間の利用を開始する。別の構文が補間に 追加され、これは最初の二つの引数を関数名で置き代えるものである。関数は一つの 引数の関数をTRANSLATEで変換したものか、COMPILEで翻訳したものでなければなら ない。結果の検証は実行されないので、関数が返す浮動小数を確認せよ。

F(X):=(MODE_DECLARE(X,FLOAT),SIN(X)-X/2.0); INTERPOLATE(SIN(X)-X/2,X,0.1,%PI) time= 60 msec INTERPOLATE(F(X),X,0.1,%PI); time= 68 msec TRANSLATE(F); INTERPOLATE(F(X),X,0.1,%PI); time= 26 msec INTERPOLATE(F,0.1,%PI); time= 5 msec

Newton法補間ルーチンもあり、DESCRIBE(NEWTON);を実行せよ。

変数: INTPOLABS

デフォルト値:[0.0] - INTERPOLATE命令の精度はINTPOLABS[0.0]とINTPOLERL[0.0] によって制御され、それらは非負の浮動小数でなければならない。INTERPOLATEは 最初の引数の評価がINTPOLABS以下の時と、続く根の近似解が INTPOLREL*<近似解の一つ>よりも違わなければ停止する。INTPOLABSとINTPOLRELの デフォルト値は0.0で、INTERPOLATEは単精度演算で可能な限りの良い答えを得る。

変数: INTPOLERROR

デフォルト値:[TRUE] - INTERPOLATEの挙動を支配する。INTERPOLATEが呼出された時、 補間されるべき関数が補間区間の末端での関数値の符号が逆になる条件を満すか満さ ないかを決定する。もし、それらが逆の符号であれば補間が実行される。もし、 それらが同じ符号で、INTPOLERRORがTRUEであれば、エラーが出力される。もし、 それらが同じ符号で、INTPOLERRORがTRUEでなければ、INTPOLERRORの値が返される。 そこで、グラフ表示には、INTPOLERRORは0.0に設定すべきである。

変数: INTPOLREL

デフォルト値:[0.0] - INTERPOILATE命令の精度はINTPOLABS[0.0]とINTPOLREL[0.0]で 制御され、これらは非負の浮動点小数でなければならない。INTERPOLATEは最初の引数 がINTPOLABS以下であった時や、続く根の近似解がINTPOLREL*<近似解の一つ>よりも 違わなければ停止する。INTPOLABSとINTPOLRELのデフォルト値は0.0で、INTERPOLATE は単精度演算で可能な限りの良い答えを得る。

関数: NEWTON (exp,var,X0,eps)

SHAREディレクトリ上のファイルNEWTON 1はNewton法を用いた補間を実行する関数を 含む。LOAD(NEWTON);で使える。NEWTON法は、FLONUMに対して全てが評価される事を 要求する為にINTERPORATEが適用出来ないものに使える。そこで、 NEWTON(x^2-a^2,x,a/2,a^2/100); は、FLONUM*a^2<a^2/100であるかどうかが判らないと文句を言うが、ASSUME(a>0); (訳者注：a>0を仮定する)を実行し、それからNEWTONを再び実行すれば動作し、 x=a+<小さなflonum>*aが全てのその記号的処理で得られる。 INTERPOLATE(X^2-a^2,x,a/2,2*a);は.5*aがflonumではないと文句を言う…。 Newton-Cotesの8panel quadrature則を用いる適合型積分器はSHARE1:QQ FASLで利用 可能である。詳細はDESCRIBE(QQ)を実行せよ。

関数: POLARTORECT (magnitude-array, phase-array)

大きさと位相の(極座標)形式から、実と虚の(複素)形式に変換する。ここで、実部と 虚部は大きさと位相の配列で

<実部>=<大きさ>*COS(<位相>) ==> <虚部>=<大きさ>*SIN(<位相>)

と置かれる。

この関数はFFTパッケージの一部である。使う為にはLOAD(FFT);を実行せよ。FFTと IFTの様に、この関数は1次元、又は2次元の配列が使える。しかし、配列の次元が 2の羃乗である必要も2D配列が平方(正方行列)である必要も無い。

関数: RECTTOPOLAR (real-array, imag-array)

POLARTORECTを無効にする。相は-%PIから%PIの範囲で与えられる。この関数はFFT パッケージの一部である。使う為にはLOAD(FFT);を実行せよ。FFTとIFTの様に、 この関数は1次元、又は2次元の配列が使える。しかし、配列の次元は2の羃乗である 必要も2D配列が平方(正方行列)となる必要も無い。