第一章 绪 论
复习题
例1计算a(21)6,取21.4,采用下列算式计算: (1)
1(21)6;
(2)99702; (3)(322)3; (4)
1(322)3。
问哪一个得到的结果最好?
解 显然
a(21)6(21)6(21)6(21)61(21)6
(21)6[(21)2]3(322)399702 (21)61(21)61[(21)]231(322)3
所以(1)(2)(3)(4),这4个算式是恒等的,但当取21.4计算时,因为(2),(3)都涉及到两个相近数相减,使有效数字丢失,而(1)在分母算式上的乘幂数比算式(4)大,所以算式(4)最好。事实上,当取21.4时,有x0.015,再由f(x)的误差
|f(xx)f(x)||f'(1.4)||x|也可直接估计出每个算式的误差,显然,算式(4)误
差最小。
具体计算可得: (1)
1(21)65.2103;
.2 . 实用数值分析解题指导
(2)997021.0; (3)(322)38.0103; (4)
1(322)35.1103。
比较可得用第(4)个算式所得的结果更接近于a。
例1.8 建立积分In解 由
10xndxn0,1,,20的递推关系式,并研究它的误差传递。 x51xn5xn11dxxn1dx 005xn1In5In1和
I0可建立下列递推公式
1dxln6ln50x51
I0ln6ln51I5In1nn(n1,2,,20) (*)
计算出I0后,由递推关系式可逐次求出I1,I2,,I20的值。但在计算I0时有舍入误差,因
*此在使用递推关系式中,实际算出的都是近似值In(n1,2,,20)。即
I0*I0e0*1*In5In1n现在来研究误差是如何传递的。
*(n1,2,,20)
设I0有误差e0,假设计算过程中不产生新的舍入误差,则由(*)式可得
**enInIn5In15In(n1,2,) 15en1从而有
en(5)ne0
第一章 绪 论 .3 .
即原始数据I0的误差e0对第n步的影响使该误差扩大了5倍。当n较大时,误差将淹没真值,因此递推公式(*)是数值不稳定的。
现在从另一方向使用这一公式
*n11In1In(n20,19,,1) (**)
5n****只要给出I20的一个近似值I20,即可递推得到I19,类似于上面的推导可得 ,I18,,I0en1每递推一步误差缩小到原值的
由于x[0,1]时,
11en,e0en
55n1,所以递推公式(**)是数值稳定的。 5xnxnxn 65x5所以有估计式
11In
6(n1)5(n1)于是
11I20 621521取
11 I201261050.00873015872I200.0087301587可得另一算法: 11I(I)(n20,19,,1)n1n5n由此可见,对于同一数学问题,使用的算法不同,效率也大不相同,只有选用数值稳定性好的算法,才能求得较准确的结果。
.4 . 实用数值分析解题指导
基于Mathematica的数值计算实例
例1 计算e,有n(n1,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(n10,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!(nk)!例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)9x88x77x66x55x44x33x22x1
在x2处的值,并验证之。
解 Mathematica程序: a[k_]:=k+1; s[8]=a[8];
s[k_]:=x*s[k+1]+a[k]; x=2; s[0]
运行结果:4097。
第二章 插值与拟合
复习题
例2.1 插值函数作为被插函数的逼近,可以用作函数值的近似计算。已知
3273,34,31255,构造二次拉格朗日插值多项式。
.6 . 实用数值分析解题指导
(1)计算3100;
(2)估计误差并与实际误差相比较。 解
(1)以插值点(27,3), (,4), (125,5)代入插值公式,得
2xxj2(x)yili(x)yii0i00xixjjji22 =
(x)(x125)(x27)(x125)(x27)(x)345
(27)(27125)(27)(125)(12527)(125)31002(100) (100)(100125)(10027)(100125)(10027)(100)345(27)(27125)(27)(125)(12527)(125)4.68782
(2) 由误差公式有
f(3)()R(x)(x27)(x)(x125)
3!记f(x)3x,f(3)108(x)x3,f(3)(x)在[27,125]上是单调递减函数。
27f(3)(x)f(3)(27)5.503105
f(3)()R(100)(10027)(100)(100125)0.618131
6实际误差:31002(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 .
解
用线性插值计算,x00.32,x10.34,则
sin0.336 7≈P=y010.3367y1y00.3367x0=0.314 567+
x1x00.3334870.3145670.33670.320.330365
0.340.32其截断误差限为
R1x1sinx\"2\"xxx0xx11M2xx0xx1 2其中M2maxsinxsinx10.3335,于是
x0xx1R10.3367sin0.3367P10.3367
10.3335 20.33670.320.33670.340.92105
/x2x0 0.939350.946
0.04用抛物插值计算,有
y2y1y1y0P20.3367P10.3367xxxx1102 0.3367x00.3367x10.3303650.330374 0.0167×0.0033其截断误差限为
R2xx0xx21M3xx0xx1xx2 3!其中M3maxfxcosx00.950,于是
R20.3367sin0.3367P20.33670.0233<0.20310
610.9500.01670.0033 6=0.330374与具有六位有效数字的正弦函数表完全一样,说明用二次注 P20.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)(xx0)(xx2)(xx0)(xx1)(xx1)(xx2)f(x1)+f(x2) f(x0)+
(xx2)(x0x2)(x1x0)(x1x2)(x2x0)(x2x1)(x3)(x5)(x2)(x5)(x2)(x3)4925
(32)(35)(52)(53)(23)(25)4925(x3)(x5)(x2)(x5)(x2)(x3) 326=
L(x)代入可得L(3.5)12.25。
(2)做出插值点(2, 4)(3, 9)(5, 25)的差商表:
i 0 1 2 xi f[xi] f[xi1,xi] f[xi2,xi1,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](xx0)f[x0,x1,x2](xx0)(xx1)
45(x2)(x2)(x3)
代入可得N(3.5)12.25。
例2.4 设f(x)x2x1.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)(x0)(x1)[x(1)](x1)[x(1)](x0)4.21.20.2
(10)(11)[0(1)](01)[1(1)](10)2 =x2x1.2
以(-1, 4.2),(-0.5, 0.45),(0, 1.2)和(1, 0.2)为插值点,则
L3(x)x22x1.2
以(-1, 4.2),(-0.5, 2.45),(0, 1.2),(0.5, 0.45)和(1, 0.2)为插值点,则
L4(x)x22x1.2
注意到f(3)(x)0,所以n2时Rn(x)0,从而有f(x)Ln(x)(n2)。
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已知函数fx的函数表如下:
求四次牛顿插值多项式,并由此求f0.596的近似值。 分析
表中给出六对数据,故最高可构造五次多项式。但由于0.596接近于
x00.40,因此可取前五对数据来做差商表。
解 构造差商表如下:
.10 . 实用数值分析解题指导
xi 0.40 0.55 0.65 0.80 0.90 fxi 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 故四次牛顿插值多项式为
x0.40.28000x0.4x0.55 P4x0.410751.11600x0.4x0.55x0.650.03134x0.4x0.55 0.19733 x0.65x0.80 于是f0.596≈P40.596=0.631 95。
例2.6
利用差分性质证明:
⑴1+2+3+…+n=
1nn1 2231⑵1+2+3+…+n=nn1
2333分析 本题结论是知道的。但这里要用差分性质来证,那就必须对某个函数构造差分。因此关键要寻找一个函数,使其利用差分的某个性质,经过推导,得出题断。
下面仅对⑴的结论给出两种证法,⑵的结论类似地可证。 证法一
容易看出,当令gn1+2+…+n时,gngn1gn1n,
为此可得差分表如下:
第一章 绪 论 .11 . i 0 1 2 ┇ gi 0 1 1+2 ┇ 1 2 ┇ 2 1 ┇ 1 1 3 0 0 4 n1 gn1 n1 n gn nn 利用差分公式ginnnkkgi,取i=0,则 k0
nknn2gng0g0g0k12g0k0110n1nn11nn1
22即
1+2+3+…+n=
证法二
由于
1nn1 2gn gn1gn1gn1gn2gn2
gn1gn2+gn3g1g1
1nn1,只要证明gi1i(i=1,2,…,n-1)且 211g11即可证得结论。事实上,gi1gigi1ii1i1i
2211ii1i1i2i(i=1,2,…,n-1),g11是显然的。 221故1+2+3+…+n=nn1
2所以定义函数gn例2.7 试对下列数据做出形如abx的拟合曲线。
2.12 . 实用数值分析解题指导
xi yi 2 1 3 6 5 22 7 46 8 61 解 写出(x)abx2的矛盾方程组:
abx12y1abx2y22 2abx5y5法方程为:
1x2112x21x12211x2a122x5bx121x512x2y11y2 2x5y5114619a11111111114929125b492922
4614916155x2ii15yxiai1i1 554b2xixiyii1i152i151a13651517219b6766 解得:a3.00,b1.00。所以有
(x)3.001.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)absinx的拟合函数。
解 矛盾方程为:
absinx1y1absinxy22 absinx5y5法方程组为:
1sinx11sinx21sinx1a11sinx21sinx5bsinx11sinx51sinx2y11y2 sinx5y555sinxii15ysinxiiai1i1 552bsinx(sinx)yiiii1i152.767135a7.660282.767131.662b11.7839
解方程组得:a1.2,b0.6。所以有
(x)1.20.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 试对数据做出形如yae的拟合函数。
解 对函数yae两边取自然对数:
bxlnabxlny
得到数据表:
.14 . 实用数值分析解题指导
xi lnyi 得矛盾方程组:
1.25 3.03815 1.37 3.19793 1.45 3.35575 1.69 3.68.66 1.77 3.83341 lnabx1lny1nabxlny22 nabx5lny5法方程组为:
lnaAAATZ
bT1x11x21x1lna11x21bxx511x51x2lny11lny2
x5lny555xii157.535lnyxiilnai1i1 55bxi2xilnyii1i157.53lna17.1059b26.0501
11.53091.141853.13256,b1.5135。所以有 解方程组可得:lna1.14185,aey(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插值多项式,并计算x1.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.208333x41.75x34.29167x24.75x16
插值为L4(1.5)13.6797。
例3 已知14个点的坐标如下表,根据表中数据构造一插值多项式,并由此计算
x4.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 已知函数ysinx的函数表如下
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.029380.9001x 一次插值多项式
0.018621.1161x0.24x2 二次插值多项式
0.000421.00387x0.0125x20.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.259107
(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.0915t0.0184t20.0142t30.00487t40.000834t50.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(1x)2x(12(1x))x22(1x)x2
化简为: 2x3x22x3
第一章 绪 论 .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求形如y1的拟合函数。 abx解 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]
运行得到结果:a026.2461 ,a124.686。
从而得到拟合曲线方程为
y11
a0a1x26.246124.686x拟合曲线为
例9 设有一形为
第一章 绪 论 .21 .
yy0eax
的函数,现测得函数值表如下:
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]
a01.18014运行结果为:-1.18014 + 3.0968 x。即a13.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]
运行得结果为:
y0ea00.307236,aa13.0968
3.0968x则原函数可近似为:y0.307236e拟合曲线为:
.22 . 实用数值分析解题指导
第三章 线性方程组的解法
复习题
例3.1 用高斯消元法解方程组
2x1x20.1x3x42.7 0.4x0.5x4x8.5x21.91234 0.3xxx5.2x3.9 1234x10.2x22.5x3x49.9 解 取小数位4位对(A,b)做高斯消元:
(A,b)=20.40.3110.510.2.1412.512.78.521.9 5.23.919.9 第一章 绪 论 .23 .
210.112.700.34.028.721.3601.151.0155.054.305
00.32.551.58.55210.112.700.34.028.721.360016.42528.377.575
006.5710.229.91210.112.700.34.028.721.360016.42528.377.575 0001.121.12(A,b)
AXb作回代过程有:
x11.00Xx22.00x
33.x4001.00例3.2 用列主元素法解方程组:
4x1x2x35 18x13x2x315 x1x2x36 解 对(A,b)做选主元及消元过程:
4115(A,b)18311518314111116111155
6由同解方程组.24 . 实用数值分析解题指导
1831151831151757173100
393618617571731003936186180037601151731(A,b)
1862222217由同解方程组AXb作回代过程有:
x11.0000x2.0000 2x3.00003例3.3 分别用Doolittle和Crout分解法求解方程组
2x1x2x34x13x22x36 x2x2x52312解 方程系数矩阵A11(1) Doolittle分解: 对等式
132142,常数项b6。
522A11132112l212l311l32u111u12u22u13u23LU u33按U的第一行、L的第一列,按U的第二行、L的第二列,U的第三行的次序,根据矩阵乘法与矩阵相等,对比等号两边矩阵元素,易得矩阵L,U元素:
第一章 绪 论 .25 .
1L0.50.510.62,U112.511.5 0.6求解LYb,得(y1,y2,y3)(4,4,0.6)。 求解UXY,得(x1,x2,x3)(1,1,1)。
(2) Crout分解
2对等式 A111321l112l212l31l22l321l33u122u13u23LU 3根据矩阵乘法与矩阵相等,对比等号两边矩阵元素,易得矩阵L,U元素:
2L0.512.51.51,U0.60.510.50.6 1求解LYb,得(y1,y2,y3)(2.0,1.6,1.0)。 求解UXY,得(x1,x2,x3)(1.0,1.0,1.0)。
例3.4
用矩阵的直接三角分解法解方程组
10 10分析 解
设
020x15101x23 243x317103x47这是常规计算题,只需按矩阵的三角分解过程认真计算即可。
10100201101l21l24331103l411l32l421l4310u2212u23u330u24 u34u44由矩阵乘法可逐行,逐列分别求出uij和lij:
.26 . 实用数值分析解题指导
1l21 l31l411l32l421l432u23u331
01=121 1010101020
u24101
= u3421u442
10u221y1501y23121y=17
30101y74有y15,y23,y36,y44
再解上三角方程组
1020x15101x23 = 621x32x44得x42,x32,x21,x11。
注
由于三角分解过程十分规律,也可按紧凑格式直接得到
101001212040051130317137001212020011253 64
其中折线右上部的元素(竖线左侧)为上三角阵u的元素,最后一列为y的元素,故
x4 例3.5
6x442, x32, x23x41, x152x31 22用平方根法(Cholesky分解)解方程组
第一章 绪 论 .27 .
323x15 220x2=3
3012x73分析
由于系数矩阵对称正定,故一定有分解形式ALL,其中L为下三角阵,
~~~T~然后由矩阵乘法即可求出L的元素。
解
设
323l11 220=l213012l31由矩阵乘法得 l113,l21由
l22l32l33l11l21l22l31l32 l3323,l313,l2223,l326,l333
l11 l21l31解得y1再由
l22l32l33y15y2=3 y7353,y216,y313。
l11 得x3l21l22l31l32l33x1y1x2=y2 xy3311,x2,x11。 321
例3.11 计算向量X2的1范数、2范数和无穷范数。
1.5
解 由向量范数定义得
.28 . 实用数值分析解题指导
X1=
XK133K=1+2+1.5=4.5
X X=2XK12K=142.25=7.25=2.69258
=MAXXK=2
1K3例3.12 计算矩阵A解 由矩阵范数定义:
1.12.52的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 AFi1j1nnaij21222.523.524.86929
10.957.46AA10.95 16.25
T
ATA特征值λ1=23.61,λ2=0.05591
2 A14.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 用雅可比方法求解方程组
x1x23x33 4x1x2x35 2x5x2x4231方程的准确解为(2,2,1)T。
解 (1) 由系数矩阵
1A42构造雅可比迭代矩阵:
11531 20B41102.5331,g5
20X(k1)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(k1) 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 . 实用数值分析解题指导
λ13.55707,λ21.7782.23374i,λ31.7782.23374i
所以有
(B)1
(2) 将方程的次序调换可以得到
421构造迭代格式:
1511x152x24 3x33(k1)1(k)(k)x(0xx5) 1234(k1)1(k)(k)x2(2x102x34)
5x(k1)1(x(k)x(k)03) 3123写成迭代矩阵形式:
0x1(k1)(k1)2x2x(k1)53131413015(k)4x142(k)5x2 5(k)4x310迭代收敛,取X(0)(1,1,1)T。结果见下表:
k 1 2 3 4 5 6 7 8 9 (k)(k)(k) x1 x1 x1X(k)X(k1) 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.7510 1.99916 -1.99865 1 1.40625103 1.99966 -1.99966 0.999831 1.012510 1.99987 -1.9998 1 2.1093810 433 第一章 绪 论 .31 . 10 11 12 1.99995 -1.99995 0.999975 1.51875104 1.99998 -1.99997 1 3.106105 51.99999 -1.99999 0.999996 2.2781210
例3.17 用高斯—塞德尔方法求解方程组
4x1x2x352x15x22x34 xx3x3231其中方程的准确解为(2,2,1)T。
解 高斯—塞德尔迭代格式:
(k1)1(k)(k)x(0xx5)1234(k1)1(k)(2x1(k1)02x34) x25(k1)1(k1)(x1(k1)x203)x13 写成迭代矩阵形式:
10x1(k1)4(k1)1 x2010x(k1)30120154x1(k)41(k)13x2 2(k)10x13611260,(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(k1) 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 如何对方程组
x18x20x37x10x29x38 9x x x7123进行调整,使得用高斯-塞德尔方法求解时收敛?并取初始向量x求近似解x(k1)(0)(0,0,0)T,用该方法
,使x(k1)x(k)103。
解 将第三个方程调到第一行后有
9x1 x2 x37x18x20x37 x0x9x8231这是主对角线严格占优方程组,故用高斯-塞德尔迭代法求解一定收敛。迭代格式为
x(k1)1(x(k)x(k)7)3192(k1)1(k1)x28(x17) 1(k1)x3(x1(k1)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.0001104103,故取最后结果为
x(4)(1.0000,1.0000,1.0000)T
.34 . 实用数值分析解题指导
例3.19 给定方程组
2111111x111x21
2x31证明:雅可比方法发散,而高斯-塞德尔收敛。 证明 对雅可比方法,迭代矩阵为
0B112设其特征值为,则
12012121 0IB3,10,2,351,故雅可比方法发散。 215i 2(B)对高斯-塞德尔方法,迭代矩阵为
1B1=101011120112001200120121 0121 212112102100012120显然,其特征值为10,2311,(B)1,故高斯-塞德尔收敛。 22T例3.20 设A为正交矩阵,B2IA,求证:线性方程组BBxb用高斯-塞德
尔方法求解必收敛。
T证明:因A为正交矩阵,故AAI,设是A的特征值,则有x0,使
第一章 绪 论 .35 .
Axx
两边与Ax作内积,有
(Ax,Ax)(Ax,x),(x,ATAx)2(x,x)
从而有(x,x)2(x,x),(21)(x,x)0,T又B的特征值为2,故B1,
的特征值不为零,从而B非奇异。故BB正定,所以高斯-塞德尔方法收敛。
例3.21 对方程组x1x2b1
x12x2b2(1) 给出解方程组的雅可比迭代矩阵,并讨论迭代收敛条件; (2) 给出解方程组的高斯-塞德尔迭代矩阵,并讨论迭代收敛条件。 解
(1) 雅可比迭代矩阵:
0 Bρ2其特征多项式为:
ρ 0ρρ22λ λ2λPB(λ)ρ/2故A的谱半径:
(B)即雅可比迭代收敛的充要条件为ρ(2) 高斯-塞德尔迭代矩阵:
2
2。
1SρS的特征多项式为:
021000ρ00ρρ2 2λPS(λ)0ρρ22ρλ(λ) λ22.36 . 实用数值分析解题指导
故S的谱半径:
(S)22
即高斯-塞德尔迭代收敛的充要条件为ρ2。
基于Mathematica的数值计算实例
3x14x25x37x402x3x3x2x01234例1 求解线性方程组
4x11x13x16x023417x12x2x33x40解 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}}。
可得方程组的通解为:
Xk1{13,20,0,17}k2{3,19,17,0}{13k13k2,20k119k2,17k2,17k1}
1x13x22x3x403x5x3x2x01234例2 求解线性方程组
x13x27x35x407x13x25x3x40解 Mathematica程序:
A={{-1,3,2,-1},{3,-5,-3,2},{1,3,7,5},{7,-3,-5,-1}}; MatrixForm[%]
第一章 绪 论 .37 .
Det[A] NullSpace[A]
运行结果为:行列式的值为360,方程组只有零解。
2x1x24x3x486x2x2x2x41234例3 求解线性方程组
5x14x3x42344x13x210x3x48解 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}。
则原方程组的通解为:
1112 Xk1130例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 . 实用数值分析解题指导 运行结果:
630411156 特征多项式
6304111560 特征方程
{116.412, -5.41182, 0, 0, 0, 0} 特征根
450003400123010121001.127960.7023660.2767740.1488170.5744090.1279570.3023660.4767740.6511830.825591100 011 特征向量
2x3yz6例5 用高斯-若当(Gauss-Jordan)消元法解线性方程组x4yz4
x2yz1解 Mathematica程序: A={{2,3,1,6},{1,4,-1,4},{1,-2,1,1}}; MatrixForm[%] RowReduce[A]; MatrixForm[%]
运行结果:
1510085010
8300181553则求得原方程组的解为:x,y,z。
888T例6 设向量x(2,5,1,3,0),求向量范数xp,p1,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 设矩阵
3A122求矩阵A的范数xp20130114 1315,p1,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的性态
11Hn21n解 Mathematica通用程序: Clear[B]
B=Table[1/(i+j),{i,1,k},{j,0,k-1}]; MatrixForm[%] BN=Inverse[B];
12131n11n1n1 12n1.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
7x14x22x318例9 用Jacobi迭代法求解方程组 4x110x2x3406x3x9x29231解 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
2x1x2例10 用Gauss迭代法求解方程组 x31x12x2x31
1x22x31解 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}吻合的十分好。
5x12x2x312例11 用逐次超松弛迭代法求解方程组 x14x22x3202x3x10x3231(取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}。
因篇幅问题不能全部显示,请点此查看更多更全内容
Copyright © 2019- efsc.cn 版权所有 赣ICP备2024042792号-1
违法及侵权请联系:TEL:199 1889 7713 E-MAIL:2724546146@qq.com
本站由北京市万商天勤律师事务所王兴未律师提供法律服务