Romberg integration by Mathematica

mathematica;教學;pdf mathematica 6教學 mathematica基本教學
In[1]: f[x_] :=2/Sqrt[Pi]*Exp[-x^2];
In[2]: R[0, 0] = (b - a)/2*(f[a] + f[b]);
In[3]: R[n_,0]:= R[n,0]=R[n-1,0]/2+(b-a)/2^n*Sum[f[a+(2k-1)*(b-a)/2^n],{k,1,2^(n-1)}];
In[4]: R[n_,m_]:=R[n,m]=R[n,m-1]+(b-a)/(4^m-1)*(R[n,m-1]-R[n-1,m-1]);
In[5]: a=0;b=1;Table[N[R[i,j]],{i,0,4},{j,0,i}]//TableForm


In[6]: err = 10; i = 1;
While[err > 10^(-6),
err = Abs[N[R[i, i] - R[i, i - 1]]];
i++];
Print[N[R[i, i]]]

最後我們可以把上面的程序Module起來，這樣子以後使用會更方便。

In[8]: romberg[fun_,a_,b_]:= Module[{err=10,i=1},
R[0,0]=(b-a)/2*(fun[a]+fun[b]);
R[n_,0]:=R[n,0]=R[n-1,0]/2+(b-a)/2^n*Sum[fun[a+(2k-1)*(b-a)/2^n],
{k,1,2^(n-1)}];
R[n_,m_]:=R[n,m]=R[n,m-1]+(b-a)/(4^m-1)*(R[n,m-1]-R[n-1,m-1]);
While[err>10^(-6),err=Abs[N[R[i,i]-R[i,i-1]]];i++];
N[R[i,i]]]

In[9]: g[x_]:=2/Sqrt[Pi]*Exp[-x^2];
In[10]: romberg[g,0,1]
Out[10]: 0.842701

下面的寫法也是可以


Mathematica 教學 Mathematica 簡易教學 數學軟體 mathematica;教學;pdf mathematica 6教學 mathematica基本教學 mathematica mathematica下載
mathematica 6 mathematica 6.0 mathematica 6下載 mathematica 5.2 mathematica指令

Tags: Teaching

• 編輯文章
• 此文章網址
• 轉寄此文章
• 目前有0 意見