maximaで定積分を漸化式と台形から解きたい

質問

maxima初心者です。

f(x):=%e^a*x^n*%e^(-x);
a:0.5;
n:20;
を4つの方法で0からaまで定積分し、その誤差を比較し正しい答え(誤差が少ない?)を出す計算式をみつけるという課題が出ました。
integrate(f(x), x, 0, a);
romberg(f(x), x, 0, a);
で出してはみたのですが、漸化式を使った解き方と、台形から解く方法、そして何が正しいのかの見極め方がわかりません。

漸化式は資料に
I[0]:%e^a-1;
for n:1 thru 20 do(
I[n]:n*I[n-1]-a^n);
print(i[20]);
というヒントがあったのですが、正直これに関しては何が何なのかさっぱりわかっていません。

そして台形に関しては
(f(0)+f(0.5))*(0.5)/2;で大きなひとつの台形として計算できると思ったので、
(f(0)+f(0.5)+2*(f(0.1)+f(0.2)+f(0.3)+f(0.4)))
*(0.5/5)/2;
や
(f(0)+f(0.5)+2*(f(0.05)+f(0.1)+f(0.15)+f(0.2)+f(0.25)+f(0.3)+f(0.35)+f(0.4)+f(0.45)))
*(0.5/10)/2;
と分割数を大きくしてみましたが、これだと埒があかないですよね..
どうにもこうにもいかなくなり、質問させていただきました。
せっかくの高性能を全く使いこなせていなくて情けないです。
どなたか超初心者の私にわかりやすく解説していただけないでしょうか。
よろしくお願いします。

回答No.1

kannnami_cho · ベストアンサー

a=1/2,n=20　のとき、式を
　ex:x^20*exp(1/2-x)
とする。そのとき
　ans:integrate(ex,x,0,1/2); の厳密解は

　2432902008176640000*sqrt(%e)-4206024238468833958514521/1048576
　となる。
　しかし得られた答えを一般的方法で小数点表示で表すと、
　%,numer;
521.0
となるがこの答えは間違っている。なぜなら積分のx=0のときの値と
　x=1/2のときの値は精度が16の範囲内では同じであって、17以上では誤差が
　あるからである。
　そのため精度を200,表示数を30とする。
　 fpprec:200;
fpprintprec:30;

ans,bfloat;
2.32340460032264176094970523337b-8

　さて、exをromberg関数で数値積分してみると。
　　romberg(ex,x,0,1/2);　
　2.3234046136998887*10^-8
を得る。
　　romberg関数の精度を15にして数値積分する。
　　rombergtol:1e-15;
　　romberg(ex,x,0,1/2),bfloat;
　　　2.32340460032264175507294877084b-8
　　　を得て精度が向上する。
　　　しかし、精度を16以上にするとromberg関数は収束せず
　　　精度15が限界である。

　台形公式を使う方法はmaximaのヘルプには記載されていないが
　maximaには存在する。
　台形公式を使うためにはsimpsn fileを呼び出す必要がある。
　load(simpsn);
台形関数はtraprule(f,a,b,n) fは関数名、a,bは数値積分の範囲、nは分割数
　関数名が必要なので関数を定義する。

　define(f(x),ex)
あるいは
　f(x):=''ex;
で定義する。

　分割数n=100で計算すると
　　　traprule(f,0,0.5,100);
2.3311490443355299*10^-8
精度はかなり悪い。

　分割数n=10000で計算すると
　　　traprule(f,0,0.5,10000);
2.3234053751829867*10^-8
　精度は向上したがromberg関数に劣る。

　分割数をn=100000で計算すると精度は向上するが
　私のパソコンで約10秒かかる。
　　traprule(f,0,0.5,100000);
　　　2.3234046080712618*10^-8

　以上から、分割数を増大すれば精度が向上するが、計算時間が
　飛躍的に増大すると考えられる。

　漸化式の資料にはａを含むがａを1/2とおく。

　　漸化式は部分積分を使って得られる。
　　つまり、
integrate(exp(1/2-x),x)=-exp(1/2-x)に注意すると
　　　'integrate(x^n*exp(1/2-x),x,0,1/2)=
-(1/2)^n+n*'integrate(x^(n-1)*exp(1/2-x),x,0,1/2);
の等式が成り立つ。
　　　そこで、
　　　　I[n]='integrate(x^n*exp(1/2-x),x,0,1/2)
　　　とすると、次の漸化式が成り立つ。

　　　　I[n]=n*I[n-1]-(1/2)^n

I[0]=integrate(x^0*exp(1/2-x),x,0,1/2);
I[0]=exp(1/2)-1

漸化式を繰り返し使ってn=20まで計算する。
　　　　つまりdo関数を使う。

I[0]:exp(1/2)-1;
　　　for i thru 20 do I[i]:i*I[i-1]-(1/2)^i;

答えは複雑なので小数点表示する。
　　　　I[20],bfloat;
2.32340460032264176094970523337b-8

小数点表示の範囲で厳密解での小数点表示と同じである.
　　　問題はすべて漸化式で表されるとは限らないことである。

お礼

お礼が遅れてしまい申し訳ありません。
丁寧な回答、ありがとうございました。