数值分析作业 - 插值

Problem 1

Lagrange 插值:

Newton 插值:

插值的结果是相同的.

Problem 2

Problem 3

Problem 4

下面的 Mathematica 程序给定插值点和函数值, 计算 Newton 差商矩阵的第一行并借此求出插值函数

newtonCoeffcient[xi_?ListQ, fi_?ListQ] := Module[{a, n, c}, (
    n = Length[xi];
    a = fi;
    c = ConstantArray[0, n];
    c[[1]] = a[[1]];
    Do[(
      Do[
       a[[j]] = (a[[j + 1]] - a[[j]])/(xi[[i + j]] - xi[[j]]), {j, 1,
        n - i}];
      c[[i + 1]] = a[[1]];
      ), {i, 1, n - 1}];
    Return[c];
    )];
newtonValue[xi_?ListQ, a_?ListQ, x_] :=
  Fold[#2[[2]] + (x - #2[[1]]) #1 &, 0, Reverse[Transpose[{xi, a}]]];

先使用 Mathematica 内置的 InterpolatingPolynomial 计算插值多项式

xi = {0, \[Pi]/6, \[Pi]/3, \[Pi]/4, \[Pi]/2};
fi = Sin[xi];
pStd = InterpolatingPolynomial[Transpose[{xi, fi}], x] // Expand;
pStd // TeXForm

再使用 Newton 法计算

1 2	cNewton = newtonCoeffcient[xi, fi] // Simplify; cNewton // TeXForm

1 2	pNewton = newtonValue[xi, cNewton, x] // Expand; pNewton // TeXForm

得到的结果与内置函数完全一致.

将给定点直接代入插值函数求近似值, 在以外的区域显然存在巨大的误差

1 2	x = N[{1, 2, 3, 4, 14, 1000}]; Map[newtonValue[xi, cNewton, #] &, x] // TeXForm

利用三角函数的周期性

pSin[x_] := Piecewise[{
    {newtonValue[xi, cNewton, Mod[x, 2 \[Pi]]],
     Mod[x, 2 \[Pi]] < \[Pi]/2},
    {newtonValue[xi, cNewton, \[Pi] - Mod[x, 2 \[Pi]]],
     Mod[x, 2 \[Pi]] < \[Pi]},
    {-newtonValue[xi, cNewton, Mod[x, 2 \[Pi]] - \[Pi]],
     Mod[x, 2 \[Pi]] < (3 \[Pi])/2},
    {-newtonValue[xi, cNewton, 2 \[Pi] - Mod[x, 2 \[Pi]]], True}
    }];
Map[pSin, x] // TeXForm

得到的结果是合理的

绝对误差

1	Map[pSin, x] - Sin[x] // TeXForm

可以注意到绝对误差均在小数点后三位.

1	Plot[pSin[x] - Sin[x], {x, 0, \[Pi]/2}]

Problem 5

定义一些辅助的变量和函数

uniformSample[a_, b_, n_] := Range[a, b, (b - a)/(n - 1)];
chebshevSample[a_, b_, n_] :=
  Table[(b - a)/2 Cos[((2 i + 1) \[Pi])/(2 (n + 1)) + (b + a)/2], {i,
    0, n}];
f[x_] := 1/(1 + 12 x^2);
newton[f_, xi_, x_] :=
  newtonValue[xi, newtonCoeffcient[xi, f[xi]], x] // Expand;
l = -2.;
r = 2.;

取进行插值并绘制图像

Clear[x];
GraphicsGrid[ParallelMap[{
    Plot[
     {f[x], newton[f, uniformSample[l, r, #], x],
      newton[f, chebshevSample[l, r, #], x]},
     {x, l, r},
     PlotLegends -> Placed[{"f[x]", "uniform", "chebshev"}, Below],
     PlotLabel -> "n=" <> TextString[#]],
    Plot[
     {f[x] - newton[f, uniformSample[l, r, #], x],
      f[x] - newton[f, chebshevSample[l, r, #], x]},
     {x, l, r},
     PlotLegends -> Placed[{"uniform", "chebshev"}, Below],
     PlotLabel -> "Error n=" <> TextString[#]]
    } &, {11, 21, 31}]]

注意到, 在中间的范围内, 两种插值方法得到的函数的形态都接近原函数, 由于均匀方式选取的插值节点更加密集, 在原点附近均匀插值的绝对误差小于 Chebshev 插值. 然而在接近边界区域, 均匀插值方法表现出了明显的震荡现象, 完全不能反映原函数的走向, 而 Chebshev 插值的结果与原函数相比误差较小.

Problem 6

实现三次自然样条插值

splineCoefficient[xi_?ListQ, fi_?ListQ] :=
  Module[{n, a, \[Delta], \[CapitalDelta], m, c, \[Beta], d, b}, (
    n = Length[xi];
    a = Drop[fi, -1];
    \[Delta] = Table[xi[[i + 1]] - xi[[i]], {i, 1, n - 1}];
    \[CapitalDelta] = Table[fi[[i + 1]] - fi[[i]], {i, 1, n - 1}];
    m = SparseArray[Flatten[{
        {1, 1} -> 1, {n, n} -> 1,
        Table[{{i, i + 1} -> \[Delta][[i]], {i, i} ->
           2 \[Delta][[i - 1]] + 2 \[Delta][[i]], {i,
            i - 1} -> \[Delta][[i - 1]]}, {i, 2, n - 1}]
        }], {n, n}];
    \[Beta] =
     SparseArray[
      Table[{i} ->
        3 (\[CapitalDelta][[i]]/\[Delta][[i]] - \[CapitalDelta][[i -
              1]]/\[Delta][[i - 1]]), {i, 2, n - 1}], {n}];
    c = Normal[LinearSolve[m, \[Beta]]];
    d = Table[(c[[i + 1]] - c[[i]])/(
      3 \[Delta][[i]]), {i, 1, n - 1}];
    b = Table[\[CapitalDelta][[i]]/\[Delta][[i]] - \[Delta][[i]]/
        3 (2 c[[i]] + c[[i + 1]]), {i, 1, n - 1}];
    c = Drop[c, -1];
    Return[Transpose[{Drop[xi, -1], a, b, c, d}]];
    )];
splineValue[data_, x_] := Module[{a, b, c, d, xi}, (
    {xi, a, b, c, d} =
     SelectFirst[Reverse[data], x >= #[[1]] &, First[data]];
    Return[a + (b + (c + d (x - xi)) (x - xi)) (x - xi)];
    )];

并绘制图形

GraphicsGrid[ParallelMap[Module[{xi, data},
    xi = uniformSample[l, r, #];
    data = splineCoefficient[xi, f[xi]];
    {Plot[
      {f[x], newton[f, xi, x], splineValue[data, x]},
      {x, l, r},
      PlotLegends -> Placed[{"f[x]", "uniform", "spline"}, Below],
      PlotLabel -> "n=" <> TextString[#]],
     Plot[
      {f[x] - newton[f, xi, x], f[x] - splineValue[data, x]},
      {x, l, r},
      PlotLegends -> Placed[{"uniform", "spline"}, Below],
      PlotLabel -> "Error n=" <> TextString[#]]
     }] &, {11, 21, 31}]]