数值分析作业 - 数值积分

数值分析

Problem 1

对一次多项式

二次多项式

则代数精度是 1.

Problem 2

使用待定系数法:

Problem 3

则需要确定二次多项式使其在与任意次数少于一的多项式正交.

则的零点为

则二点 Gauss 积分直接得到了精确值.

Problem 4

定义四种积分算法

centerIntergral[f_, a_, b_, m_] := Module[{\[CapitalDelta]h},
   \[CapitalDelta]h = (b - a)/m;
   Return[\[CapitalDelta]h Sum[
      f[a + \[CapitalDelta]h/2 + i \[CapitalDelta]h], {i, 0,
       m - 1}]];
   ];
trapezoidIntergral[f_, a_, b_, m_] := Module[{\[CapitalDelta]h},
   \[CapitalDelta]h = (b - a)/m;
   Return[\[CapitalDelta]h/
     2 (f[a] + f[b] +
       2 Sum[f[a + i \[CapitalDelta]h], {i, 1, m - 1}])];
   ];
simpsonIntergral[f_, a_, b_, m_] := Module[{\[CapitalDelta]h},
   \[CapitalDelta]h = (b - a)/m;
   Return[\[CapitalDelta]h/
     6 (f[a] + f[b] +
       2 Sum[f[a + i \[CapitalDelta]h], {i, 1, m - 1}] +
       4 Sum[f[a + \[CapitalDelta]h/2 + i \[CapitalDelta]h], {i, 0,
          m - 1}])];
   ];
gauss3IntegralSegment[f_, a_, b_] := (b - a)/
   2 (5/9 f[(a + b)/2 - (a - b)/2 Sqrt[3/5]] + 8/9 f[(a + b)/2] +
     5/9 f[(a + b)/2 + (a - b)/2 Sqrt[3/5]]);
gauss3Integral[f_, a_, b_, m_] :=
  Sum[gauss3IntegralSegment[f, a + (b - a)/m i,
    a + (b - a)/m (i + 1)], {i, 0, m - 1}];

正确积分值是 , 保留 20 位有效数字的工作精度.

f[x_] = N[x E^x, 20];
a = 0;
b = 1;
m = {1, 2, 4, 8, 16, 32};
real = Integrate[f[x], {x, a, b}]

用四种方法和不同积分区间数量求解积分

TableForm[
  Transpose[
   Map[Table[
      ScientificForm[1 - #[f, a, b, 2^i], 5], {i, 0,
       5}] &, {centerIntergral, trapezoidIntergral, simpsonIntergral,
     gauss3Integral}]],
  TableHeadings -> {Table[2^i, {i, 0, 5}], {"center", "trapezoid",
     "simpson", "gauss"}}] // TeXForm

能注意到不同方法的误差都随着区间划分增多而减少, 相同积分区间数量时, 四种方法的误差依次减小.

Problem 5

romberg[f_, a_, b_, n_] := Module[{\[CapitalDelta]h, r},
   \[CapitalDelta]h = b - a;
   r = {{\[CapitalDelta]h/2 (f[a] + f[b])}};
   Do[
    \[CapitalDelta]h /= 2;
    AppendTo[r, ConstantArray[0, j]];
    r[[j, 1]] =
     1/2 r[[j - 1, 1]] + \[CapitalDelta]h Sum[
        f[a + (2 i - 1) \[CapitalDelta]h], {i, 1, 2^(j - 2)}];
    Do[
     r[[j, k]] = (4^(k - 1) r[[j, k - 1]] - r[[j - 1, k - 1]])/(
       4^(k - 1) - 1);
     , {k, 2, j}]
    , {j, 2, n}];
   Return[r];
   ];