| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
3.1 積分について 3.2 積分に関する諸定義
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
MACSYMAは積分処理の為に幾つかのルーチンを持っている。INTEGRATE命令はその内で 最も使われるものである。ANTIDパッケージもあり、これは指定していない関数の操作 (そしてその微分も勿論の事)を行う。数値計算の利用ではROMBERG関数とIMSL版の Romberg積分となるDCADREがある。適応的積分法もあり、これは Newton-Cotes 8 panel quadrature則を用いたQUANC8である。双曲関数は現在作成中で、 その詳細はDESCRIBE(SPECINT)を実行せよ。
(* 訳者注： IMSL版の積分関数はMaximaには含まれていない。 *) |
全般的にMACSYMAでは"初等関数"(有理関数、三角関数、対数、指数、根等)と多少の 拡張(error関数、dilogarithm)で可積分なものの積分処理のみとし、g(x)やh(x)の 様な未知の関数の積分は扱わない。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
(C1) 'INTEGRATE(%E**SQRT(A*Y),Y,0,4);
4
/
[ SQRT(A Y)
(D1) I %E dY
]
/
0
(C2) CHANGEVAR(D1,Y-Z^2/A,Z,Y);
Is A positive, negative, or zero?
pos;
0
/
[ ABS(Z)
2 I Z %E dZ
]
/
- 2 SQRT(A)
(D2) - ----------------------------
A
(C3) ev(%,'risch);
Is A positive, negative, or zero?
pos;
2 SQRT(A) 2 SQRT(A)
2 (- 2 SQRT(A) %E + %E - 1)
(D3) - ---------------------------------------------
A
|
CHANGEVARは和や積の添字の変更にも使って良い。しかし、和や積で変更を行った時、 この変更はシフト、つまり、I=J+...,であって、より次数の高い関数にはならない。 例えば、
(C3) SUM(A[I]*X^(I-2),I,0,INF);
INF
====
\ I - 2
(D3) > A X
/ I
====
I = 0
(C4) CHANGEVAR(%,I-2-N,N,I);
INF
====
\ N
(D4) > A X
/ N + 2
====
N = - 2
|
/B /S(X) | | | | F(X,Y) DY DX . | | /A /R(X) |
関数F(X,Y)は二変数関数でtranslateで変換されたものか、compileで翻訳されたもの でなければならない。更に、R(X)とS(X)も各々が1変数関数で変換か翻訳された関数で あり、AとBは浮動点小数でなければならない。このルーチンは二つの大域変数を持ち、 それらはxとyの区間の分割数:DBLINT_X,DBLINT_Yを決める。その両方の初期値は10で、 他の整数値(2*DBLINT_X+1点がx方向に計算され、y方向は2*DBLINT_Y+1となる)と独立 して変更する事か可能である。
このルーチンはX軸を分割し、Xの各値に対して最初にR(X)とS(X)を計算する。そして、 R(X)とS(X)の間のY軸を分割し、Y軸に沿って積分をSimpson則を用いて計算する。 それから、X軸に沿った積分を関数の値をYの積分でSimpson則を用いて計算する。 この手順は色々な理由で数値的に不安定であるが、それなりに速いものである。 しかし、高周波成分を持った関数や特異点(領域に極や分岐点)を持つ関数に対する 適用は避ける事。
Yの積分はR(X)とS(X)がどれだけ離れているかに依存し、距離S(X)-R(X)がXで急速に 変化すれば、様々なY積分で異なった刻幅を持つ切捨てによって重大な誤差が発生する かもしれない。
領域での収束性を改善する為にDBLINT_XとDBLINT_Yを増す事も可能であるが、計算 時間が増大する。関数値は保存されず、その関数が非常に多くの時間を費すもので あれば、何かを変更すると再計算で待つ必要がある(御免)。関数F,RとSの両方は DBLINTを呼出す為に、最初にtranslateで変換されたものかcompileで翻訳されたもの でなければならない。これは多くの場合で翻訳されたコードで最大の速度向上を目指 した結果である！ これらの補助的な数値の利用に関してはLPH(またはGJC)に尋ねよ。
(* 訳者注：
DEFINTを簡単に説明すると通常はINTEGRATE(LISPの命令としては$INTEGRATE)で用いら
れているSININTを内部で呼出して、その結果に上限と下限の値を代入している。
このSININTでもRISCH積分を行うRISCHINTを呼出す事もあるが、EVを用いてRISCHINTを
用いる様に指定する事が出来ない。この点を修正する為にはdefint.lispのANTIDERIVの
修正が最低でも必要となる。
DEFINTを用いた定積分で下限や上限に記号や式が含まれている場合、その正負を尋ねて
くる事がある。この場合、正であればpos;、負であればneg;、零であればzero;と入力
する。更に、assumeを用いて正負の指定を予め行っていれば、Maximaはこの様な質問
を行わない。
以下に例を示す；
(C1) integrate(sqrt(2*x-x^2),x,0,a);
Is a positive, negative, or zero?
pos;
2
ASIN(a - 1) + (a - 1) SQRT(2 a - a ) %PI
(D1) ------------------------------------ + ---
2 4
(C2) assume(a>0 and a<1);
(D2) [a > 0, a < 1]
(C3) integrate(sqrt(2*x-x^2),x,0,a);
2
(a - 1) SQRT(2 a - a ) - ASIN(1 - a) %PI
(D3) ------------------------------------ + ---
2 4
数値計算の手法を用いる場合には、ROMBERGの他にQUANC8がある。QUANC8の方が計算
速度、精度と安定性でROMBERGに勝っている事が多いが、ROMBERGの方が精度が良好な
場合もある。以下のsqrt(2*x-x^2)の積分の例ではROMBERGが精度で勝り、
exp(-x)*sin(x)の積分では計算速度でROMBERGが勝っている。
(C8) romberg(sqrt(2*x-x^2),x,0,1);
Evaluation took 0.21 seconds (0.21 elapsed)
(D8) .7853897937007632
(C9) quanc8(sqrt(2*x-x^2),x,0,1);
Evaluation took 0.04 seconds (0.04 elapsed)
(D9) .7849358178522697
(C10)
(C10) integrate(sqrt(2*x-x^2),x,0,1);
Evaluation took 0.13 seconds (0.13 elapsed)
%PI
(D10) ---
4
(C11) bfloat(%);
Evaluation took 0.00 seconds (0.00 elapsed)
(D11) 7.853981633974483B-1
(C12) quanc8(exp(-x)*sin(x),x,0.,1.);
Evaluation took 0.02 seconds (0.02 elapsed)
(D12) .2458370070002374
(C13) romberg(exp(-x)*sin(x),x,0.,1.);
Evaluation took 0.00 seconds (0.00 elapsed)
(D13) .2458370426035679
(C14) integrate(exp(-x)*sin(x),x,0,1);
Evaluation took 0.14 seconds (0.14 elapsed)
- 1 - 1
%E SIN(1) %E COS(1) 1
(D14) - ------------ - ------------ + -
2 2 2
(C15) bfloat(%);
Evaluation took 0.01 seconds (0.01 elapsed)
(D15) 2.458370070002374B-1
尚、QUANC8を用いる為には、Maxima-5.6ではqq.lspの修正が必要である。
*)
|
EXP(A*X+B)*COS(C*X)^N*SIN(C*X) |
(* 訳者注： このフラグとINTSCEルーチンはMaxima-5.6には含まれていない。 *) |
(C1) 'INTEGRATE(SINH(A*X)*F(T-X),X,0,T)+B*F(T)=T**2;
T
/
[ 2
(D1) I (SINH(A X) F(T - X)) dX + B F(T) = T
]
/
0
(C2) LAPLACE(%,T,S);
A LAPLACE(F(T), T, S)
(D2) ---------------------
2 2
S - A
2
+ B LAPLACE(F(T), T, S) = --
3
S
(C3) LINSOLVE([%],['LAPLACE(F(T),T,S)]);
SOLUTION
2 2
2 S - 2 A
(E3) LAPLACE(F(T), T, S) = --------------------
5 2 3
B S + (A - A B) S
(D3) [E3]
(C4) ILT(E3,S,T);
IS A B (A B - 1) POSITIVE, NEGATIVE, OR ZERO?
POS;
2
SQRT(A) SQRT(A B - B) T
2 COSH(------------------------)
B
(D4) F(T) = - --------------------------------
A
2
A T 2
+ ------- + ------------------
A B - 1 3 2 2
A B - 2 A B + A
|
(* 訳者注： Laplace変換に関してはLAPLACEの項を参照せよ。尚、分母を一次と二次の多項式に 因子分解する事は記号積分を確実に行う為の一つの有効な手法である。 *) |
(* 訳者注：
尚、integrateでRISCHの積分を強制的に行わせる事が可能である。これは以下の様に
EVを用いると良い。
(C1) ev(integrate(3^log(x),x),'risch);
LOG(3) LOG(x)
x %E
(D1) -----------------
LOG(3) + 1
但し、定積分の場合はこの手法は使えない。定積分でも同様の処理を行いたければ
defint.lisp内部で定義されたANTIDERIV関数を修正する必要がある。
*)
|
INTEGRATEはDEPENDS命令で設定されるDEPENDENCIES集合によって影響されない。 INTEGRATE(exp,var,low,high)はvarに対するlowからhighまでのexpの定積分を行うが、 計算出来なければ名詞型で返す。極限にはvarを含んでいてはならない。幾つかの手法 が用いられ、不定積分や経路積分では直接代入も含む。広義の積分では、正の無限大 に関してINFを使い、負の無限大にはMINFを使っても良い。積分"形式"が操作 (例えば、幾つかの助変数に関してある数値を代入するまで計算出来ない積分)に 望ましいものであった場合、名詞型の'INTEGRATEが利用されて積分記号と一緒に 表示される(以下の注1を見よ)。
(* 訳者注： このINFとMINFに関しては、-INFや-MINFは各々MINFやINFと同値なものではない事に 注意する。Maximaの広義の積分や、その他の代入操作でINFやMINFを用いても良いが、 -INFや-MINFを用いると全く無意味な結果を得る事がある。 *) |
関数LDEFINTはLIMITを積分の下限と上限の極限の評価に用いる。時々、積分計算中 に利用者に式の符号が何であるかを尋ねて来る。ここでの適切な応答は、正の符号 であればPOS;、零であればZERO;、負の符号の場合はNEG;である。
(C1) INTEGRATE(SIN(X)**3,X);
3
COS (X)
(D1) ------- - COS(X)
3
(C2) INTEGRATE(X**A/(X+1)**(5/2),X,0,INF);
IS A + 1 POSITIVE, NEGATIVE, OR ZERO?
POS;
IS 2 A - 3 POSITIVE, NEGATIVE, OR ZERO?
NEG;
3
(D2) BETA(A + 1, - - A)
2
(C3) GRADEF(Q(X),SIN(X**2));
(D3) Q(X)
(C4) DIFF(LOG(Q(R(X))),X);
d 2
(-- R(X)) SIN(R (X))
dX
(D4) --------------------
Q(R(X))
(C5) INTEGRATE(%,X);
(D5) LOG(Q(R(X)))
|
(注1) MACSYMAで積分が出来ない事が、直ちに常にその積分が閉形式で存在しない事を
意味しない。以下に示す例でintegrateは名詞型を返すが、その原始関数を簡単に見つ
ける事が可能である。例えば、X^3+X+1=0の根で被積分関数を以下の形式に
書換えると計算が出来る。
1/((X-A)*(X-B)*(X-C)) |
(C6) INTEGRATE(1/(X^3+X+1),X);
/
[ 1
(D6) I ---------- dX
] 3
/ X + X + 1
|
(* 訳者注： この様に有理式の積分で分母を一次式の式の積に分解しておく事は非常に有効である。 何故なら、この様な形式が積分のアルゴリズムで考慮されている為である。 *) |
(*訳者注：
Maxima-5.6では3^log(x)やsqrt(x+1/x-2)の積分が正常に行えない問題点を持っている。
最初の3^log(x)に関しては、RISCHを用いると正しい答を計算する。又、新しいもので
あれば虫は修正されている。ここでx^log(3)の積分であれば通常のintegrateで良い。
この3^log(x)の積分に関してはDEFINTで定積分を実行すると誤った答を返すが、これは
DEFINTの節で述べた様に、内部的にSININTが用いられ、SININTではこの種の関数の積分
が上手く扱えない為である。この問題は前述の様にMaximaの新しいバージョンで修正
されている。
(C10) integrate(3^log(x),x);
2 LOG(x)
3
(D10) ---------
2 LOG(3)
(C11) risch(3^log(x),x);
LOG(3) LOG(x)
x %E
(D11) -----------------
LOG(3) + 1
(C12) integrate(x^log(3),x);
LOG(3) + 1
x
(D12) -----------
LOG(3) + 1
(C13) trigsimp(d11-d12);
(D13) 0
この様にrischを用いた積分と関数をx^log(3)に変形したものの積分結果は一致する。
実際の積分もrischで計算したものが正しい。参考迄にMaxima-5.9preの結果を載せて
置く。
(C8) integrate(3^log(x),x);
1
(------ + 1) LOG(x)
LOG(3)
3
(D8) --------------------
1
(------ + 1) LOG(3)
LOG(3)
(C9) risch(3^log(x),x);
LOG(3) LOG(x)
x %E
(D9) -----------------
LOG(3) + 1
(C10)
次に、sqrt(x+1/x-2)の積分に関してはsqrt(factor(x+1/x-2))として積分を行えば
正しい答を返す。こちらの虫はMaxima-5.9preでも修正されていない。
ここの例では、xの値域を開区間(0,1)で考えているが、値域を指定せずに積分を行う。
(C14) integrate(sqrt(x+1/x-2),x);
3/2
2 x - 6 SQRT(x)
(D14) ------------------
3
(C15) integrate(sqrt(factor(x+1/x-2)),x);
/
[ ABS(x - 1)
(D15) I ---------- dx
] SQRT(x)
/
この様に、factorを用いなければ誤った答を返す。又、D15で示す様にXの値域が明確で
ない為に名詞型で結果が返されている。次に、assumeでxの値域を指定する。
(C16) assume(x<1 and x>0);
(D16) [x < 1, x > 0]
(C17) integrate(sqrt(factor(x+1/x-2)),x);
3/2
2 x - 6 SQRT(x)
(D17) - ------------------
3
Xの値域をASSUMEで仮定するとD17で示す正しい結果を得る。
これらの例で示す様に、Maximaの積分は正しい答を返すとは限らないが、内部処理を
適切に行う事で正しい答を得る事も可能な場合がある。これはMaximaに限った話では
無く、数式処理一般で記号積分の結果は面倒でも確認した方が良い。
*)
|
(C86) integrate(sin(x)=0,x);
/
[
(D86) - COS(x) = INTEGRATIONCONSTANT1 + I 0 dx
]
/
(C87) INTEGRATION_CONSTANT_COUNTER;
(D87) 1
|
EXP(A*X+B)*COS(C*X)^N*SIN(C*X)^M |
(* 訳者注： INTSCEはMaxima-5.6には含まれていない。 *) |
(* 訳者注：
minfとinfの扱いは注意しなければならない。特に、-infや-minfの様に符号を付けて
用いると本来の結果と同値でも無意味な結果を得るので注意が必要である。
(C37) ldefint(exp(-x)*sin(x),x,0,-minf);
MINF MINF
%E SIN(- MINF) %E COS(- MINF) 1
(D37) - ------------------ - ------------------ + -
2 2 2
(C38) ldefint(exp(-x)*sin(x),x,0,inf);
1
(D38) -
2
(C39) ldefint(exp(-x)*sin(x),x,minf,0);
- x - x
limit - %E SIN(x) - %E COS(x)
x -> MINF 1
(D39) - --------------------------------------- - -
2 2
(C40) ldefint(exp(-x)*sin(x),x,-inf,0);
INF INF
%E SIN(INF) %E COS(INF) 1
(D40) - -------------- + -------------- - -
2 2 2
尚、LDEFINTは内部的には極の判別を一切行わずに、記号積分した結果に上限と下限を
LIMITで代入するだけである。一応、右極限と左極限を下限と上限で取る様にしている
が、単純に上限に対しては'MINUS、下限に対しては'PLUSを内部的に付けるだけなので、
上限や下限で不連続になる関数の場合、上限と下限の大小関係を逆にしてLDEFINTを
計算すれば無意味になる可能性がある。
尚、DEFINTやINTEGRATEで定積分を行う場合は区間内での極限の判別も行っている。
而し、LDEFINTではSININT関数を用いて機械的に記号積分を行うだけである。
その為、関数をいきなりLDEFINTを用いて積分したり、結果の検証を省く事は薦めら
れない。
以下に安易な例を示す。
(C1) ldefint(1/x^2,x,-1,1);
(D1) - 2
(C2) ldefint(1/x^2,x,0,1);
1
(D2) (limit -) - 1
x -> 0 x
(C3) ldefint(1/x^2,x,-1,0);
1
(D3) - (limit -) - 1
x -> 0 x
(C4) d2+d3;
(D4) - 2
(C5) defint(1/x^2,x,-1,1);
Integral is divergent
-- an error. Quitting. To debug this try DEBUGMODE(TRUE);)
(C6) integrate(1/x^2,x,-1,1);
Integral is divergent
-- an error. Quitting. To debug this try DEBUGMODE(TRUE);)
この例で示す様に、D1のLDEFINTの結果は-2となっている。これはMaximaでは1/x^2を
パターンマッチングで安易に記号積分し、区間の上限と下限の極限を取っている為
である。これに対し、DEFINTとINTEGRATEでは極が存在する為にエラーを返している。
この様に、LDEFINTを用いる場合には被積分関数の連続性に関してDEFINTやINTEGRATE
以上に注意が必要となる。
尚、LDEFINTにはZEROA,ZEROBが使える。各々が0の右極限と0の左極限を現わす。
(C10) ldefint(1/x^2,x,zeroa,1);
1
(D10) (limit -) - 1
x -> 0+ x
(C11) ldefint(1/x^2,x,zerob,-1);
1
(D11) (limit -) + 1
x -> 0- x
但し、1+'zeroaの様な使い方は出来ないので注意する。
尚、LDEFINTはDEFINTと同様に積分でRISCH積分を用いる事が出来ない。但し、DEFINTと
同様の修正を加えれば、EVを用いてRISCH積分を用いる事が可能である。
*)
|
POTENTIALZEROLOC[0] |
[indeterminatej=expressionj, indeterminatek=expressionk, ...] |
前者は後者の全ての右手側に対する非リスト式に同値であり、指定された右手側は 定積分の下限として用いられる。積分の成功はそれらの変数と順序に依存する。 POTENTIALZEROLOCはデフォルト値として0が設定されている。
3引数版ではQUANC8('関数名,lo,hi)にて最初の引数で指定された関数の積分を区間lo からhiで計算する。第一引数は'関数名の様に関数名に引用符'を付けたものでなれば ならない。
4引数版ではQUANC8(<f(x)またはxの式>,x,lo,hi)にて、関数か式(最初の引数)の変数 (二番目の引数)で区間loからhiで積分を計算したものとなる。
使われる手法はNewton-Cotesの8次多項式による求積法で、ルーチンは適応型である。 区間の分割に時間を費すのは、大域変数QUANC8_RELERR(デフォルト値=1.0e-8)と QUANC8_ABSERR(デフォルト値=1.0e-8)で指定されたエラー条件に達した時のみであり、 エラー条件は
(* 訳者注： share1/qq.lspの$QUANC8を定義している個所で、DEFUNをDEFMTRFUNに修正する必要が ある。 *) |
(C1) RESIDUE(S/(S**2+A**2),S,A*%I);
1
(D1) -
2
(C2) RESIDUE(SIN(A*X)/X**4,X,0);
3
A
(D2) - --
6
|
(C1) RISCH(X^2*ERF(X),X);
2 2
- X X 3 2
%E (%E SQRT(%PI) X ERF(X) + X + 1)
(D1) ------------------------------------------
3 SQRT(%PI)
(C2) DIFF(%,X),RATSIMP;
2
(D2) X ERF(X)
|
(* 訳者注： INTEGRATEで自動的にRISCHを用いるとあるが、Maxima-5.6でこれが正確に動作してい るか非常に疑わしい。商用のMacsymaでは問題無く動作しているらしいが、Maxima 5.6 ではevで意図的にrischを用いる様にした結果とintegrateのみの結果が異なっている。 integrateでRISCH積分を強制的に利用するにはEVを用いて、 ev(integrate(3^log(x),x),'risch); の様にすると良い。尚、DEFINT等ではこの指定が出来ない為、defint.lispで定義され ている関数ANTIDERIVの修正が必要になる。 *) |
ROMBERG(<被積分関数>,<積分変数>,<下限>,<上限>);
例:
ROMBERG(SIN(Y),Y,1,%PI);
TIME= 39 MSEC. 1.5403023
F(X):=1/(X^5+X+1);
ROMBERG(F(X),X,1.5,0);
TIME= 162 MSEC. - 0.75293843
|
第二のものは効率的なもので、次の様に使われる：
ROMBERG(<関数名>,<下限>,<上限>); |
例：
F(X):=(MODE_DECLARE([FUNCTION(F),X],FLOAT),1/(X^5+X+1));
TRANSLATE(F);
ROMBERG(F,1.5,0);
TIME= 13 MSEC. - 0.75293843
|
最初の引数はTRANSLATE関数で変換された関数か、COMPILE関数で翻訳された関数でな ければならない(COMPILE関数で翻訳されていれば、FLONUMを返す為に宣言されていな ければならない)。最初の引数が未だに変換されたものでなければ、ROMBERGはそれを TRANSLATEで変換せずにエラーを返す。積分の精度は大域変数ROMBERGTOL(デフォルト値 1.E-4)とROMBERGIT(デフォルト値11)で操作される。もし、続く近似で相対誤差が ROMBERGTOLよりも小さければROMBERGは結果を返す。諦める前にROMGERGIT倍の刻幅を 半分にして試みる。ROMBERGが実行する反復と関数評価の数はROMGERGABSとROMBERGMIN で制御される。詳細はDESCRIBE(ROMBERGABS,ROMBERGMIN);を実行せよ。ROMBERGは 再帰的に呼び出されていても良く、それ故、二重、三重積分が実行可能である。
例：
INTEGRATE(INTEGRATE(X*Y/(X+Y),Y,0,X/2),X,1,3);
13/3 (2 LOG(2/3) + 1)
%,NUMER;
0.81930233
DEFINE_VARIABLE(X,0.0,FLOAT,"Global variable in function F")$
F(Y):=(MODE_DECLARE(Y,FLOAT), X*Y/(X+Y) )$
G(X):=ROMBERG('F,0,X/2)$
ROMBERG(G,1,3);
0.8193023
|
この方法の長所は関数Fが他の目的、例えば、プロットの為に使う事が可能な事である。 短所は関数Fとその自由変数Xの両方に対する名前を考慮しなければならない事である。 即ち、大域変数無しであれば：
G1(X):=(MODE_DECLARE(X,FLOAT), ROMBERG(X*Y/(X+Y),Y,0,X/2))$
ROMBERG(G1,1,3);
0.8193023
|
となる。
ここでの長所は簡潔さにある。
Q(A,B):=ROMBERG(ROMBERG(X*Y/(X+Y),Y,0,X/2),X,A,B)$
Q(1,3);
0.8193023
|
この方法はより簡潔なもので、変数はROMBERGの文脈に含まれる為に宣言される必要が 無い。多重積分にROMBERGを利用する事は非常に大きな短所を残念ながら持っている。 多重積分を表現する事で幾何学的情報が欠落する為に膨大な特別な計算が必要とされ る為、この方法は信頼出来ない。利用者はROMBERGTOLとROMBERGITスイッチを正確に 理解して使わなければならない(Romberg積分のIMSL版はMacsymaで利用可能である。 詳細はDESCRIBE(DCADRE);を実行せよ)。
(* 訳者注： rombergの上限と下限に関しては、内部計算でdouble-floatを用いている為、bigfloat で変換した浮動小数は使えない。 *) |
TRUEであれば、0.0B0を返す。それ故、ROMBERGABSが0.0(デフォルト値)であれば相対 誤差検証が得られる。この追加変数は小さな値域で積分計算を行いたい時に便利で ある。そこで、小さな主要な値域で最初に積分する事が相対的精度検証を用いて行え、 後に続く残りの値域上の積分は絶対的精度の検証で用いる。
例： Integral(exp(-x),x,0,50) |
を(数値的に)10000000分の1の相対精度で計算したいとする。Nをカウンターとして、 どれだけの関数評価が必要とされたかを見られる様に次の関数を定義する。
F(X):=(MODE_DECLARE(N,INTEGER,X,FLOAT),N:N+1,EXP(-X))$
TRANSLATE(F)$
/* 事の初めに一度に全体の積分を試みる */
BLOCK([ROMBERGTOL:1.E-6,ROMBERABS:0.],N:0,ROMBERG(F,0,50));
==> 1.00000003
N; ==> 257 /* 関数評価の回数 */
|
それから、先ず最初にIntegral(exp(-x),x,0,10)を実行してROMBERGABSを 1.E-6*(この部分積分)に設定すると云う風に積分を実行する。
BLOCK([ROMBERGTOL:1.E-6,ROMBERGABS:0.,SUM:0.],
N:0,SUM:ROMBERG(F,0,10),ROMBERGABS:SUM*ROMBERGTOL,ROMBERGTOL:0.,
SUM+ROMBERG(F,10,50)); ==> 1.00000001 /* 前の結果と同様 */
N; ==> 130
|
デフォルト値:[11] - ROMBERG積分命令の精度は大域変数のROMBERGTOL[1.E-4]と ROMBERGIT[11]で支配される。ROMBERGは、もし、隣り合った近似解での相対差が ROMBERGTOLよりも小さければ値を返す。諦める前に刻幅のROMGERGITを半分にして試行 する。
デフォルト値:[0] - ROMBERGによる関数評価の最小数を制御する。ROMBERGは その第一引数を少なくとも2^(ROMBERGMIN+2)+1回評価する。これは通常の収束テスト が時々悪い通り方をする時に、周期的関数の積分に対して便利である。
デフォルト値:[1.E-4] - ROMBERG積分命令の精度は大域変数ROMBERGTOL[1.E-4]と ROMBERGIT[11]が支配する。ROMBERGが結果を返すのは、隣り合った近似解の相対差が ROMBERGTOL以下になった時である。諦める前に刻幅のROMBERGITを半分にして試す。
| [ << ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |