数学实验-线性代数方程组的数值解

实验5线性代数方程组的数值解法数学实验报告实验5线性代数方程组的数值解法数学实验报告5线性代数方程组的数值解法1黄浩2011011743一、实验目的MATLAB定性作初步分析；通过实例学习用线性代数方程组解决简化的实际问题。二、实验内容《数学实验》第二版（问题叙述：通过求解线性方程组𝐴1𝑥=𝑏1,𝐴2𝑥=𝑏2，理解条件数的意义和方程组性态对解的影响，其中𝐴1是n阶范德蒙矩阵，即1

𝑥20

… 𝑥𝑛10𝐴 = 1

𝑥2

… 𝑥𝑛1

𝑥=1+

0,1,…,𝑛 11 …

1 … … 𝑘[1

𝑥2𝑛1

… 𝑥𝑛1]𝑛1]𝐴2nb1，b2𝐴1𝐴2的行和。编程构造𝐴1（𝐴2可直接用命令产生）和你能预先知道方程组𝐴1𝑥=和𝐴2𝑥=𝑏2（用预先知道的解可验证程序。n=5,7,9，…，计算𝐴1和𝐴2的条件数。为观察他们是否病态，做以和𝐴2的元素𝐴1(n分别加扰动ε后求解；𝐴1和𝐴2不变，b1，b2b1(n),b2(n)分别加扰动εAbε10^-10,10^-8,10^-6。经扰动得到的解记做‖𝑥𝑥‖𝑥‖模型转换及实验过程：小题.b1,b2𝐴1,𝐴2的行和，可知方程组𝐴1𝑥=和𝐴2𝑥=𝑏2n1n=5matlab（.11,2x1,x2excel算而得）为：xx1e1x2e211.0000000000000800.0000000000000801.0000000000000200.00000000000002010.999999999999463-0.0000000000005370.999999999999558-0.00000000000044211.0000000000009700.0000000000009701.0000000000020100.00000000000201010.999999999999288-0.0000000000007120.999999999996863-0.00000000000313711.0000000000001800.0000000000001801.0000000000015600.000000000001560由上表可见，当n=5时，所得的解都接近真值，误差在10^-12的量级左右。（2）小题n=5,7,9,11,13,151和2c1c2（程序见四结果如下：n579111315c13.574E+058.739E+072.274E+106.518E+122.059E+158.291E+17c24.766E+054.754E+084.932E+115.227E+146.263E+173.675E+17𝐴1和𝐴2b)（.3，结果如下：(为使结果清晰简洁，在此仅列出n=5,9,13的情况，n=7,11,15略去)ε=10^-10时：nxx1xx250.9999999284994020.9999999370003031.0000002510851691.0000012599939300.9999996704145830.9999943300273341.0000001916677021.0000088199574730.9999999583331440.999995590021263nxx1xx290.9999975548050940.9999829278035681.0000149244303261.0012291970545280.9999603223369260.9784890657946091.0000600154317051.1577467696177110.9999435079279020.4084498576816771.0000338876411232.2304238808731590.9999873483708230.4354941236522751.0000026878626961.8788737437097660.9999997511934060.780281609027228nxx1xx2131.0007784033869290.9999998767583270.9937743987274311.0000261185116641.0226872410459610.9987794790026500.9501840256039321.0233044427317511.0734106035353810.7678410435503330.9235081059886532.3669384594613431.0577906781685654.0983444092691000.96809962810653313.4983690949393491.01276967899037019.4039025417116410.99638464083603822.9640269119150101.00068723718260013.9700648285337990.9999212432336596.8528369055604181.0000041151939490.000188901053803=10^-8nxx1xx250.9999928500302660.9999937027771001.0000251082268421.0001259444579840.9999670418066980.9994332499390881.0000191665851081.0008816112058550.9999958333510870.999559194397072nxx1xx290.9997569218861830.9999249646886321.0014841958519581.0054025423430740.9960527298593050.9054555097425781.0059726416246391.6933262585415040.9943760878476931.5999734614530311.0033746909787466.4079447885959670.9987396985502615.3092689110300421.0002678327406924.8628176967655060.9999752006605240.034295576497007nxx1xx2130.9980046593334370.9999998767282631.0164113342717761.0000261243251180.9384398540328910.9987792205553051.1392634494355081.0233092230919600.7883771891288210.7677944673667701.2275767009039882.3672082173739580.8224007796560294.0993377962873331.10134460129783313.5007796714360190.95802934170017619.4078052094758501.01230286523370722.9681988664444180.99757692077678013.9728915924515431.0002879131499591.0002879131499590.9999843910790946.8539364954676760.000001889299579ε=10^-6时：nxx1xx250.9992852977929170.9993966095201831.0025097875871361.0120678095963350.9967055393606440.9456948568165061.0019158683868581.084474667174326n0.999583506872445xx10.957762666412834xx290.9757501420061600.9999223275065541.1480658167117691.0055924194987990.6062136691657740.9021326589030921.5958409910913531.7176938349075140.4389487716335661.6913518835480071.3366650248743926.5980119249310320.8742701611529635.5310139213940691.0267194461452354.998579957366803n0.997525977218788xx10.000355009332673xx2130.7071952366081760.9999998767279653.3904929638964631.0000261243832617.8992555052967080.99877921797036220.9772809601275481.02330927090318629.1195663172520300.76779400153262433.1329169063691252.36721091537023223.8745993777145214.09934773169287015.07902114649785113.5008037809264344.78290586997535419.4078442421848052.68118057389295422.9682405924371480.67162490808792313.9729198644592591.0386948092493446.8539474930660780.9979195655092280.000000018893025由上表可见：n=96n=1322ε。n<101.6n>1033（因此，我认为第一版上的参考答案不够全面）Ab大，随ε的增加，解的偏差也有增加的趋势。方程组的解表格在此略去，程序也类似，其具体影响及分析见下一小问。（3）小题A、b以通过计算误差‖𝑥𝑥̃‖，来与条件数估计的误差比较，来进一步得出结论。‖𝑥‖2-norm（x）bA(n,n)以微小的扰动（程序见四.4（E1、E2分别为两方程解的误差‖𝑥𝑥̃‖）：‖𝑥‖nε=10^10ε=10^8ε=10^6条件数E1E2E1E2E1E2Cond(A1)Cond(A2)52.07E075.12E062.07E055.12E042.07E034.90E023.57E+054.77E+0573.19E062.10E033.19E041.89E013.19E021.74E+008.74E+074.75E+0893.30E057.28E013.30E033.20E+003.30E013.31E+002.27E+104.93E+11112.19E045.94E+002.56E026.02E+002.55E+006.02E+006.52E+125.23E+14132.51E024.71E+001.25E014.71E+001.59E+014.71E+002.06E+156.26E+17152.31E+007.84E+001.79E+007.51E+006.19E+015.96E+008.29E+173.68E+17由上表可得出结论：第一，随着阶数的增加，误差值有增加的趋势第二，ε的变化对误差大小也有影响，但不是绝对的，对于n=5,7,9,11,13，εn=15，找不到εE第三，实际算得的误差远远小于条件数指出的误差上限第四，当n<10时，微小扰动对范德蒙矩阵的影响较小，但当n>10时，微小n=15尔伯特阵类似。再保持𝐴1和𝐴2不变，b1，b2b1(n),b2(n)分别加扰动ε（程序见四.5）结果如下：E1E2E1E2E1E2Cond(A1)Cond(A2)52.07E-075.12E-062.07E-055.12E-042.07E-035.12E-023.57E+054.77E+0573.19E-062.10E-033.19E-042.10E-013.19E-022.10E+018.74E+074.75E+0893.33E-059.33E-013.30E-039.33E+013.30E-019.33E+032.27E+104.93E+11111.62E-044.31E+022.57E-024.31E+042.56E+004.31E+066.52E+125.23E+14133.55E-022.49E+041.31E-012.49E+061.59E+012.49E+082.06E+156.26E+17151.06E+013.73E+049.61E-013.73E+067.01E+013.73E+088.29E+173.68E+17nε=10^10εnε=10^10ε=10^8ε=10^6条件数第一，随着阶数的增加，误差值有增加的趋势第二，ε越n=5,7,9,11,13εεE第三，实际算得的误差远远小于条件数指出的误差上限n<10n>10bA是，最终范德蒙阵的病态程度也没有赶上希尔伯特阵bA造成的结果偏差更大《数学实验》第二版（问题叙述：已知方程组Ax=b，其中A∈R^(20×20)，定义为：试通过迭代法求解此方程组，认识迭代法收敛的含义以及迭代初值和方程组系数矩阵性质对收敛速度的影响。实验要求：选取不同的初始向量x(0)b用雅可比迭代法和高斯-赛德尔迭代法计算，观测得到的迭代向量序列是否均收敛？若收敛，记录迭代次数，分析计算结果并得出你的结论；bx(0)A非主对角线元素不变，每次用雅可比迭代法计算，要求迭代误差满足‖x(k+1)x(k)‖ <比较收敛速度，分析现象并得出你的结论。∞模型转换与实验过程：小题matlab设计雅克比迭代法和高斯-赛德尔迭代法的函数，设为JacDD1GSDD1（.6、四.7Abx0mtol，xJ（xG）k。首先，保持常数项向量b=(1,1,1,…,1)T不变，分别设定初始向量x(0)(0,0,0…,0)T(5,5,5,…,5)T(100,100,100,…,100)T（程序见四.8，所得结果如下：x(0)x(0)解向量(0,0,0,…,0)T(5,5,5,…,5)T(100,100,100,…,100)TxJxGxJxGxJxG0.4816350110.4816335080.4816371820.4816409430.4816366890.4816398830.5734131580.5734116110.5734168540.5734210570.5734160170.5734194990.6327989970.6327976730.632804280.6328083480.6328030890.6328062750.6520986790.6520977760.6521053730.6521088050.6521038710.6521062630.6609418470.6609414470.6609497970.6609523360.6609480230.6609494090.6643052030.6643053610.664314220.6643156670.664312220.6643125230.6657142460.6657150040.6657241290.6657243440.6657219490.6657211990.6662729030.6662742810.666283440.6662823630.6662811270.6662794310.6664962310.6664982130.6665072060.6665048770.6665048050.6665023220.6665735250.6665760480.6665847190.666581280.6665822740.666579190.6665735250.6665764810.6665847190.6665803960.6665822740.6665787840.6664962310.6664994830.6665072060.6665022740.6665048050.6665010990.6662729030.6662762940.666283440.6662781920.6662811270.6662773790.6657142460.6657176190.6657241290.665718850.6657219490.6657183160.6643052030.6643084110.664314220.664309170.664312220.6643088380.6609418470.6609447570.6609497970.66094520.6609480230.6609450050.6520986790.6521011780.6521053730.6521014210.6521038710.6521013130.6327989970.6328009950.632804280.6328011180.6328030890.6328010640.5734131580.5734145680.5734168540.5734146220.5734160170.5734145980.4816350110.4816358440.4816371820.4816358630.4816366890.481635855171119122415迭代次数由上表可见，随着初值x(0)迭代次数但是二者的结果相近，都在误差上限之内。为了证明两种求解方法的收敛性，我们以x(0)=(5,5,5,…,5)𝑇时的高斯-赛德（程序见四似，在此略去：k=0k=1k=2k=3k=4k=5k=651.5833333330.8097993830.5867626890.5159863770.4929626410.48538750751.8472222220.9783629120.7052563650.6168153850.5877768120.57818087752.0231481481.084426440.7813237890.6819209230.6490964610.63821656752.0744598771.1206981240.8071642940.7035546140.6692025920.65779093752.0976723251.1372138040.8191053420.7135320420.6784443550.66677151852.1058170441.1436399750.8237874540.7173972820.6819885220.67019847552.1091088681.1463317770.8257794470.7190360460.683482730.67163852852.1103362321.1474109130.8265897910.7196988670.6840822740.67221365952.1108151111.1478512890.8269263590.7199736250.6843293130.67244801252.1109972041.1480285630.8270643090.7200860580.684429750.67252176652.111067461.1481001480.8271211030.7201323670.68445080.67243591352.1110943441.1481289050.8271443740.7201513680.6842499120.67203420152.1111046791.148140450.827153910.719918050.6832395860.67110569252.1111086421.1481450730.8271578080.7178314910.6805376480.66952080652.1111101641.1481469220.8242658810.7102614320.675688690.66701021152.1111107471.148147660.805458660.6951998230.6686267110.66268609152.1111109721.1134257320.7566388370.6735497460.6567218090.65311633152.1111110580.9629628860.6907658520.6445060140.6352482410.63333137351.6944444240.7185570680.600404170.5787075670.5745170960.57365217950.7916666590.5333397520.4909645160.4834934290.4820235360.481719644k=7k=8k=9k=10k=11k=120.4828815270.4820500490.4817736970.4816817550.4816511440.4816409430.5749989270.5739417540.573590110.5734730570.5734340650.5734210570.6346027310.6334008590.6330008350.632867610.6328231840.6328083480.6539954920.6527322070.6523115110.6521712730.6521244310.6521088050.6628857660.6615916860.6611604250.6610163940.6609682580.6609523360.6662714780.6649629750.6645261480.6643799480.6643313890.6643156670.6676920040.6663752910.6659342820.6657871530.6657392950.6657243440.6682561060.6669308290.666486750.666341410.666295940.6662823630.6684736390.6671370010.6666960450.6665576910.6665165820.6665048770.6685091180.6671722840.6667505750.666625920.6665908350.666581280.6683888310.6671001890.6667211670.6666159640.666587780.6665803960.6680849140.6669276810.6666119130.6665289910.6665076820.6665022740.6675508740.6666014440.6663582090.666297130.6662819480.6662781920.6666577650.6659471590.6657736310.6657315260.6657213240.665718850.6649462430.6644591990.6643443560.6643171770.6643107140.664309170.6613428570.6610368410.6609663530.6609499580.660946110.66094520.6523289070.6521531220.6521132330.6521040570.6521019220.6521014210.6329187230.6328276480.6328071450.6328024570.6328013720.6328011180.5734671370.5734264380.5734173010.5734152160.5734147340.5734146220.4816544160.4816400440.4816368120.4816360740.4816359030.481635863迭代次数增加k越来越小，是逐渐收敛的。然后，保持初始向量x(0)=(2,2,2,…,2)T不变，分别设定常数项向量b(0.01,0.01,0.01,…,0.01)T(1,1,1…,1)T(5,5,5…,5)T（xJ，b(0.01,0.01,0.01,…,0.01)T(1,1,1,…,1)Tb(0.01,0.01,0.01,…,0.01)T(1,1,1,…,1)T(5,5,5,…,5)T解向量kxJxGxJxGxJxG0.0048163670.0048164390.48163670.4816405622.4081724482.4081652010.005734160.0057342440.5734160320.5734205922.8670613492.8670550730.006328030.0063281150.6328031040.6328078783.1639886493.1639849690.0065210380.0065211130.6521038780.6521084383.2604853863.2604853060.006609480.0066095390.6609480160.6609521713.304699743.3047035850.0066431220.0066431640.6643121940.6643157723.321515273.3215231710.0066572190.0066572440.6657219010.665724723.3285594793.3285715020.0066628110.0066628190.6662810590.6662829453.3313520093.3313681320.0066650470.0066650420.6665047220.6665055613.3324681483.3324881660.0066658220.0066658060.6665821830.6665819633.3328543633.3328778060.0066658220.0066657970.6665821830.6665810013.3328543633.3328804820.0066650470.0066650170.6665047220.6665027633.3324681483.3324959760.0066628110.0066627780.6662810590.6662785563.3313520093.3313804540.0066572190.0066571860.6657219010.6657191013.3285594793.3285874180.0066431220.006643090.6643121940.6643093323.321515273.3215416260.006609480.0066094510.6609480160.6609452983.304699743.3047235330.0065210380.0065210140.6521038780.6521014763.2604853863.2605057490.006328030.0063280110.6328031040.6328011463.1639886493.1640049050.005734160.0057341460.5734160320.5734146342.8670613492.8670728070.0048163670.0048163590.48163670.4816358682.4081724482.40817921251518111510由上表可见，随着b少，而且雅克比迭代法比高斯-赛德尔迭代法所需的迭代次数要多，但是二者的结果相近，都在误差上限之内。x=（0,0,0,…，0）T1042545综上所述，我们可以得到以下几点结论：x(0)迭代次数增多。这是因为迭代收敛时，迭代值逐次逼近真值，因而初值与真值越接近，就可以在较少迭代次数的情况下满足误差上限。在迭代收敛的情况下，当b0迭代次数减少。迭代次数要多。因矩阵的谱半径不超过（任一种）范数，可求得：ρ(Bj)≤‖𝐵𝑗‖2=0.4893062,ρ (Bg)≤‖𝐵𝑔‖2=0.3292286显然，二者都是收敛的，而且后者的收敛速度快。小题我们固定常数项向量b=初始向量x(0)=(1,1,1,…,1)T，将A的主对角线元素𝑥𝑖𝑖解向量3,9,15,21,27,331,3,5,7,9,11（程序见四.10𝑥𝑖𝑖解向量39152127330.4816367890.1219637270.0703231230.0494345270.0381191020.0310207030.5734161790.1289623430.0727104420.0506266140.0388319270.0314944760.6328033050.1327656850.073965230.0512443130.039198370.0317368180.6521041190.1331791860.0740476840.0512734260.0392118290.0317441330.6609482860.1333090390.0740714350.0512814990.0392154770.0317460970.664312480.1333280570.0740736230.0512820480.0392156720.0317461930.6657221950.1333327940.0740740980.0512821620.0392157110.0317462150.6662813540.133333610.0740741530.0512821740.0392157140.0317462180.6665050150.1333337970.0740741640.0512821760.0392157140.0317462190.6665824740.1333338330.0740741660.0512821770.0392157140.031746220.6665824740.1333338330.0740741660.0512821770.0392157140.031746220.6665050150.1333337970.0740741640.0512821760.0392157140.0317462190.6662813540.133333610.0740741530.0512821740.0392157140.0317462180.6657221950.1333327940.0740740980.0512821620.0392157110.0317462150.664312480.1333280570.0740736230.0512820480.0392156720.0317461930.6609482860.1333090390.0740714350.0512814990.0392154770.0317460970.6521041190.1331791860.0740476840.0512734260.0392118290.0317441330.6328033050.1327656850.073965230.0512443130.039198370.0317368180.5734161790.1289623430.0727104420.0506266140.0388319270.0314944760.4816367890.1219637270.0703231230.0494345270.0381191020.031020703迭代迭代次数1687665由上表可得出结论幅度加快，但是减小的幅度也变得平缓。《数学实验》第二版（问题叙述：入产出关系、外部需求、初始投入等如下表所示：投入产出投入产出农业制造业服务业外部需求总产出农造业301045115200服务业2060070150初始投入3511075总投入100200150如果今年对农业、制造业和服务业的外部需求分别为50，150，100元，问这三个部门的总产出分别应为多少？如果三个部门的外部需求分别增加1增加多少？模型转换及实验过程：小题.设一定时期内第i个部门的总产出为𝑥𝑖，其中对第j个部门的投入为设消耗系数矩阵A=aij ，其中𝑎𝑖𝑗=𝑥𝑖𝑗÷则由题中表格，可作出三个部门间投入系数aij的表格：投入产出农业制造业服务业农业0.15 0.10.2制造业0.3 0.050.3服务业0.2 0.30设产出向量x= x1,x2,x3T，需求向量d= d1,d2,d3T，根据每个部门的总产出等于总投入，可以得到下列关系：x=Ax+d即：(I−A)x=dmatlab（程序见四.11果为：

d=(139.280,267.606,208.138)T50，150，100139.280、267.606、208.138小题.因为(IA)x=d，即x=(IA)−1d，因此：xx=(I−A)−1(d+J),J分别为(1,0,0)𝑇(0,1,0)𝑇(0,0,1)𝑇即：总产出的增加额dx=xx−x=(I−A)−1J需求增加产出增加农业制造业服务业matlab(I)x=J需求增加产出增加农业制造业服务业农业1.350.560.44制造业0.251.270.43服务业0.340.491.22得出结论：11.35、0.56、0.4410.25、1.27、0.4310.34、0.49、1.22三、实验总结本次数学实验是求解线性代数方程组的数值解。因素，在方程不合适时，收敛步数可能会达到上千步甚至更多。习中，通过更实际化的数模练习，相信会更加熟练地掌握这一技能。四、程序清单1.第一题（1）n=5时的高斯消元解n=5;k=0:n-1;x=1+0.1*k;fori=0:(n-1)A1(i+1,:)=x.^i;endA1=A1';A2=hilb(n);x1=A1\b1x2=A2\b22.第一题（2）输出n=5,7,9,11,13,15的条件数forj=1:6n=3+2*j;A1=ones(n);k=0:n-1;x=1+0.1*k;fori=0:(n-1)A1(i+1,:)=x.^i;C

endA1=A1';A2=hilb(n);c1=cond(A1);c2=cond(A2);C(j,:)=[nc1c2];x（n=5,9,13的情况，n=7,11,15bAAb10^-10）clear;forj=1:3n=1+4*j;A1=ones(n);k=0:n-1;x=1+0.1*k;fori=0:(n-1)A1(i+1,:)=x.^i;

endA1=A1';A11=A1;A11(n,n)=A11(n,n)+10^-10;A2=hilb(n);A22=A2;A22(n,n)=A22(n,n)+10^-10;b1=A1*ones(n,1);b2=A2*ones(n,1);x11=A11\b1;x22=A22\b2;[x11x22]第一题（3）bAclear;forj=1:6n=3+2*j;A1=ones(n);k=0:n-1;x=1+0.1*k;fori=0:(n-1)A1(i+1,:)=x.^i;endA1=A1';