《計算方法》復(fù)習(xí)題

#/43en=(-5)e°即原始數(shù)據(jù)I。*的誤差eo對第n步的影響使該誤差擴大了 5n倍。當(dāng)n較大時，誤差將淹沒真值，因此遞推公式(*)是數(shù)值不穩(wěn)定的。 文檔來自于網(wǎng)絡(luò)搜索現(xiàn)在從另一方向使用這一公式TOC\o"1-5"\h\z1(1 )In」 In(n=20,19, ,1) (**)55 丿只要給出丨20的一個近似值丨只要給出丨20的一個近似值丨；0，即可遞推得到1；9，丨；8,…，IO，類似于上面的推導(dǎo)可得每遞推一步誤差縮小到原值的由于x?每遞推一步誤差縮小到原值的由于x?[0,1]時，1丄，所以遞推公式(**)是數(shù)值穩(wěn)定的。5n n nxxx一< <一6 5x5所以有估計式6(n1).16(n1).15(n1)于是11

￡I20￡621 521丄?丄I20:126105 0.00873015872I20 0.0087301587可得另一算法： 11|In4=：(—Tn)(n=20,19,…，1)L5n由此可見，對于同一數(shù)學(xué)問題，使用的算法不同，效率也大不相同，只有選用數(shù)值穩(wěn)定性好的算法，才能求得較準(zhǔn)確的結(jié)果。 文檔來自于網(wǎng)絡(luò)搜索

基于Mathematica的數(shù)值計算實例例1計算e,二有n(n=1,2,...,10)位有效數(shù)字的近似值，并列表。解Mathematica程序：Table[{N[E,n],N[Pi,n]},{n,1,10}];TableForm[%]運行結(jié)果：3.3.2.73.12.723.142.7183.1422.71833.14162.718283.141592.7182823.1415932.71828183.14159272.718281833.141592652.7182818283.141592654例2用程序計算.1500,612有n(n=10,11，…，15)位有效數(shù)字的近似值。解Mathematica程序:Table[{N[Sqrt[1500],n],N[12A(1/6),n]},{n,10,15}];TableForm[%]運行結(jié)果：1.5130857491.51308574941.5130857491.51308574941.513085749421.5130857494231.51308574942291.513085749422938.72983346238.729833462138.7298334620738.72983346207438.7298334620742例3計算cos75.50,arctan34.7°,log579的近似值。解Mathematica程序：{Cos[75.5Degree],ArcTan[34.7Degree],Log[5,79]}〃N;文檔來自于網(wǎng)絡(luò)搜索TableForm[%]運行結(jié)果：0.250380.5445482.71489例4二項式系數(shù)定義為C(n,k) n! ，利用該定義計算C(50,36)。k!(n—k)!解Mathematica程序：CC[n_,k」:=n!/(k!*(n-k)!)CC[50,36]運行結(jié)果：937845656300例5分別給出前20個素數(shù)及第100個素數(shù)。解Mathematica程序：Table[Prime[n],{n,1,20}]Prime[100]運行結(jié)果：前20個素數(shù)：{2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71}文檔來自于網(wǎng)絡(luò)搜索第100個素數(shù)=541例6用秦九韶法計算P5(x)=9x88x77x66x55x44x33x22x1在x=2處的值，并驗證之。解Mathematica程序：a[k_]:=k+1;s[8]=a[8];s[k_]:=x*s[k+1]+a[k];x=2;s[0]運行結(jié)果：4097。第二章 插值與擬合復(fù)習(xí)題例2.1插值函數(shù)作為被插函數(shù)的逼近，可以用作函數(shù)值的近似計算。已知327=3,364=4,3,125=5，構(gòu)造二次拉格朗日插值多項式。

(1)計算3100；(2)估計誤差并與實際誤差相比較。解(1)以插值點(27,3),(64,4),(125,5)代入插值公式，得x_Xjx_Xjn j衛(wèi)Xi-Xj「2(x)-'Jyili(x)八yiiA i(100-64)(100-125)(27-64)(27-125)=(X-64)(x-125)述3+(x-27)(x-125)如十(x-27)(x-64)=(100-64)(100-125)(27-64)(27-125).(100-27)(100-125)4(100-27)(100-64)5(64-27)(64-125) (125-27)(125-64)=4.68782(2)由誤差公式有R(x)=S(x-27)(x-64)(x-125)3!記f(x)=記f(x)=3"(x)尋*(x)在［27，125］上是單調(diào)遞減函數(shù)。■R(100)<f■R(100)<f⑶(x)豈f(3)(27) 5.6450310*^(100-27)(100-64)(100-125)“6181316實際誤差:^/實際誤差:^/103―駕(100)=0.04623。例2.2已知sin0.32=0.314567,sin0.34=0.333487,sin0.36=0.352274，用線性插值文檔來自于網(wǎng)絡(luò)搜索及拋物插值計算sin0.3367的值并估計截斷誤差。文檔來自于網(wǎng)絡(luò)搜索分析題目中相當(dāng)于告訴了插值條件。考慮到0.3367位于0.32與0.34之間，根據(jù)插值法的特點，線性插值時，取 0.32和0.34作為插值節(jié)點；拋物插值時，三個點全取。由于一次、二次插值函數(shù)表達(dá)式較簡單，可采用牛頓型公式，誤差估計用拉格朗日型余項表達(dá)式。文檔來自于網(wǎng)絡(luò)搜索

用線性插值計算，xo=0.32,xi=0.34，則y.-y0sin0.3367?R0.3367=y- -0.3367-x°=0.314567+xi—x00.3367-0.32=0.3303650.333487-0.3145670.3367-0.32=0.3303650.34—0.32其截斷誤差限為1"1"R(x)=2(sinx)|x=g(x-x0x-xi冷M2x-x0x_xi其中m2=max(sinx)=sinx^0.3335，于1R(0.3367)=sin0.3367-R(0.3367廬?x0.3335x0.3367-0.32 0.3367-0.34<0.9210」用拋物插值計算，有P2(0.3367)=R(0.3367)+||y2_y1_y1_y0)/(x2_x0)漢[\X2—X1 X1—X。丿0.3367-x00.3367-x00.3367-x11=0.3303650.93935-0.9460.040.0167x-0.0033=0.330374其截斷誤差限為1R2X乞弓M3X-x°x-X1x-X2其中M3=max廠"(xD=cosx0v0.950，于是X0空魚1R20.3367I；|sinO.3367-P20.3367 0.9500.01670.003360.0233V0.20310注 P20.3367=0.330374與具有六位有效數(shù)字的正弦函數(shù)表完全一樣，說明用二次插值精度已相當(dāng)高了。例2.3用插值點(2,4),(3,9),(5,25)分別構(gòu)造拉格朗日插值函數(shù)和牛頓插值

函數(shù)，并計算L（3.5）和N（3.5）。文檔來自于網(wǎng)絡(luò)搜索解（1）以插值點（2,4）,（3,9）,（5,25）代入插值公式，得(X—X1(X—X1)(X-X2)f(xo)+(x-sx-xj(x-X2)(Xo-X2) (Xi-Xo)(Xi-X2)f(xi)+(X-Xo)(X-Xi)

(X2-Xo)(X2-Xi)f(X2)(X-3)(x-5) 4(x-2)(x-5) 9(x-2)(x-3)(2-3)(2-5) (3-2)(3-5) (5-2)(5-3)9 25l(x)=3(x_3)(x-5)-2(x-2)(x-5) $(x-2)(x-3)代入可得L(3.5)=12.25。（2）做出插值點（2,4）（3,9）（5,25）的差商表:iXif[Xi]f[Xi4，Xi]f[Xi_2，Xi』Xi]024i39(9-4)/(3-2)=52525(25_9)/(5_2)=8(8-5)/(5-2)=1N(x)二f[Xo] f[Xo,Xi](X-Xo)f[Xo,Xi,X2】(X-Xo)(X-Xi)=4 5(x-2)(x-2)(x-3)代入可得N(3.5)h2.25。例2.4設(shè)f(x)=x2-2xi.2，xi和f(xi)取值如下:Xi-1-0.500.51f(Xi)4.22.451.20.450.2分別構(gòu)造L2（x）,L3（x）,L4（x），并比較結(jié)果。解以（-i,4.2），（0,i.2）和（i,0.2）為插值點，則L2（x）=（x—°）（x「）如2+肢一（一1）］「1）“2+以一（一1）］?一叭02（_1_0）（-1-1） ［0-（-1）］0-1）［1-（-1）］1-0）=X2-2x1.2以(-1,4.2),(-0.5,0.45),(0,1.2)和(1,0.2)為插值點，則L3(x)=x2-2x1.2以(-1,4.2),(-0.5,2.45),(0,1.2),(0.5,0.45)和(1,0.2)為插值點，貝U文檔來自于網(wǎng)絡(luò)搜索L4(x)=x2-2x1.2注意到f⑶（x）=0,所以n_2時Rn（x）=0，從而有f（x）wLn（x）（n一2）。Xi0.400.550.650.800.901.05yi0.410750.578150.696750.888111.026521.25382例2.5已知函數(shù)fx的函數(shù)表如下：求四次牛頓插值多項式，并由此求 f0.596的近似值。分析 表中給出六對數(shù)據(jù)，故最高可構(gòu)造五次多項式。但由于 0.596接近于X0=0.40，因此可取前五對數(shù)據(jù)來做差商表。解構(gòu)造差商表如下：Xif(Xi)一階差商二階差商三階差商四階差商0.400.410750.550.578151.116000.650.696751.186000.280000.800.888111.275730.358930.197330.901.026521.384100.433480.213000.03134故四次牛頓插值多項式為P4x1=0.410751.11600X-0.4 0.28000x-0.4x-0.550.19733x-0.4x-0.55x-0.65 0.03134x-0.4x-0.55x-0.65x-0.80于是f0.596~P40.596=0.63195。例2.6利用差分性質(zhì)證明：1’⑴1+2+3+…+n=nn，12⑵13+23+33+…+n3=」nn112」分析本題結(jié)論是知道的。但這里要用差分性質(zhì)來證，那就必須對某個函數(shù)構(gòu)造差分。因此關(guān)鍵要尋找一個函數(shù)，使其利用差分的某個性質(zhì)，經(jīng)過推導(dǎo)，得出題斷。 文檔來自于網(wǎng)絡(luò)搜索下面僅對⑴的結(jié)論給出兩種證法，⑵的結(jié)論類似地可證。證法一 容易看出，當(dāng)令gni>1+2+-+n時，gn-gn-1=.ign斗二n,為此可得差分表如下:ig(i)A心2△400111+2121111211111n-1g(n-1)11n-11110ng(n)n10ikgi，取i=0，則nn利用差分公式gjn='、TOC\o"1-5"\h\zn'n'k 'n" G2g(n)二瓦廠$g(o)=g(o)+d勾(0)+°Q2g(0)

心卞丿 J丿 i2丿4 」 1 」=0n1nn-11nn12即1「1+2+3+…+n=nn12證法二 由于TOC\o"1-5"\h\zgn= ：gn—1gn—1=gn—1gn—2gn—2= =：gn—1；=gn—2+：gn—3； =g1g11所以定義函數(shù)gn=^nn，1，只要證明■:gi-11=i(i=1,2,…，n—1)且j “ j j 1j 1jg1=1即可證得結(jié)論。事實上， gigi-gi-1=?ii?1i-1i二\o"CurrentDocument"1 1—i(i+1—i+1)=—i2=i(i=1,2,…，n—1),g(1)=1是顯然的。\o"CurrentDocument"2 21故1+2+3+-+n=nn?12例2.7試對下列數(shù)據(jù)做出形如a'bx2的擬合曲線。

Xi23578yi16224661a+bx；.?2a+bx2解寫出“x)a+bx；.?2a+bx2.2abx5法方程為:~112X22X12X2■y/l2549642X5~112X22X12X2■y/l2549642X52549亦1MH_412X22X5254911642246641517219」〉1517219」〉[6766一解得：a解得：a=—3.00,b=1.00。所以有5P 21ZXi工yii#iaVi#5送x25送X；5送Xi2yii# -1J=1555151a136(x)n-3.001.00X2例2.8試給出下列數(shù)據(jù)Xi0.30.50.60.70.9

yi1.377311.487661.538791.586531.67的形如「（x）=a-bsinx的擬合函數(shù)。解矛盾方程為:a+bsin%=y1a+bsinX2=y2absinx5二法方程組為:'I11sin%sinx2111'I11sin%sinx21111sinx5sinxsinx2a[ -1?丨I- ]b_ sinx1sinx51sinx211y2sinx5ay一■y/l55Zsinxi5sinXji55Zsinxi5sinXji=152、sin人i=12yiy5E(sinxjyi22.76713a_7.66028[2.767131.66462」一[11.7839_解方程組得：a=—1.2,b=0.6。所以有(xH-1.20.6sinx例2.9給出下列數(shù)據(jù)Xi1.251.371.451.691.77yi20.866624.481928.66739.672646.22bx試對數(shù)據(jù)做出形如y=ae的擬合函數(shù)。解對函數(shù)y=aebx兩邊取自然對數(shù):lnabx=Iny得到數(shù)據(jù)表:

Xi1.251.371.451.691.77Inyi3.038153.197933.355753.68.663.83341得矛盾方程組:Ina+b^=Iny1

na+bx2=Iny2

I…nabx5二In法方程組為:ATATA?nal二Az"Inal乞In"Inal乞Inyi#[b-5送XiInyi1255''Xji=157Xi2i￡"11…111 x2Ina]_[1 1 …1]Iny2,1 X2…X5i aJb ,1X2…X5aJX5一_lny5一~1 xj _lnyi|-5Zxii=15 7.53Ina_17.10597.53 11.5309|tb, [26.0501解方程組可得:Ina=1.14185,a二e1.14185=3.13256,b=1.5135。所以有解方程組可得:y(x^3.13256e1.5135x基于Mathematica的數(shù)值計算實例例1已知f(0)=0,f(1)=3,f{2}=9,求f(x)的二次插值多項式。解Mathematica程序:A1={{0,0},{1,3},{2,9}};f=lnterpolatingPolynomial[A1,x];Expand]%]運行結(jié)果：Expand]%]運行結(jié)果：3x.竺例2已知函數(shù)表x12345y141310711求出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))文檔來自于網(wǎng)絡(luò)搜索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文檔來自于網(wǎng)絡(luò)搜索L4兇〃NSimplify]%]L3[1.5]//N運行結(jié)果：插值多項式為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)文檔來自于網(wǎng)絡(luò)搜索插值多項式的簡式為4 3 20.208333x-1.75x 4.29167x-4.75x16插值為L4(1.5)=13.6797。例3已知14個點的坐標(biāo)如下表，根據(jù)表中數(shù)據(jù)構(gòu)造一插值多項式，并由此計算x=4.25點的函數(shù)值的近似值(不要求給出插值多項式)。 文檔來自于網(wǎng)絡(luò)搜索解已知數(shù)據(jù)表：1 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.02.71828 4.48169 7.38906 12.1825 20.0855 33.1155 54.5982 54.5982 90.0171Mathematica程序：A=Table[{x,Exp[x]},{x,1,5,0.5}]//Ng1=ListPlot[Table[A],Prolog->AbsolutePointSize[6],DisplayFunction->Identity];lnterpolation[A,lnterpolationOrder->1]g2=Plot[%兇,{x,1,5},PlotStyle->Thickness[0.005],DisplayFunction->Identity]文檔來自于網(wǎng)絡(luò)搜索Show[g1,g2,DisplayFunction->$DisplayFunction]%%%[4.25]運行結(jié)果：所求插值f(4.25)：70.3076。例4已知函數(shù)y二sinx的函數(shù)表如下x00.10.20.30.40.50.6sinx00.099830.198670.295520.389420.4794260.56464(1)分別用線性插值、二次插值和三次插值求sin0.57891的近似值，并估計截斷誤差;(2)試用4次等距節(jié)點插值公式計算 f(0.48)的近似值，并估計誤差。解Mathematica語句：Table[a={{0.3,0.4,0.5,0.6,},{0.29552,0.38942,0.47943,0.56464}}];文檔來自于網(wǎng)絡(luò)搜索Transpose[a][[{2,3}]];y1=InterpolatingPolynomial[%,x]//ExpandTranspose[a][[{2,3,4}]];y2=InterpolatingPolynomial[%,x]//ExpandTranspose[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],Thicknessp005]},文檔來自于網(wǎng)絡(luò)搜索PlotLegend->{"y1","y2","y3"},LegendPosition->{1,0}]文檔來自于網(wǎng)絡(luò)搜索Map[NumberForm[#,10]&,{Sin[0.48],y1,y2,y3}/.x->0.48]文檔來自于網(wǎng)絡(luò)搜索

運行結(jié)果:一次插值多項式二次插值多項式一次插值多項式二次插值多項式三次插值多項式2_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.5resError[0.48,2,{0.4,0.5,0.6},0.4]/.t->0.4文檔來自于網(wǎng)絡(luò)搜索resError[0.48,3,{0.3,0.4,0.5,0.6},0.3]/.t->0.3文檔來自于網(wǎng)絡(luò)搜索以下是三種插值在計算sin0.57891時的誤差的運行結(jié)果0.00038354-0.0000294744.2554910z(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]}, 文檔來自于網(wǎng)絡(luò)搜索{0,0,f2[0],f2[1],f2[2],f2[3],f2[4]},{0,0,0,f3[0],f3[1],f3[2],f3[3]},文檔來自于網(wǎng)絡(luò)搜索{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]}}；文檔來自于網(wǎng)絡(luò)搜索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]; 文檔來自于網(wǎng)絡(luò)搜索NN[x]=Sum[a[k]*Product[(t-j),{j,0,k-1}],{k,0,6}]〃N 文檔來自于網(wǎng)絡(luò)搜索N[Expand[%],3]〃Chop%/.t->0.48運行結(jié)果：所得差分表為：0000000文檔來自于網(wǎng)絡(luò)搜索0.099830.0998300000文檔來自于網(wǎng)絡(luò)搜索0.19870.09884-0.000990000文檔來自于網(wǎng)絡(luò)搜索0.29550.09685-0.00199-0.001000文檔來自于網(wǎng)絡(luò)搜索0.38940.0939-0.00295-0.000960.0000400文檔來自于網(wǎng)絡(luò)搜索0.47940.09001-0.003894-0.0009440.000016-0.0000240文檔來自于網(wǎng)絡(luò)搜索0.56460.08521-0.004792-0.0008980.0000460.000030.000054文檔來自于網(wǎng)絡(luò)搜索插值多項式為：0.09983t-0.00099(-1.+t)t-0.001(-2.+t)(-1.+t)t+0.00004(-3.+t)(-2.+t)(-1.+t)t文檔來自于網(wǎng)絡(luò)搜索-0.000024(-4.+t)(-3.+t)(-2.+t)(-1.+t)t+0.000054(-5.+t)(-4.+t)(-3.+t)(-2.+t)(-1.+t)t文檔來自于網(wǎng)絡(luò)搜索化簡得：0.0915t 0.0184t2-0.0142t3 0.00487t4-0.000834t5 0.000054t6則f(0.48) 0.468457。例5過0,1兩點構(gòu)造一個三次插值多項式，滿足條件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]))DI[x_,k_]:=D[l[x,k],x]/.x->x[k]h[x_,k_]:=(1-2(x-x[k])(DI[x,k]/.x->x[k]))(I[x,k]A2);文檔來自于網(wǎng)絡(luò)搜索H[x_,k_]:=(x-x[k])(I[x,k]A2);TabIe[h[x,k],{k,0,1}];TabIe[H[x,k],{k,0,1}];HH[x_,n_]:=Sum[y[k]*h[x,k]+y'[k]H[x,k],{k,O,n}];HH[x,1]Expand]%]運行結(jié)果： 2(1—x)2x(1-2(—1 x))x2 2(—1x)x2化簡為： 2x-3x22x3例6求數(shù)據(jù)表123456789一0.0.0.x-10.80.60.40.20246y-0.0.0.1.2.2.3.3.3.320942959826239200004445033456017787的最小二乘二次擬合多項式。解最小二乘求解的Mathematica程序：CIear兇X=TabIe[-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.1564};文檔來自于網(wǎng)絡(luò)搜索q[a_,b_,c_]:=Sum[(a+b*X[[k]]+c*X[[k]]A2-Y[[k]])A2,{k,1,10}] 文檔來自于網(wǎng)絡(luò)搜索SoIve[{D[q[a,b,c],a]==0,D[q[a,b,c],b]==0,D[q[a,b,c],c]==0},{a,b,c}]文檔來自于網(wǎng)絡(luò)搜索運行結(jié)果：擬合多項式系數(shù){{a->2.49261,b->2.42838,c->-0.351241}}擬合多項式22.49261+2.42838X+0.351241XLL=TabIe[{X[[k]],Y[[k]]},{k,1,9}]g1=ListPIot[LL,ProIog->AbsoIutePointSize[10]]f=Fit[LL,{1,x,xA2},x]g2=Plot[f,{x,-1,1}]Show[g1,g2]運行得擬合效果：例7（雙曲擬合）給定數(shù)據(jù)xi.01.11.21.31.41.5000000y.971.772.597.429.168.0651求形如、二的擬合函數(shù)。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}};文檔來自于網(wǎng)絡(luò)搜索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]運行得到結(jié)果：a°一-26.2461，ai=24.686。從而得到擬合曲線方程為

的函數(shù)，現(xiàn)測得函數(shù)值表如下:X.1c.20.300000.80.90.4.5.6.7y.34c.570.8901.34.4911.80.552.483.005試用最小二乘法確定y0和a。擬合曲線為例9設(shè)有一形為y=aQ擬合曲線為例9設(shè)有一形為y=aQ+a1x-26.2461 24.686X目二y°e-ax解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.89,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]-1.18014運行結(jié)果為：-1.18014+3.0968x。即丿-&=3.0968A=-1.18014;b=3.0968;a=Exp[A];y=a*Exp[b*x]

g仁ListPlot[L,Prolog->AbsolutePointSize[5]]g2=Plot[y,{x,0,0.9}]Show[g1,g2]運行得結(jié)果為： y=e%=0.307236,a=-a1=J.0968則原函數(shù)可近似為： y=0.307236e3.0968x擬合曲線為：第三章線性方程組的解法復(fù)習(xí)題例3.1用高斯消元法解方程組x2-0.1x3X4二2.70.5x24x3-8.5&=21.9

-x2x35.2x4二-3.9

x10.2x22.5x3-x4=9.9解取小數(shù)位4位對（A,b）做高斯消元：21-.112.70.40.54-8.521.90.3-115.2-3.910.22.5-19.9(A,b)=21-0.112.7 '00.34.02-8.721.360-1.151.0155.05-4.3050-0.32.55-1.58.551丿f21-0.112.7 、00.34.02-8.721.360016.425-28.377.575k006.57-10.229.91』21-0.112.7 '00.34.02-8.721.360016.425-28.377.5750001.12-1.12丿=(A,b)由同解方程組AX二b作回代過程有:X1'5.00、X22.00X33.001X4丿「1.00」X二例3.2用列主元素法解方程組:4xi-X2X3二4xi*-18xi+3x2—X3=-15、Xi+X2+X3=6解對(A,b)做選主元及消元過程:

廣4-115、f_183-1-仗(A,b)=—183-1-15T4-115<1116」116>-183-1-150175_3930-183-1-150175_393071731I6186f-183-1-150717316186002222k217丿=(A,b)-18-11718-1531‘1.0000、X2=2.000060000』例3.3分別用Doolittle和Crout分解法求解方程組2%x2x3=4*%+3x2+2x3=6當(dāng)+2x2+2x3=521P解方程系數(shù)矩陣A=132，常數(shù)項b=6J22>(1)Doolittle分解:對等式41P『1f、2U11U12U13A=132=l211U22U23=LU<122」31l32bkU33丿按U的第一行、L的第一列，按U的第二行、L的第二列，U的第三行的次序，根據(jù)矩

陣乘法與矩陣相等，對比等號兩邊矩陣元素，易得矩陣 L,U元素：文檔來自于網(wǎng)絡(luò)搜索‘2 11'L=0.5 1，U=2.5 1.5卩.5 0.6 1‘< °6丿求解LY二b，得(y1,y2,y3)=(4,4,0.6)。求解UX二Y，得(x1,x2,x3^(1,1,1)。(2)Crout分解z211、『1U12U13對等式A=132=l21l222U23=LUJ22」a31l32133丿3丿根據(jù)矩陣乘法與矩陣相等，對比等號兩邊矩陣元素，易得矩陣 L,U元素:2、10.50.5"L=0.52.5，u=10.6U1.5 0.6j1丿求解LY=b，得(y1,y2,y3^(2.0,1.6,1.0)。求解UX-Y，得(x「X2，X3)二(1.0,1.0,1.0)。例3.4 用矩陣的直接三角分解法解方程組*1020]"x1、01011X231243X317Q1 03丿'JJ分析 這是常規(guī)計算題，只需按矩陣的三角分解過程認(rèn)真計算即可。解設(shè)1020、「1\‘1020、0101l21 1U22 U23 U241243131 l32 1U33 U342103」1J41 142 G 1丿< U44丿

廣1\■1l21101l31I32112141I42I43b1°10b102002U22U23U2410U33U342U44■■I)由矩陣乘法可逐行，逐列分別求出uij和lij:0112)■1、A、廣5、01y2_3121y317?101』M丿I7」有yi=5,y=3,y3=6,y4=4得X得X4=2, X3=2，X2=1，*1020"101X2321X36<2丿gXi—1O注 由于三角分解過程十分規(guī)律，也可按緊湊格式直接得到彳0205、，Z10205、0101301013T124 31712216Q1 0 37丿少1024丿46一x42,X3 -2,x2=3-x4-1,X122其中折線右上部的元素（豎線左側(cè)）為上三角陣u的元素，最后一列為X4y的元素，故=5_2x3=1例3.5用平方根法（Cholesky分解）解方程組‘32<3U丿分析 由于系數(shù)矩陣對稱正定，故一定有分解形式A=LL，其中L為下三角陣,然后由矩陣乘法即可求出L的元素?！?3、l11l11l21J220=l21l22l22l32V3012」J31l32l33丿l33丿設(shè)由矩陣乘法得l11I21一3’l3122「23,G—6，I33=3解得再由<ll11l21l22y2=3Q31l32 l33}A丿l7丿511v'3'y246，。■1M1l21 l31p1l22 l32X2=y2l33J區(qū)3丿3yi=1。1X2 ，X12例3.111.5丿的1范數(shù)、2范數(shù)和無窮范數(shù)。解由向量范數(shù)定義得3|x||i訂|XkK=1=1+2+1.5=4.5X2=.'X2k=.1 4 2.25―7.25=2.69258x廣max風(fēng)鼻2[11_2Frobenius范數(shù)。例3.12計算矩陣A= 的1范數(shù)、2范數(shù)、sFrobenius范數(shù)。(25 —3.5丿解由矩陣范數(shù)定義：A1=max{1.1+2.5, -2+3.5}=5.5\A=max{1.1+-2,2.5+3.5}=6rnn 2" . IAf=巨Z|aij +22+2.52+3.52=4.86929Vyj絲T 7.46 -10.95、AA=1-10.9516.25/ata特征值入=23.6541,h=0.05591A2=7人=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]]A2,{i,1,2},{j,1,2}]]〃NT=Transpose[A].A;MatrixForm[%];Eigenvalues[A];A2=N[Sqrt[Max[Eigenvalues[T]]],10]運行結(jié)果為：5.56.4.869294.86354706例3.16用雅可比方法求解方程組x1x23x3=3*4捲+x2_x3=5

2%+5x2+2x3=-4方程的準(zhǔn)確解為(2,-2,1)T。解(1)由系數(shù)矩陣TOC\o"1-5"\h\z'1 1 3A= 4 1 -1<2 5 2構(gòu)造雅可比迭代矩陣：f0_1-3、『3、B=—401,g=5<_1-2.50」廠2」X(k1)=BX(k) g，取X(0)=(1,-1,1)T。計算結(jié)果見下表:kx(k)x1x(k)x1x(k)x1|x(k)-x(鬥仁11.50.5-87.529.83333-9-4.759.5315.3056-39.083310.666730.0833413.5648-45.555680.402869.73615-55.174431.143598.324176.69916-144.488324.022-24.6844292.8787-178.648558.268-667.566642.8818363.37852.0271-1219.02551.45791444.93-2667.53-495.4452719.56

可以看到迭代發(fā)散。B的特征值。入二-3.55707,k=1.77854 2.23374i,入=1.77854-2.23374i所以有(B)1(2)將方程的次序調(diào)換可以得到TOC\o"1-5"\h\z'4 12 5d 1構(gòu)造迭代格式:x嚴(yán)」(0—x2k)+x3k)+5)4視宀=丄(—2x1k)+0_2x3k)_4)5x3“」(-膚-x2k)+0+3)I3寫成迭代矩陣形式:f(k*)f(k*)\X1x2?虻)1、、4■‘X1(k)、42x2k)+5一5(k)40丿嚴(yán)丿1L J迭代收斂，取X(0)=(1,-1,1)T。結(jié)果見下表:k(k)x1(k)x1x1k) |x(k)-x(kV11.75-1.610.7521.9-1.90.950.331.9625-1.9410.062541.985-1.9850.99250.04551.99438-1.99110.00937561.99775-1.997750.9988756.75X10」71.99916-1.9986511.40625x10^81.99966-1.999660.9998311.0125"0‘

91.99987-1.999812.10938x10,101.99995-1.999950.9999751.51875X10°111.99998-1.9999713.16406x10」121.99999-1.999990.9999962.2781200'例3.17用高斯一塞德爾方法求解方程組4x-|X2_X3二52x「5x22x3=-4x1x23x3二3其中方程的準(zhǔn)確解為(2,-2,1)T。解高斯一塞德爾迭代格式:Xi("J(0_x2k)+x3k)+5)4+O-2x3k)-4)5Xi(fJ(-Xi(f-x2?+0+3)3寫成迭代矩陣形式:(k4t)(k4t)X2111-5144_X1(k)1411(k)13—— X2+— 102(k)1011N一612012一[.60一因為r(S)二0.182574,?(SM?(B)。所以，高斯一塞德爾迭代格式收斂速度比雅可比迭代快。取X(0)=(1,1,1)T為初始值，下表列出迭代計算結(jié)果：kX1k)x2k)v(k)X3l|x(k)-x(k%11.25-1.71.152.721.9625-2.0451.02750.712532.01813-2.018251.000040.05562542.00457-2.001850.9990910.016404252.00023-1.999730.9998320.0043387261.99989-1.999890.999999-43.43695X10471.99997-1.999991.000019.966X10-582.-21-62.64189X10Mathematica程序：A={{4,1,-1},{2,5,2},{1,1,3}};MatixForm[%];b={5,-4,3};l1=ldentityMatix[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=lnverse[DD-L];G=I1-DDLN.A;MatrixForm[%];p=Max[Abs[Eigenvalues[N[G]]]]R=p-1f=DDLN.b;x[0]={1,1,1};x[n_]:=G.x[n-1]+f;Table[x[n],{n,1,8}];MatrixForm[%]//NTable[Max[Abs[x[n]-x[n-1]]],{n,1,8}];MatrixForm[%]//NLinearSolve[A,b]例3.18如何對方程組-Xi 8X2 OX3=7?—Xi+OX2+9x3=89xi-X2-X3=7進(jìn)行調(diào)整，使得用高斯-塞德爾方法求解時收斂？并取初始向量 x(0)=(0,0,0)T，用該方法求近似解x(k1)，使x(kJx(k\.<10"。解將第三個方程調(diào)到第一行后有9Xi—X2—X3=7?-為+8x2+0x3=7一捲+0x2十9x3=8這是主對角線嚴(yán)格占優(yōu)方程組，故用高斯 -塞德爾迭代法求解一定收斂。迭代格式為Xi(k—1(x2k)+x3k)+7)銳宀=8(嚴(yán)+7)才出=1(xi燦+8)由x(0)=(0,0,0)丁，得x⑴=(0.7778,0.9722,0.9753)Tx(2)=(0.9942,0.9993,0.9994)tx(3)=(0.9999,0.9999,0.9999)tx(4)=(1.0000,1.0000,1.0000)T因為 x(4)-x(3)二=0.0001=10,"0“，故取最后結(jié)果為

x(4)=(1.0000,1.0000,1.0000)T例3.佃給定方程組證明：雅可比方法發(fā)散，而高斯-11證明：雅可比方法發(fā)散，而高斯-11丫“:9111JX2★1-2人X3丿-塞德爾收斂。證明對雅可比方法，迭代矩陣為1

"1

"2

-1-1

丄-1

丄012設(shè)其特征值為’則>3“5>It,‘J4 1(BH-2 1，故雅可比方法發(fā)散。對高斯-塞德爾方法,迭代矩陣為-111~2-10J)f 、對高斯-塞德爾方法,迭代矩陣為-111~2-10J)f 、101-1、(11、0 2P-11220-1 0-1-1220110112丿1)0 0 飛<2丿'\、227故高斯-塞德爾收斂。31，孫2顯然，其特征值為 \=0,'2例3.20設(shè)A為正交矩陣,B=21A，求證：線性方程組BTBx=例3.20設(shè)A為正交矩陣,爾方法求解必收斂。

證明：因A為正交矩陣，故ATA=I，設(shè)■是A的特征值，則有x=0，使Ax=x兩邊與Ax作內(nèi)積，有(Ax,Ax)=(Ax,x),(x,ATAx)=2(x,x)從而有(x,x)=&2(x,x)，(丸2—1)(x,x)=0,岡=1，又B的特征值為2-丸，故B的特征值不為零，從而B非奇異。故BtB正定，所以高斯-塞德爾方法收斂。例3.21例3.21對方程組片「x2=b|Xi2x2=b2給出解方程組的雅可比迭代矩陣，并討論迭代收斂條件；給出解方程組的高斯-塞德爾迭代矩陣，并討論迭代收斂條件。解(1)雅可比迭代矩陣:-p'「0-p'B=p<~2其特征多項式為:PB(PB(莎二p2故A的譜半徑:P(B)=占v2即雅可比迭代收斂的充要條件為 p<J2。(2)高斯-塞德爾迭代矩陣:of0-of02p2丿S的特征多項式為:

故S的譜半徑:(S)r故S的譜半徑:(S)r2~2P匸)即高斯-塞德爾迭代收斂的充要條件為 p<<2?；贛athematica的數(shù)值計算實例TOC\o"1-5"\h\z3% +4x2 _5x3+7x4 =0，‘ 亠——2為—3x2 +3x3—2x4 =0例1求解線性萬程組/ 1 2 3 44x<h11x2一13x3+16x4=0了為—2x2+x3+3%=0解Mathematica程序：Det[A]=0,A={{3,4,-5,7},{2,-3,3,-2},{4,11,-13,16},{7,-2,1,3}};文檔來自于網(wǎng)絡(luò)搜索MatrixForm[%]Det[A]r=Sum[RowReduce[A][[i,i]],{i,1,4}]NullSpace[A]運行結(jié)果為：行列式等于0,系數(shù)矩陣的秩為2，有非零解，基礎(chǔ)解系為{{-13,-20,0,17},{3,19,17,0}}文檔來自于網(wǎng)絡(luò)搜索可得方程組的通解為：X*{-13,-20,0,17}k2{3,19,17,0}={-13k13k2^20k119k2,17k2,17k1}|-1x13x22x3-X4=0例2求解線性方程組3咅_5x2_3x32x4二0

x-t3x27x35x例2求解線性方程組I7X--3x?-5X3-X4=0解Mathematica程序:A={{-1,3,2,-1},{3,-5,-3,2},{1,3,7,5},{7,-3,-5,-1}};文檔來自于網(wǎng)絡(luò)搜索MatrixForm[%]Det[A]NullSpace[A]運行結(jié)果為：行列式的值為例3求解線性方程組_36嚴(yán)0，方程組只有零解。2片+x2_4x3_x4=86片_2x2_2x3+2x4=415x2+14x3+3x4=—4

4x--3x2_IOX3-x^—8解Mathematica程序:Clear[A,b]A={{2,1,-4,-1},{6,-2,-2,2},{0,5,14,3},{4,-3,-10,-1}}; 文檔來自于網(wǎng)絡(luò)搜索MatrixForm[%];b={8,4,-4,8};Det[A]r=Sum[RowReduce[A][[i,i]],{i,1,4}]NullSpace[A]LinearSolve[A,b]運行結(jié)果為：行列式等于o,系數(shù)矩陣的秩為3,有非零解，基礎(chǔ)解系為{{-1,1,-1,3}}，一個特解為文檔來自于網(wǎng)絡(luò)搜索{1,2,-1,0}。則原方程組的通解為:X特征值與特征例4方陣A由前36個相鄰整數(shù)構(gòu)成，求A的特征多項式、特征方程、向量。

特征值與特征解Mathematica程序:A=Partition[Range[36],6];MatrixForm[%]CharacteristicPolynomial[A,-]CharacteristicPolynomial[A,■]==0Eigenvalues[N[A]]〃ChopEigenvectors[A]〃N;MatrixForm[%]運行結(jié)果：特征多項式特征方程特征根文檔來自—630?！?11特征多項式特征方程特征根文檔來自-630-4-1115 6=00,0,0}0001、001001001000-0.2767740.1488170.57440910.4767740.6511830.8255911i特征向量{116.412,-5.41182,0,特征向量4-54-53-42-31-27.12796-0.7023662.1279570.3023662x3y例5用高斯-若當(dāng)(Gauss-Jordan)消元法解線性方程組*x+4y-x_2y+解Mathematica程序:A={{2,3,1,6},{1,4,-1,4},{1,-2,1,1}};MatrixForm[%]RowReduce[A];MatrixForm[%]運行結(jié)果：1001580105001838」

則求得原方程組的解為:15則求得原方程組的解為:15x二8例6設(shè)向量x=(2,-5,1,3,0)t，求向量范數(shù)xp,p=1,2,二。解Mathematica程序：x={2,-5,1,3,0};x1=Sum[Abs[x[[i]]],{i,1,5}]運行結(jié)果：11x2=Sqrt[x.x]運行結(jié)果：39x:=Max[Table[Abs[x[[i]]],{i,1,5}]]運行結(jié)果：5例7設(shè)矩陣TOC\o"1-5"\h\z\o"CurrentDocument"2 0 1\o"CurrentDocument"A= -1 0 -1 4\o"CurrentDocument"2 1 1 3<2 3 15」求矩陣A的范數(shù)||Xp,p=1,2嚴(yán)。解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}]]文檔來自于網(wǎng)絡(luò)搜索T=Transpose[A].A;MatrixForm[%];Eigenvalues[N[T]];A2=N[Sqrt[Max[Eigenvalues[N[T]]]]]運行結(jié)果： 13 11 8.25543例8考察n階希爾伯特(Hilbert)矩陣Hn的性態(tài)11…11TOC\o"1-5"\h\z2 n1 1 … 1Hn=2 3 n十1\o"CurrentDocument"a a + 3丄丄…1』n+1 2n-仁解Mathematica通用程序：Clear[B]B=Table[1/(i+j),{i,1,k},{j,0,k-1}];MatrixForm[%]BN=lnverse[B];B1=Max[Sum[Abs[B[[n]]],{n,1,k}]];BN1=Max[Sum[Abs[BN[[n]]],{n,1,k}]];condA仁B1*BN1〃NB100=Max[Table[Sum[Abs[B[[n,m]]],{m,1,k}],{n,1,k}]];文檔來自于網(wǎng)絡(luò)搜索BN100=Max[Table[Sum[Abs[BN[[n,m]]],{m,1,k}],{n,1,k}]];文檔來自于網(wǎng)絡(luò)搜索condA100=B100*BN100〃N當(dāng)k取4時，運行結(jié)果： 28375 283757捲-4x2 2x3=18例9用Jacobi迭代法求解方程組 』4捲+10x2—x3=40?Xr+3x2+9x3=29解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]]}}文檔來自于網(wǎng)絡(luò)搜索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}]Print[k,xb,,z];xa=xb,{k,1,10}]文檔來自于網(wǎng)絡(luò)搜索運行結(jié)果：0.454927即Jacobi迭代法迭代矩陣的譜半徑小于 1,迭代法收斂。n 近似解 誤差{2.57143,4.,3.22222} 4.{3.93651,3.29365,0.174603} 3.04762{4.40363,2.44286,-0.5} 0.8507944{4.1102,2.18855,-0.527816}0.2934245{3.97283,2.30314,-0.24743}0.2803856{3.9582,2.38612,-0.194045