《计算方法》复习题-1

第一章 绪 论

复习题

例1计算a(21)6，取21.4，采用下列算式计算： （1）

1(21)6；

（2）99702； （3）(322)3； （4）

1(322)3。

问哪一个得到的结果最好？

解 显然

a(21)6(21)6(21)6(21)61(21)6

(21)6[(21)2]3(322)399702 (21)61(21)61[(21)]231(322)3

所以(1)(2)(3)(4)，这4个算式是恒等的，但当取21.4计算时，因为（2），（3）都涉及到两个相近数相减，使有效数字丢失，而（1）在分母算式上的乘幂数比算式（4）大，所以算式（4）最好。事实上，当取21.4时，有x0.015，再由f(x)的误差

|f(xx)f(x)||f'(1.4)||x|也可直接估计出每个算式的误差，显然，算式（4）误

差最小。

具体计算可得： （1）

1(21)65.2103；

．2 ． 实用数值分析解题指导

（2）997021.0； （3）(322)38.0103； （4）

1(322)35.1103。

比较可得用第（4）个算式所得的结果更接近于a。

例1.8 建立积分In解 由

10xndxn0,1,,20的递推关系式，并研究它的误差传递。 x51xn5xn11dxxn1dx 005xn1In5In1和

I0可建立下列递推公式

1dxln6ln50x51

I0ln6ln51I5In1nn(n1,2,,20) (*)

计算出I0后，由递推关系式可逐次求出I1,I2,,I20的值。但在计算I0时有舍入误差，因

*此在使用递推关系式中，实际算出的都是近似值In(n1,2,,20)。即

I0*I0e0*1*In5In1n现在来研究误差是如何传递的。

*(n1,2,,20)

设I0有误差e0，假设计算过程中不产生新的舍入误差，则由(*)式可得

**enInIn5In15In(n1,2,) 15en1从而有

en(5)ne0

第一章 绪 论 ．3 ．

即原始数据I0的误差e0对第n步的影响使该误差扩大了5倍。当n较大时，误差将淹没真值，因此递推公式(*)是数值不稳定的。

现在从另一方向使用这一公式

*n11In1In(n20,19,,1) (**)

5n****只要给出I20的一个近似值I20，即可递推得到I19，类似于上面的推导可得 ,I18,,I0en1每递推一步误差缩小到原值的

由于x[0,1]时，

11en,e0en

55n1，所以递推公式(**)是数值稳定的。 5xnxnxn 65x5所以有估计式

11In

6(n1)5(n1)于是

11I20 621521取

11 I201261050.00873015872I200.0087301587可得另一算法： 11I(I)(n20,19,,1)n1n5n由此可见，对于同一数学问题，使用的算法不同，效率也大不相同，只有选用数值稳定性好的算法，才能求得较准确的结果。

．4 ． 实用数值分析解题指导

基于Mathematica的数值计算实例

例1 计算e,有n(n1,2,...,10)位有效数字的近似值，并列表。 解 Mathematica程序：

Table[{N[E,n],N[Pi,n]},{n,1,10}]; TableForm[%]

运行结果：

3. 3. 2.7 3.1 2.72 3.14 2.718 3.142 2.7183 3.1416 2.71828 3.14159 2.718282 3.141593 2.7182818 3.1415927 2.71828183 3.14159265 2.718281828 3.1415926

例2 用程序计算1500,解 Mathematica程序：

612有n(n10,11,,15)位有效数字的近似值。

Table[{N[Sqrt[1500],n],N[12^(1/6),n]},{n,10,15}]; TableForm[%]

运行结果：

38.72983346 1.513085749 38.729833462 1.5130857494 38.7298334621 1.51308574942 38.72983346207 1.513085749423 38.729833462074 1.5130857494229 38.7298334620742 1.5130857494229

例3 计算cos75.5,arctan34.7,log579的近似值。 解 Mathematica程序：

{Cos[75.5 Degree],ArcTan[34.7 Degree],Log[5,79]}//N; TableForm[%]

运行结果：

0.25038 0.48 2.714

00 第一章 绪 论 ．5 ．

例4 二项式系数定义为C(n,k)解 Mathematica程序： CC[n_,k_]:=n!/(k!*(n-k)!) CC[50,36]

运行结果：

937845656300

n!，利用该定义计算C(50,36)。

k!(nk)!例5 分别给出前20个素数及第100个素数。 解 Mathematica程序： Table[Prime[n], {n, 1, 20}] Prime[100]

运行结果：

前20个素数：

{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71} 第100个素数＝1

例6 用秦九韶法计算

P5(x)9x88x77x66x55x44x33x22x1

在x2处的值，并验证之。

解 Mathematica程序： a[k_]:=k+1; s[8]=a[8];

s[k_]:=x*s[k+1]+a[k]; x=2; s[0]

运行结果：4097。

第二章 插值与拟合

复习题

例2.1 插值函数作为被插函数的逼近，可以用作函数值的近似计算。已知

3273,34,31255，构造二次拉格朗日插值多项式。

．6 ． 实用数值分析解题指导

（1）计算3100；

（2）估计误差并与实际误差相比较。 解

（1）以插值点(27,3), (,4), (125,5)代入插值公式，得

2xxj2(x)yili(x)yii0i00xixjjji22 =

(x)(x125)(x27)(x125)(x27)(x)345

(27)(27125)(27)(125)(12527)(125)31002(100) (100)(100125)(10027)(100125)(10027)(100)345(27)(27125)(27)(125)(12527)(125)4.68782

(2) 由误差公式有

f(3)()R(x)(x27)(x)(x125)

3!记f(x)3x,f(3)108(x)x3,f(3)(x)在[27，125]上是单调递减函数。

27f(3)(x)f(3)(27)5.503105

f(3)()R(100)(10027)(100)(100125)0.618131

6实际误差：31002(100)0.04623。

例2.2已知sin0.32＝0.314 567，sin0.34＝0.333 487，sin0.36＝0.352 274，用线性插值及抛物插值计算sin0.336 7的值并估计截断误差。

分析

题目中相当于告诉了插值条件。考虑到0.336 7位于0.32与0.34之间，根据

插值法的特点，线性插值时，取0.32和0.34作为插值节点；抛物插值时，三个点全取。由于一次、二次插值函数表达式较简单，可采用牛顿型公式，误差估计用拉格朗日型余项表达式。

第一章 绪 论 ．7 ．

解

用线性插值计算，x00.32，x10.34，则

sin0.336 7≈P＝y010.3367y1y00.3367x0＝0.314 567＋

x1x00.3334870.3145670.33670.320.330365

0.340.32其截断误差限为

R1x1sinx\"2\"xxx0xx11M2xx0xx1 2其中M2maxsinxsinx10.3335，于是

x0xx1R10.3367sin0.3367P10.3367

10.3335 20.33670.320.33670.340.92105

/x2x0 0.939350.946

0.04用抛物插值计算，有

y2y1y1y0P20.3367P10.3367xxxx1102 0.3367x00.3367x10.3303650.330374 0.0167×0.0033其截断误差限为

R2xx0xx21M3xx0xx1xx2 3!其中M3maxfxcosx00.950，于是

R20.3367sin0.3367P20.33670.0233＜0.20310

610.9500.01670.0033 6=0.330374与具有六位有效数字的正弦函数表完全一样，说明用二次注 P20.3367插值精度已相当高了。

例2.3 用插值点（2，4），（3，9），（5，25）分别构造拉格朗日插值函数和牛顿插值

．8 ． 实用数值分析解题指导

函数，并计算L(3.5)和N(3.5)。

解

（1）以插值点（2，4），（3，9），（5，25）代入插值公式，得

L(x)(xx0)(xx2)(xx0)(xx1)(xx1)(xx2)f(x1)+f(x2) f(x0)+

(xx2)(x0x2)(x1x0)(x1x2)(x2x0)(x2x1)(x3)(x5)(x2)(x5)(x2)(x3)4925

(32)(35)(52)(53)(23)(25)4925(x3)(x5)(x2)(x5)(x2)(x3) 326=

L(x)代入可得L(3.5)12.25。

（2）做出插值点（2, 4）（3, 9）（5, 25）的差商表：

i 0 1 2 xi f[xi] f[xi1,xi] f[xi2,xi1,xi] 2 3 5 4 9 25 (9-4)/(3-2)=5 (25-9)/(5-2)=8 (8-5)/(5-2)=1

N(x)f[x0]f[x0,x1](xx0)f[x0,x1,x2](xx0)(xx1)

45(x2)(x2)(x3)

代入可得N(3.5)12.25。

例2.4 设f(x)x2x1.2，xi和f(xi)取值如下：

xi 2-1 4.2 -0.5 2.45 0 1.2 0.5 0.45 1 0.2 f(xi) 分别构造L2(x),L3(x),L4(x)，并比较结果。

解 以（-1, 4.2），（0, 1.2）和（1, 0.2）为插值点，则

第一章 绪 论 ．9 ．

L2(x)(x0)(x1)[x(1)](x1)[x(1)](x0)4.21.20.2

(10)(11)[0(1)](01)[1(1)](10)2 ＝x2x1.2

以（-1, 4.2），（-0.5, 0.45），（0, 1.2）和（1, 0.2）为插值点，则

L3(x)x22x1.2

以（-1, 4.2），（-0.5, 2.45），（0, 1.2），（0.5, 0.45）和（1, 0.2）为插值点，则

L4(x)x22x1.2

注意到f(3)(x)0，所以n2时Rn(x)0，从而有f(x)Ln(x)(n2)。

xi yi 0.40 0.410 75 0.55 0.578 15 0.65 0.696 75 0.80 0.888 11 0.90 1.026 52 1.05 1.253 82 例2.5已知函数fx的函数表如下：

求四次牛顿插值多项式，并由此求f0.596的近似值。 分析

表中给出六对数据，故最高可构造五次多项式。但由于0.596接近于

x00.40，因此可取前五对数据来做差商表。

解 构造差商表如下：

．10 ． 实用数值分析解题指导

xi 0.40 0.55 0.65 0.80 0.90 fxi 0.410 75 0.578 15 0.696 75 0.888 11 1.026 52 一阶差商 1.116 00 1.186 00 1.275 73 1.384 10 二阶差商 0.280 00 0.358 93 0.433 48 三阶差商 0.197 33 0.213 00 四阶差商 0.031 34 故四次牛顿插值多项式为

x0.40.28000x0.4x0.55 P4x0.410751.11600x0.4x0.55x0.650.03134x0.4x0.55 0.19733 x0.65x0.80 于是f0.596≈P40.596＝0.631 95。

例2.6

利用差分性质证明：

⑴1＋2＋3＋…＋n＝

1nn1 2231⑵1＋2＋3＋…＋n＝nn1

2333分析 本题结论是知道的。但这里要用差分性质来证，那就必须对某个函数构造差分。因此关键要寻找一个函数，使其利用差分的某个性质，经过推导，得出题断。

下面仅对⑴的结论给出两种证法，⑵的结论类似地可证。 证法一

容易看出，当令gn1＋2＋…＋n时，gngn1gn1n，

为此可得差分表如下：

第一章 绪 论 ．11 ． i 0 1 2 ┇ gi 0 1 1+2 ┇  1 2 ┇ 2 1 ┇ 1 1 3 0 0 4 n1 gn1 n1 n gn nn 利用差分公式ginnnkkgi，取i＝0，则 k0

nknn2gng0g0g0k12g0k0110n1nn11nn1

22即

1＋2＋3＋…＋n＝

证法二

由于

1nn1 2gn gn1gn1gn1gn2gn2

gn1gn2+gn3g1g1

1nn1，只要证明gi1i（i＝1，2，…，n－1）且 211g11即可证得结论。事实上，gi1gigi1ii1i1i

2211ii1i1i2i（i＝1，2，…，n－1），g11是显然的。 221故1＋2＋3＋…＋n＝nn1

2所以定义函数gn例2.7 试对下列数据做出形如abx的拟合曲线。

2．12 ． 实用数值分析解题指导

xi yi 2 1 3 6 5 22 7 46 8 61 解 写出(x)abx2的矛盾方程组：

abx12y1abx2y22 2abx5y5法方程为：

1x2112x21x12211x2a122x5bx121x512x2y11y2 2x5y5114619a11111111114929125b492922

4614916155x2ii15yxiai1i1 554b2xixiyii1i152i151a13651517219b6766 解得：a3.00，b1.00。所以有

(x)3.001.00x2

例2.8 试给出下列数据

xi 0.3 0.5 0.6 0.7 0.9 第一章 绪 论 ．13 ． yi 1.37731 1.48766 1.53879 1.58653 1.67 的形如(x)absinx的拟合函数。

解 矛盾方程为：

absinx1y1absinxy22 absinx5y5法方程组为：

1sinx11sinx21sinx1a11sinx21sinx5bsinx11sinx51sinx2y11y2 sinx5y555sinxii15ysinxiiai1i1 552bsinx(sinx)yiiii1i152.767135a7.660282.767131.662b11.7839

解方程组得：a1.2，b0.6。所以有

(x)1.20.6sinx

例2.9 给出下列数据 xi yi 1.25 20.8666 bx1.37 24.4819 1.45 28.667 1.69 39.6726 1.77 46.22 试对数据做出形如yae的拟合函数。

解 对函数yae两边取自然对数：

bxlnabxlny

得到数据表：

．14 ． 实用数值分析解题指导

xi lnyi 得矛盾方程组：

1.25 3.03815 1.37 3.19793 1.45 3.35575 1.69 3.68.66 1.77 3.83341 lnabx1lny1nabxlny22 nabx5lny5法方程组为：

lnaAAATZ

bT1x11x21x1lna11x21bxx511x51x2lny11lny2

x5lny555xii157.535lnyxiilnai1i1 55bxi2xilnyii1i157.53lna17.1059b26.0501

11.53091.141853.13256，b1.5135。所以有 解方程组可得：lna1.14185，aey(x)3.13256e1.5135x

基于Mathematica的数值计算实例

例1 已知f(0)0,f(1)3,f{2}9,求f(x)的二次插值多项式。

第一章 绪 论 ．15 ．

解 Mathematica程序： A1 = {{0, 0}, {1, 3}, {2, 9}}; f = InterpolatingPolynomial[A1, x]; Expand[%]

3x3x2运行结果： 22例2 已知函数表

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.306

2.5 12.1825

3.0 20.0855

3.5 33.1155

4.0

4.5

5.0 90.0171

.5982 .5982

Mathematica程序：

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

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

．16 ． 实用数值分析解题指导

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]

运行结果：

所求插值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.342 0.5 0.479426 0.6 0.5 sinx (1) 分别用线性插值、二次插值和三次插值求sin0.571的近似值，并估计截断误差； (2) 试用4次等距节点插值公式计算f(0.48)的近似值，并估计误差。 解 Mathematica语句：

Table[a = {{0.3,0.4, 0.5, 0.6,}, {0.29552, 0.342, 0.47943, 0.5}}]; 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 -> {\"y1\Map[NumberForm[#, 10] &, {Sin[0. 48], y1, y2, y3} /. x -> 0.48]

第一章 绪 论 ．17 ． 运行结果：

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 }

最后估计误差的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.571时的误差的运行结果

0.000383 -0.000029474

4.259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.342;y[5]=0.479426;y[6]=0.5; 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

．18 ． 实用数值分析解题指导

N[Expand[%],3]//Chop %/.t->0.48

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

0 0 0 0 0 0 0 0.09983 0.09983 0 0 0 0 0 0.1987 0.09884 -0.00099 0 0 0 0 0.2955 0.09685 -0.00199 -0.001 0 0 0 0.34 0.0939 -0.00295 -0.00096 0.00004 0 0 0.4794 0.09001 -0.0034 -0.000944 0.000016 -0.000024 0 0.56 0.08521 -0.004792 -0.0008 0.000046 0.00003 0.0000 插值多项式为：

0.09983 t - 0.00099 (-1. + t) t - 0.001 (-2. + t) (-1. + t) t + 0.00004 (-3. + t) (-2. + t) (-1. + t) t -0.000024 (-4. + t) (-3. + t) (-2. + t) (-1. + t) t + 0.0000 (-5. + t) (-4. + t) (-3. + t) (-2. + t) (-1. + t) t 化简得： 则

0.0915t0.0184t20.0142t30.00487t40.000834t50.0000t6

f(0.48)0.468457。

例5 过0,1两点构造一个三次插值多项式，满足条件

f(0)0,f'(0)2,f(1)0,f'(1)2。

解 Mathematica程序：

x[0]=0;x[1]=1;y[0]=0;y[1]=1; y'[0]=2;y'[1]=2;

w[x_,1_]:=Product[(x-x[k]),{k,0,1}];

l[x_,k_]:=w[x,1]/((x-x[k])(D[w[x,1],x]/.x->x[k])) Dl[x_,k_]:=D[l[x,k],x]/.x->x[k]

h[x_,k_]:=(1-2(x-x[k])(Dl[x,k]/.x->x[k]))(l[x,k]^2); H[x_,k_]:=(x-x[k])(l[x,k]^2); Table[h[x,k],{k,0,1}]; Table[H[x,k],{k,0,1}];

HH[x_,n_]:=Sum[y[k]*h[x,k]+y'[k]H[x,k],{k,0,n}]; HH[x,1] Expand[%]

运行结果：

2(1x)2x(12(1x))x22(1x)x2

化简为： 2x3x22x3

第一章 绪 论 ．19 ．

例6 求数据表

1

x2 -0.8 0.4295

3 -0.6 0.9826

4 -0.4 1.2392

5 －0.2 2.0000

6 7 0.2 3.0334

4

8 0.

6

3.5601

9

-1 -0.3209

0 2.4445

0.

y3.7787

的最小二乘二次拟合多项式。

解 最小二乘求解的Mathematica程序： Clear[x]

X=Table[-1+0.2*i,{i,0,9}];

Y={-0.3209,0.4295,0.9826,1.2392,2.0000,2.4445,3.0334,3.5601,3.7787,4.15}; q[a_,b_,c_]:=Sum[(a+b*X[[k]]+c*X[[k]]^2-Y[[k]])^2,{k,1,10}] Solve[{D[q[a,b,c],a]==0,D[q[a,b,c],b]==0,D[q[a,b,c],c]==0},{a,b,c}]

运行结果：

拟合多项式系数

{{a -> 2.49261, b -> 2.42838, c -> -0.351241}}

拟合多项式

2.49261 + 2.42838

x + 0.351241x2

LL=Table[{X[[k]],Y[[k]]},{k,1,9}]

g1=ListPlot[LL,Prolog->AbsolutePointSize[10]] f=Fit[LL,{1,x,x^2},x] g2=Plot[f,{x,-1,1}] Show[g1,g2]

运行得拟合效果：

．20 ． 实用数值分析解题指导

例7（双曲拟合）给定数据

x1.0 0.971 .1 0.772 1.2 0.597 1.3 0.429 1.4 0.168 1.5 0.065 1y求形如y1的拟合函数。 abx解 Mathematica程序： Clear[X,Y,f,LL,g1,g2]

L={{1.0,0.971},{1.1,0.772},{1.2,0.597},{1.3,0.429},{1.4,0.168},{1.5,0.065}}; X={1.0,1.1,1.2,1.3,1.4,1.5};

Y={0.971,0.772,0.597,0.429,0.168,0.065}; YY=1/Y;

LL=Table[{X[[n]],YY[[n]]},{n,1,6}] f=Fit[LL,{1,x},x]

g1=ListPlot[L,Prolog->AbsolutePointSize[15]] g2=Plot[1/f,{x,0,3}] Show[g1,g2]

运行得到结果：a026.2461 ，a124.686。

从而得到拟合曲线方程为

y11

a0a1x26.246124.686x拟合曲线为

例9 设有一形为

第一章 绪 论 ．21 ．

yy0eax

的函数，现测得函数值表如下：

x0.1 .2 0.34 .57 0.3 0. 0.4 0.34 0.5 1.49 0.6 1.80 0.7 1.55 0.8 2.48 0.9 3.00 0 y5 试用最小二乘法确定y0和a。

解 Mathematica程序： Clear[A,X,Y,LL,f,g1,g2]

X={0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9}; Y={0.34,0.57,0.,1.34,1.49,1.80,2.55,3.48,5.00}; L=Table[{X[[n]],Y[[n]]},{n,1,9}] LL=Table[{X[[n]],Log[Y[[n]]]},{n,1,9}] f=Fit[LL,{1,x},x]

a01.18014运行结果为：-1.18014 + 3.0968 x。即a13.0968A=-1.18014; b=3.0968; a=Exp[A]; y=a*Exp[b*x]

g1=ListPlot[L,Prolog->AbsolutePointSize[5]] g2=Plot[y,{x,0,0.9}] Show[g1,g2]

运行得结果为：

y0ea00.307236,aa13.0968

3.0968x则原函数可近似为：y0.307236e拟合曲线为：

．22 ． 实用数值分析解题指导

第三章 线性方程组的解法

复习题

例3.1 用高斯消元法解方程组

2x1x20.1x3x42.7 0.4x0.5x4x8.5x21.91234 0.3xxx5.2x3.9 1234x10.2x22.5x3x49.9 解 取小数位4位对(A,b)做高斯消元：

(A,b)=20.40.3110.510.2.1412.512.78.521.9 5.23.919.9 第一章 绪 论 ．23 ．

210.112.700.34.028.721.3601.151.0155.054.305

00.32.551.58.55210.112.700.34.028.721.360016.42528.377.575

006.5710.229.91210.112.700.34.028.721.360016.42528.377.575 0001.121.12(A,b)

AXb作回代过程有：

x11.00Xx22.00x

33.x4001.00例3.2 用列主元素法解方程组：

4x1x2x35 18x13x2x315 x1x2x36 解 对(A,b)做选主元及消元过程：

4115(A,b)18311518314111116111155

6由同解方程组．24 ． 实用数值分析解题指导

1831151831151757173100

393618617571731003936186180037601151731(A,b)

1862222217由同解方程组AXb作回代过程有：

x11.0000x2.0000 2x3.00003例3.3 分别用Doolittle和Crout分解法求解方程组

2x1x2x34x13x22x36 x2x2x52312解 方程系数矩阵A11(1) Doolittle分解： 对等式

132142，常数项b6。

522A11132112l212l311l32u111u12u22u13u23LU u33按U的第一行、L的第一列，按U的第二行、L的第二列，U的第三行的次序，根据矩阵乘法与矩阵相等，对比等号两边矩阵元素，易得矩阵L,U元素：

第一章 绪 论 ．25 ．

1L0.50.510.62，U112.511.5 0.6求解LYb，得(y1,y2,y3)(4,4,0.6)。 求解UXY，得(x1,x2,x3)(1,1,1)。

(2) Crout分解

2对等式 A111321l112l212l31l22l321l33u122u13u23LU 3根据矩阵乘法与矩阵相等，对比等号两边矩阵元素，易得矩阵L,U元素：

2L0.512.51.51，U0.60.510.50.6 1求解LYb，得(y1,y2,y3)(2.0,1.6,1.0)。 求解UXY，得(x1,x2,x3)(1.0,1.0,1.0)。

例3.4

用矩阵的直接三角分解法解方程组

10 10分析 解

设

020x15101x23 243x317103x47这是常规计算题，只需按矩阵的三角分解过程认真计算即可。

10100201101l21l24331103l411l32l421l4310u2212u23u330u24 u34u44由矩阵乘法可逐行，逐列分别求出uij和lij：

．26 ． 实用数值分析解题指导

1l21 l31l411l32l421l432u23u331

01＝121 1010101020

u24101

＝ u3421u442

10u221y1501y23121y＝17

30101y74有y15，y23，y36，y44

再解上三角方程组

1020x15101x23 ＝ 621x32x44得x42，x32，x21，x11。

注

由于三角分解过程十分规律，也可按紧凑格式直接得到

101001212040051130317137001212020011253 64

其中折线右上部的元素（竖线左侧）为上三角阵u的元素，最后一列为y的元素，故

x4 例3.5

6x442， x32， x23x41， x152x31 22用平方根法（Cholesky分解）解方程组

第一章 绪 论 ．27 ．

323x15 220x2=3

3012x73分析

由于系数矩阵对称正定，故一定有分解形式ALL，其中L为下三角阵，

~~~T~然后由矩阵乘法即可求出L的元素。

解

设

323l11 220＝l213012l31由矩阵乘法得 l113，l21由

l22l32l33l11l21l22l31l32 l3323，l313，l2223，l326，l333

l11 l21l31解得y1再由

l22l32l33y15y2＝3 y7353，y216，y313。

l11 得x3l21l22l31l32l33x1y1x2＝y2 xy3311，x2，x11。 321



例3.11 计算向量X2的1范数、2范数和无穷范数。

1.5

解 由向量范数定义得

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

X1=

XK133K=1+2+1.5=4.5

X X=2XK12K=142.25=7.25=2.69258

=MAXXK=2

1K3例3.12 计算矩阵A解 由矩阵范数定义：

1.12.52的1范数、2范数、∞范数和Frobenius范数。 3.5 A1=max{1.1+2.5，2+3.5}=5.5 A∞=max{1.1+2，2.5+3.5}=6 AFi1j1nnaij21222.523.524.86929

10.957.46AA10.95 16.25

T

ATA特征值λ1＝23.61，λ2＝0.05591

2 A14.86355

用Mathematica编写的程序： A={{1.1，-2}，{2.5，-3.5}};

MatrixForm[%];

A1=Max[Sum[Abs[A[[n]]]，{n，1，2}]] B=Transpose[A];

AI=Max[Sum[Abs[B[[m]]]，{m，1，2}]]

AF=Sqrt[Sum[A[[i，j]]^2，{i，1，2}，{j，1，2}]]//N T=Transpose[A].A ; MatrixForm[%]; Eigenvalues[A];

A2=N[Sqrt[Max[Eigenvalues[T]]]，10]

第一章 绪 论 ．29 ． 运行结果为：

5.5 6. 4.86929 4.863706

例3.16 用雅可比方法求解方程组

x1x23x33 4x1x2x35 2x5x2x4231方程的准确解为(2,2,1)T。

解 (1) 由系数矩阵

1A42构造雅可比迭代矩阵：

11531 20B41102.5331,g5

20X(k1)BX(k)g，取X(0)(1,1,1)T。计算结果见下表：

k 1 2 3 4 5 6 7 8 9 (k)(k)(k) x1 x1 x1X(k)X(k1) 1.5 0.5 -8 7.5 9.83333 -9 -4.75 9.5 15.3056 -39.0833 10.6667 30.0833 13.58 -45.5556 80.4028 69.7361 -55.1744 31.1435 98.3241 76.6991 -144.488 324.022 -24.6844 292.878 -178.8 558.268 -667.566 2.881 363.378 52.0271 -1219.02 551.457 1444.93 -2667.53 -495.445 2719.56 可以看到迭代发散。B的特征值。

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

λ13.55707,λ21.7782.23374i,λ31.7782.23374i

所以有

(B)1

(2) 将方程的次序调换可以得到

421构造迭代格式：

1511x152x24 3x33(k1)1(k)(k)x(0xx5) 1234(k1)1(k)(k)x2(2x102x34)

5x(k1)1(x(k)x(k)03) 3123写成迭代矩阵形式：

0x1(k1)(k1)2x2x(k1)53131413015(k)4x142(k)5x2 5(k)4x310迭代收敛，取X(0)(1,1,1)T。结果见下表：

k 1 2 3 4 5 6 7 8 9 (k)(k)(k) x1 x1 x1X(k)X(k1) 1.75 -1.6 1 0.75 1.9 -1.9 0.95 0.3 1.9625 -1.94 1 0.0625 1.985 -1.985 0.9925 0.045 1.99438 -1.991 1 0.009375 1.99775 -1.99775 0.998875 6.7510 1.99916 -1.99865 1 1.40625103 1.99966 -1.99966 0.999831 1.012510 1.99987 -1.9998 1 2.1093810 433 第一章 绪 论 ．31 ． 10 11 12 1.99995 -1.99995 0.999975 1.51875104 1.99998 -1.99997 1 3.106105 51.99999 -1.99999 0.999996 2.2781210

例3.17 用高斯—塞德尔方法求解方程组

4x1x2x352x15x22x34 xx3x3231其中方程的准确解为(2,2,1)T。

解 高斯—塞德尔迭代格式：

(k1)1(k)(k)x(0xx5)1234(k1)1(k)(2x1(k1)02x34) x25(k1)1(k1)(x1(k1)x203)x13 写成迭代矩阵形式：

10x1(k1)4(k1)1 x2010x(k1)30120154x1(k)41(k)13x2 2(k)10x13611260,(S)(B)。所以，高斯—塞德尔迭代格式收敛速度比雅可比因为(S)0.182574迭代快。

取X(0)(1,1,1)为初始值，下表列出迭代计算结果：

T．32 ． 实用数值分析解题指导

k X1(k) 1.25 1.9625 2.01813 2.00457 2.00023 1.999 1.99997 2. (k) X2(k) X3X(k)X(k1) 1 2 3 4 5 6 7 8 -1.7 -2.045 -2.01825 -2.00185 -1.99973 -1.999 -1.99999 -2 1.15 1.0275 1.00004 0.999091 0.999832 0.999999 1.00001 1 2.7 0.7125 0.055625 0.01042 0.00433872 3.43695×10-4 9.966×10-5 2.1×10-6 Mathematica程序：

A={{4，1，-1}，{2，5，2}，{1，1，3}}; MatixForm[%]; b={5，-4，3}; I1=IdentityMatix[3];

DD={{A[[1，1]]，0，0}，{0，A[[2，2]]，0}，{0，0，A[[3，3]]}}; L=-{{0，0，0}，{A[[2，1]]，0，0}，{A[[3，1]]，A[[3，2]]，0}}; DDLN=Inverse[DD-L]; G=I1-DDLN.A; MatrixForm[%];

p=Max[Abs[Eigenvalues[N[G]]]] R=p-1 f=DDLN.b; x[0]={1，1，1}; x[n_]：=G.x[n-1]+f; Table[x[n]，{n，1，8}];

第一章 绪 论 ．33 ． MatrixForm[%]//N Table[Max[Abs[x[n]-x[n-1]]]，{n，1，8}]; MatrixForm[%]//N LinearSolve[A，b]

例3.18 如何对方程组

x18x20x37x10x29x38  9x x x7123进行调整，使得用高斯-塞德尔方法求解时收敛？并取初始向量x求近似解x(k1)(0)(0,0,0)T，用该方法

，使x(k1)x(k)103。

解 将第三个方程调到第一行后有

9x1 x2 x37x18x20x37 x0x9x8231这是主对角线严格占优方程组，故用高斯-塞德尔迭代法求解一定收敛。迭代格式为

x(k1)1(x(k)x(k)7)3192(k1)1(k1)x28(x17) 1(k1)x3(x1(k1)8) 9由x(0)(0,0,0)T，得

x(1)(0.7778,0.9722,0.9753)T x(2)(0.9942,0.9993,0.9994)T x(3)(0.9999,0.9999,0.9999)T x(4)(1.0000,1.0000,1.0000)T

因为 x(4)x(3)0.0001104103，故取最后结果为

x(4)(1.0000,1.0000,1.0000)T

．34 ． 实用数值分析解题指导

例3.19 给定方程组

2111111x111x21

2x31证明：雅可比方法发散，而高斯-塞德尔收敛。 证明 对雅可比方法，迭代矩阵为

0B112设其特征值为，则

12012121 0IB3,10,2,351，故雅可比方法发散。 215i 2(B)对高斯-塞德尔方法，迭代矩阵为

1B1=101011120112001200120121 0121 212112102100012120显然，其特征值为10,2311,(B)1，故高斯-塞德尔收敛。 22T例3.20 设A为正交矩阵，B2IA，求证：线性方程组BBxb用高斯-塞德

尔方法求解必收敛。

T证明：因A为正交矩阵，故AAI，设是A的特征值，则有x0，使

第一章 绪 论 ．35 ．

Axx

两边与Ax作内积，有

(Ax,Ax)(Ax,x),(x,ATAx)2(x,x)

从而有(x,x)2(x,x)，(21)(x,x)0,T又B的特征值为2，故B1，

的特征值不为零，从而B非奇异。故BB正定，所以高斯-塞德尔方法收敛。

例3.21 对方程组x1x2b1

x12x2b2(1) 给出解方程组的雅可比迭代矩阵，并讨论迭代收敛条件； (2) 给出解方程组的高斯-塞德尔迭代矩阵，并讨论迭代收敛条件。 解

(1) 雅可比迭代矩阵：

0 Bρ2其特征多项式为：

ρ 0ρρ22λ λ2λPB(λ)ρ/2故A的谱半径：

(B)即雅可比迭代收敛的充要条件为ρ(2) 高斯-塞德尔迭代矩阵：

2

2。

1SρS的特征多项式为：

021000ρ00ρρ2 2λPS(λ)0ρρ22ρλ(λ) λ22．36 ． 实用数值分析解题指导

故S的谱半径：

(S)22

即高斯-塞德尔迭代收敛的充要条件为ρ2。

基于Mathematica的数值计算实例

3x14x25x37x402x3x3x2x01234例1 求解线性方程组

4x11x13x16x023417x12x2x33x40解 Mathematica程序： Det[A]＝0，

A={{3,4,-5,7},{2,-3,3,-2},{4,11,-13,16},{7,-2,1,3}}; MatrixForm[%]

Det[A] r=Sum[RowReduce[A][[i,i]],{i,1,4}] NullSpace[A]

运行结果为：行列式等于0，系数矩阵的秩为2，有非零解，基础解系为{{-13, -20, 0, 17}, {3, 19, 17, 0}}。

可得方程组的通解为：

Xk1{13,20,0,17}k2{3,19,17,0}{13k13k2,20k119k2,17k2,17k1}

1x13x22x3x403x5x3x2x01234例2 求解线性方程组

x13x27x35x407x13x25x3x40解 Mathematica程序：

A={{-1,3,2,-1},{3,-5,-3,2},{1,3,7,5},{7,-3,-5,-1}}; MatrixForm[%]

第一章 绪 论 ．37 ．

Det[A] NullSpace[A]

运行结果为：行列式的值为360，方程组只有零解。

2x1x24x3x486x2x2x2x41234例3 求解线性方程组 

5x14x3x42344x13x210x3x48解 Mathematica程序： Clear[A,b]

A={{2,1,-4,-1},{6,-2,-2,2},{0,5,14,3},{4,-3,-10,-1}}; MatrixForm[%]; b={8,4,-4,8};

Det[A] r=Sum[RowReduce[A][[i,i]],{i,1,4}] NullSpace[A] LinearSolve[A,b]

运行结果为：行列式等于0，系数矩阵的秩为3，有非零解，基础解系为{{-1, 1, -1, 3}}，一个特解为 {1, 2, -1, 0}。

则原方程组的通解为：

1112 Xk1130例4 方阵A由前36个相邻整数构成，求A的特征多项式、特征方程、特征值与特征向量。

解 Mathematica程序： A=Partition[Range[36], 6]; MatrixForm[%]

CharacteristicPolynomial[A,] CharacteristicPolynomial[A,]＝＝0 Eigenvalues[N[A]]//Chop Eigenvectors[A]//N; MatrixForm[%]

．38 ． 实用数值分析解题指导 运行结果：

630411156 特征多项式

6304111560 特征方程

{116.412, -5.41182, 0, 0, 0, 0} 特征根

450003400123010121001.127960.7023660.2767740.1488170.5744090.1279570.3023660.4767740.6511830.825591100 011 特征向量

2x3yz6例5 用高斯-若当(Gauss-Jordan)消元法解线性方程组x4yz4

x2yz1解 Mathematica程序： A={{2,3,1,6},{1,4,-1,4},{1,-2,1,1}}; MatrixForm[%] RowReduce[A]; MatrixForm[%]

运行结果：

1510085010

8300181553则求得原方程组的解为：x，y，z。

888T例6 设向量x(2,5,1,3,0)，求向量范数xp,p1,2,。

解 Mathematica程序： x={2,-5,1,3,0};

x1=Sum[Abs[x[[i]]],{i,1,5}] 运行结果：11 x2=Sqrt[x.x] 运行结果：39

第一章 绪 论 ．39 ．

x

=Max[Table[Abs[x[[i]]],{i,1,5}]] 运行结果：5

例7 设矩阵

3A122求矩阵A的范数xp20130114 1315,p1,2,。

解 Mathematica程序：

A={{3,2,0,1},{-1,0,-1,4},{2,1,1,3},{2,3,1,5}}; MatrixForm[%];

A1=Max[Sum[Abs[A[[n]]],{n,1,4}]]

A100=Max[Table[Sum[Abs[A[[n,m]]],{m,1,4}],{n,1,4}]] T=Transpose[A].A; MatrixForm[%]; Eigenvalues[N[T]];

A2=N[Sqrt[Max[Eigenvalues[N[T]]]]]

运行结果： 13 11 8.253

例8 考察n阶希尔伯特(Hilbert)矩阵Hn的性态

11Hn21n解 Mathematica通用程序： Clear[B]

B=Table[1/(i+j),{i,1,k},{j,0,k-1}]; MatrixForm[%] BN=Inverse[B];

12131n11n1n1 12n1．40 ． 实用数值分析解题指导

B1=Max[Sum[Abs[B[[n]]],{n,1,k}]]; BN1=Max[Sum[Abs[BN[[n]]],{n,1,k}]]; condA1=B1*BN1//N

B100=Max[Table[Sum[Abs[B[[n,m]]],{m,1,k}],{n,1,k}]]; BN100=Max[Table[Sum[Abs[BN[[n,m]]],{m,1,k}],{n,1,k}]]; condA100=B100*BN100//N

当k取4时, 运行结果： 28375 28375

7x14x22x318例9 用Jacobi迭代法求解方程组 4x110x2x3406x3x9x29231解 Mathematica程序：

A = {{7, -4, 2}, {4, 10, -1}, {6, 3, 9}}; MatrixForm[%]; b = {18, 40, 29}; I1 = IdentityMatrix[3];

DD = {{A[[1, 1]], 0, 0}, {0, A[[2, 2]], 0}, {0, 0, A[[3, 3]]}} DDN = Inverse[DD]; J = I1 - DDN.A; MatrixForm[%];

R = Max[Abs[Eigenvalues[N[J]]]] f = DDN.b; xa = {0, 0, 0};

Do[xb = J.xa + f // N; z = Max[Abs[xb - xa]] // N;

Print[k , \" \" , xb , \" \" , z]; xa = xb, {k, 1, 10}]

运行结果：0.4927

即Jacobi迭代法迭代矩阵的谱半径小于1，迭代法收敛。

n 近似解 误差

1 {2.57143, 4., 3.22222} 4. 2 {3.93651, 3.29365, 0.174603} 3.04762 3 {4.40363, 2.44286, -0.5} 0.850794 4 {4.1102, 2.18855, -0.527816} 0.293424 5 {3.97283, 2.30314, -0.24743} 0.280385 6 {3.9582, 2.38612, -0.194045} 0.0829873 7 {3.99037, 2.39732, -0.211953} 0.0321684 8 {4.00188, 2.38266, -0.237129} 0.0251759

第一章 绪 论 ．41 ．

9 {4.0007, 2.37553, -0.239917} 0.0071222 10 {3.99742, 2.37573, -0.236755} 0.00327318 11 {3.99663, 2.37735, -0.234637} 0.00211733 12 {3.99696, 2.37788, -0.234651} 0.00052876 13 {3.99726, 2.37775, -0.235043} 0.000392207 14 {3.9973, 2.37759, -0.235203} 0.000161619 15 {3.99725, 2.37756, -0.235174} 0.0000465373

2x1x2例10 用Gauss迭代法求解方程组 x31x12x2x31

1x22x31解 Mathematica程序：

A={{2,1,1},{1,2,1},{0,1,2}}; MatrixForm[%] b={-1,-1,-1};

I1=IdentityMatrix[3];;

DD={{A[[1,1]],0,0},{0,A[[2,2]],0},{0,0,A[[3,3]]}}; L=-{{0,0,0},{A[[2,1]],0,0},{A[[3,1]],A[[3,2]],0}}; DDLN=Inverse[DD-L]; G=I1-DDLN.A; MatrixForm[%];

R=Max[Abs[Eigenvalues[N[G]]]] f=DDLN.b; xa={0,0,0};

Do[xb=G.xa+f//N;z=Max[Abs[xb-xa]]//N;

Print[k , \" \" , xb , \" \" , z];xa=xb,{k,1,15}] LinearSolve[A,b] 运行结果：0.375

即迭代矩阵的谱半径小于1，迭代法收敛。

n 近似解 误差

1 {-0.5, -0.25, -0.375} 0.5 2 {-0.1875, -0.21875, -0.390625} 0.3125 3 {-0.195312, -0.207031, -0.3984} 0.0117188 4 {-0.198242, -0.202637, -0.398682} 0.00439453 5 {-0.199341, -0.2009, -0.399506} 0.001795 6 {-0.199753, -0.200371, -0.399815} 0.000617981 7 {-0.199907, -0.200139, -0.39993} 0.000231743 8 {-0.199965, -0.200052, -0.399974} 0.0000869036 9 {-0.199987, -0.20002, -0.39999} 0.0000325888

；．42 ． 实用数值分析解题指导

10 {-0.199995, -0.200007, -0.399996} 0.0000122208 11 {-0.199998, -0.200003, -0.399999} 4.58281×10-6 12 {-0.199999, -0.200001, -0.399999} 1.71855×10-6 13 {-0.2, -0.2, -0.4} 6.44457×10-7 14 {-0.2, -0.2, -0.4} 2.41671×10-7 15 {-0.2, -0.2, -0.4} 9.06268 ×10-8 {-0.2, -0.2, -0.4}

由此可见，迭代结果与真解{02,0.2,0.4}吻合的十分好。

5x12x2x312例11 用逐次超松弛迭代法求解方程组 x14x22x3202x3x10x3231（取1.25,x0{1,1,1}T）。 解 Mathematica语句： A={{4,3,0},{3,4,-1},{0,-1,4}}; MatrixForm[%]; Det[A];

b={24,30,-24};w=1.25; I1=IdentityMatrix[3];;

DD={{A[[1,1]],0,0},{0,A[[2,2]],0},{0,0,A[[3,3]]}}; L=-{{0,0,0},{A[[2,1]],0,0},{A[[3,1]],A[[3,2]],0}}; DDLN=Inverse[DD-w*L]; SOR=I1-w*DDLN.A; MatrixForm[%];

R=Max[Abs[Eigenvalues[N[SOR]]]] g=w*DDLN.b; f[x_]:=SOR.x+g

FixedPointList[f,{1.,1.,1.}]; TableForm[%];

Table[FixedPoint[f,{1.,1.,1.},n],{n,1,15}]; TableForm[%]

运行结果：0.25

即迭代矩阵的谱半径=0.32625小于1，迭代法收敛。

运行得到的迭代解为：

第一章 绪 论 ．43 ．

6.3125 3.51953 -6.65015 2.62231 3.95853 -4.60042 3.1333 4.01026 -5.09669 2.95705 4.00748 -4.97349 3.00372 4.00292 -5.00571 2.99633 4.00093 -4.99828 3.00005 4.00026 -5.00035 2.99975 4.00007 -4.999 3. 4.00001 -5.00002 2.99999 4. -4.99999 3. 4. -5. 3. 4. -5. 3. 4. -5. 3. 4. -5. 3. 4. -5.

该方程组的精确解为：{3,4,-5}。