插值与拟合总复习

第二章 插值与拟合 .27 .

?lna?bx1?lny1??na?bx2?lny2 ?????na?bx5?lny5法方程组为：

?lna?TATA???AZ

?b??1??x11x2???1x1???1?1x2?lna??1????????x5?????b??x1??1x5??51x2???lny1???1?lny2?? ???x5????lny5????5?5?x?i??i?1?5??7.53?5??lnyi??xi??lna???i?1i?1? ??5??5??b??xlny??xi2?i???i?i?1??i?1?7.53??lna??17.1059???????

11.5309??b??26.0501?1.14185解方程组可得：lna?1.14185，a?e?3.13256，b?1.5135。所以有

y(x)?3.13256e1.5135x

三、基于Mathematica的数值计算实例

例1 已知f(0)?0,f(1)?3,f{2}?9,求f(x)的二次插值多项式。 解 Mathematica程序： A1 = {{0, 0}, {1, 3}, {2, 9}}; f = InterpolatingPolynomial[A1, x]; Expand[%]

运行结果：

3x2?3x22

例2 已知函数表

. 28 . 实用数值分析解题指导

x y 1 14 2 13 3 10 4 7 5 11 求出Lagrange插值多项式，并计算x?1.5处的y的近似值。

解 Mathematica程序：

r0[x_,x0_,a_,b_,c_,d_]:=(x-a)(x-b)(x-c)(x-d)/((x0-a)(x0-b)(x0-c)(x0-d))

L4[x_]:=r0[x,1,2,3,4,5]*14+r0[x,2,1,3,4,5]*13+r0[x,3,1,2,4,5]*10+r0[x,4,1,2,3,5]*7+r0[x,5,1,2,3,4]*11

L4[x]//N Simplify[%] L3[1.5]//N

运行结果：

插值多项式为

0.583333 (-5. + x) (-4. + x) (-3. + x) (-2. + x) - 2.16667 (-5. + x) (-4. + x) (-3. + x) (-1. + x) + 2.5 (-5. + x) (-4. + x) (-2. + x) (-1. + x) - 1.16667 (-5. + x) (-3. + x) (-2. + x) (-1. + x) +0.458333 (-4. + x) (-3. + x) (-2. + x) (-1. + x)

插值多项式的简式为

0.208333x4?1.75x3?4.29167x2?4.75x?16

插值为L4(1.5)?13.6797。

例3 已知14个点的坐标如下表，根据表中数据构造一插值多项式，并由此计算

x?4.25点的函数值的近似值（不要求给出插值多项式）。

解 已知数据表：

1 2.71828

1.5 4.48169

2.0 7.38906

2.5 12.1825

3.0 20.0855

3.5 33.1155

4.0

4.5

5.0 90.0171

54.5982 54.5982

Mathematica程序：

A=Table[{x,Exp[x]},{x,1,5,0.5}]//N

g1=ListPlot[Table[A],Prolog->AbsolutePointSize[6], DisplayFunction->Identity];

Interpolation[A,InterpolationOrder->1]

g2=Plot[%[x],{x,1,5},PlotStyle->Thickness[0.005],DisplayFunction->Identity] Show[g1,g2,DisplayFunction->$DisplayFunction] %%%[4.25]

运行结果：

第二章 插值与拟合 .29 .

所求插值f(4.25)?70.3076。

例4 已知函数y?sinx的函数表如下

x 0 0 0.1 0.09983 0.2 0.19867 0.3 0.29552 0.4 0.38942 0.5 0.479426 0.6 0.56464 sinx (1) 分别用线性插值、二次插值和三次插值求sin0.57891的近似值，并估计截断误差； (2) 试用4次等距节点插值公式计算f(0.48)的近似值，并估计误差。 解 Mathematica语句：

Table[a = {{0.3,0.4, 0.5, 0.6,}, {0.29552, 0.38942, 0.47943, 0.56464}}]; Transpose[a][[{2, 3}]];

y1 = InterpolatingPolynomial[%, x] // Expand Transpose[a][[{2, 3, 4}]];

y2 = InterpolatingPolynomial[%, x] // Expand Transpose[a][[{1, 2, 3, 4}]];

y3 = InterpolatingPolynomial[%, x] // Expand << Graphics`Legend` Plot[{y1, y2, y3}, {x, 0, 1.8},

PlotStyle -> {Thickness[0.01], Thickness[0.008], Thickness[0.005]}, PlotLegend -> {\Map[NumberForm[#, 10] &, {Sin[0. 48], y1, y2, y3} /. x -> 0.48]

运行结果：

0.02938?0.9001x 一次插值多项式

?0.01862?1.1161x?0.24x2 二次插值多项式

?0.00042?1.00387x?0.0125x2?0.151667x3 三次插值多项式

下面是sin0.48的真值与三种插值的比较：

{0.4617791755, 0.461428, 0.461812, 0.46178288 }

. 30 . 实用数值分析解题指导

最后估计误差的Mathematica程序为： resError[x_, n_, xp_List, k_] := D[Sin[t], {t, n + 1}]/((n + 1)!)*

Product[(x - xp)[[j + 1]] /. t -> k, {j, 0, n}] resError[0.48, 1, {0.4, 0.5}, 0.5] /. t -> 0.5 resError[0.48, 2, {0.4, 0.5, 0.6}, 0.4] /. t -> 0.4 resError[0.48, 3, {0.3, 0.4, 0.5, 0.6}, 0.3] /. t -> 0.3

以下是三种插值在计算sin0.57891时的误差的运行结果

0.00038354 -0.000029474

4.25549?10?7

(2) Mathematica程序： x[k]=k*0.1;

y[0]=0;y[1]=0.09983;y[2]=0.19867;y[3]=0.29552; y[4]=0.38942;y[5]=0.479426;y[6]=0.56464; f1[i_]:=y[i+1]-y[i] f2[i_]:=f1[i+1]-f1[i] f3[i_]:=f2[i+1]-f2[i] f4[i_]:=f3[i+1]-f3[i] f5[i_]:=f4[i+1]-f4[i] f6[i_]:=f5[i+1]-f5[i]

A={{y[0],y[1],y[2],y[3],y[4],y[5],y[6]},{0,f1[0],f1[1],f1[2],f1[3],f1[4],f1[5]}, {0,0,f2[0],f2[1],f2[2],f2[3],f2[4]},{0,0,0,f3[0],f3[1],f3[2],f3[3]}, {0,0,0,0,f4[0],f4[1],f4[2]},{0,0,0,0,0,f5[0],f5[1]},{0,0,0,0,0,0,f6[0]}}; Transpose[A]//N; N[MatrixForm[%],4]

a[0]=y[0];a[1]=f1[0];a[2]=f2[0];a[3]=f3[0];a[4]=f4[0];a[5]=f5[0];a[6]=f6[0]; NN[x]=Sum[a[k]*Product[(t-j),{j,0,k-1}],{k,0,6}]//N N[Expand[%],3]//Chop %/.t->0.48

运行结果： 所得差分表为：

0 0 0 0 0 0 0 0.09983 0.09983 0 0 0 0 0