《數(shù)學實驗 第4版》課件 第二章 方程及方程組的解

數(shù)學實驗國家“十二五”規(guī)劃教材山東省優(yōu)秀教材一等獎2第二章方程及方程組的解實驗2.1線性方程組的解

實驗2.2非線性方程的解數(shù)學實驗3若系數(shù)矩陣的秩，則方程組有唯一解；

線性方程組的通解=對應齊次方程組的通解+非齊次方程組的一個特解,其特解的求法屬于解的第一類問題，通解部分屬第二類問題．我們將線性方程組的求解分為兩類：一類是求方程組的唯一解或求特解;

另一類是求方程組的無窮解即通解．可以通過系數(shù)矩陣的秩來判斷：若系數(shù)矩陣的秩，則方程組可能有無窮解；數(shù)學實驗4二、應用舉例一、解線性方程組的MATLAB實現(xiàn)實驗2.1線性方程組的解數(shù)學實驗5一、解線性方程組的MATLAB實現(xiàn)1．求線性方程組的唯一解或特解

若是線性方程組的矩陣形式，則可以直接用矩陣的除法或初等變換法求線性方程組的唯一解或特解。例1

求方程組的解．A=[56000;15600;01560;00156;00015];B=[10001]';R_A=rank(A)%求秩X=A\B↙%求解解

方法一實驗2.1線性方程組的解6R_A=5X=2.2662-1.72181.0571-0.59400.3188運行后結果如下方法二

用函數(shù)rref求解：C=[A,B]%由系數(shù)矩陣和常數(shù)列構成增廣矩陣CR=rref(C)↙%將C化成行最簡形R=1.000000002.266201.0000000-1.7218001.0000001.05710001.00000-0.594000001.00000.3188則R的最后一列元素就是所求之解．7例2

求方程組的一個特解．A=[11-3-1;3-1-34;15-9-8];B=[140]';X=A\B↙%由于系數(shù)矩陣不滿秩，該解法可能存在誤差解警告:秩虧，秩=2，tol=8.837264e-15.X=[00-0.53330.6000]’

可以驗證所求X是方程組一個特解的近似值.這是因為系數(shù)矩陣奇異，方程組的解不存在或者不唯一，通過下述方式判斷方程組是否有解.實驗2.1線性方程組的解8R=1.00000-1.50000.75001.250001.0000-1.5000-1.7500-0.250000000因此，該方程組有無窮解.要求得一個特解，可使用函數(shù)rrefA=[11-3-1;3-1-34;15-9-8];B=[140]';C=[A,B];%構成增廣矩陣

R=rref(C)↙C=[A,B];R_A=rank(A)R_C=rank(C)R_A=2R_C=2輸出結果為令

，得方程組的一個特解X=(1.2500–0.250000)'．由此可見，不同的求解命令適用的條件不同，需要先判定適用的命令再進行求解.92.求線性齊次方程組的通解

在Matlab中，函數(shù)null用來求解零空間，即滿足A·X=0的解空間，實際上是求出解空間的一組基（基礎解系）．基本格式：z=null(A)%z的列向量為方程組的正交規(guī)范基，滿足z=null(A,′r′)%z的列向量是方程AX=0的有理基實驗2.1線性方程組的解10例3

求方程組的通解：解A=[1221;21-2-2;1-1-4-3];formatrat%指定有理式格式輸出B=null(A,'r')↙%求解空間的有理基運行后顯示結果如下：B=25/3-2-4/31001進一步寫出通解，在命令窗口輸入：實驗2.1線性方程組的解symsk1k2X=k1*B(:,1)+k2*B(:,2)%寫出方程組的通解pretty(X)↙11運行后結果如下：X=[2*k1+5/3*k2][-2*k1-4/3*k2][k1][k2]實驗2.1線性方程組的解或通過行最簡形得到基：B=rref(A)↙B=1.00000-2.0000-1.666701.00002.00001.33330000123．求非齊次線性方程組的通解

非齊次線性方程組需要先判斷方程組是否有解，若有解，再去求通解．因此，步驟為：第四步：的通解為的通解加上的一個特解．是否有解，若有解則進行第二步第一步：判斷的一個特解第二步：求第三步：求的通解實驗2.1線性方程組的解13例4

求解方程組解在Matlab中建立M文件如下：A=[1-23-1;3-15-3;212-2];b=[123]';B=[Ab];n=4;R_A=rank(A)R_B=rank(B)formatratifR_A==R_B&R_A==n%判斷有唯一解

X=A\belseifR_A==R_B&R_A<n↙%判斷有無窮解實驗2.1線性方程組的解14X=A\b%求特解

C=null(A,'r')%求AX=0的基礎解系elseX='equitionnosolve'%判斷無解end運行后結果顯示：R_A=2R_B=3X=equitionnosolve說明該方程組無解.15例5

求方程組的通解：解法一在Matlab編輯器中建立M文件如下：A=[11-3-1;3-1-34;15-9-8];b=[140]';B=[Ab];n=4;R_A=rank(A)R_B=rank(B)formatrat ifR_A==R_B&&R_A==nX=A\b elseifR_A==R_B&&R_A<nC=rref(B)%求增廣矩陣的行最簡形elseX='Equationhasnosolves'end實驗2.1線性方程組的解16對應齊次方程組的基礎解系為：非齊次方程組的特解為：所以，原方程組的通解為：運行后結果顯示為：R_A=2R_B=2C=10-3/23/45/401-3/2-7/4-1/40000017二、應用舉例

１．投入產(chǎn)出綜合平衡分析

國民經(jīng)濟各個部門之間存在著相互依存的關系，每個部門在運轉(zhuǎn)中將其他部門的產(chǎn)品或半成品經(jīng)過加工（稱為投入）變成自己的產(chǎn)品（稱為產(chǎn)出），問題是如何根據(jù)各部門之間的投入產(chǎn)出關系，確定各部門的產(chǎn)出水平，以滿足社會的需要，下面考慮一個簡化的問題：

設國民經(jīng)濟僅有農(nóng)業(yè)、制造業(yè)和服務業(yè)三個部門構成，已知某年它們之間的產(chǎn)出關系、外部需求、初始投入等如表2.1所示（數(shù)字表示產(chǎn)值，單位為億元）．實驗2.1線性方程組的解18

假定每個部門的產(chǎn)出與投入是成正比的，由表2.1能夠確定這三個部門的投入產(chǎn)出表，如表2.2農(nóng)業(yè)制造業(yè)服務業(yè)外部需求總產(chǎn)出農(nóng)造業(yè)301045115200服務業(yè)2060/70150初始投入3511075總投入100200150表2.1國民經(jīng)濟各個部門間的關系產(chǎn)出投入表2.2投入產(chǎn)出表農(nóng)業(yè)制造業(yè)服務業(yè)農(nóng)業(yè)0.150.100.20制造業(yè)0.300.050.30服務業(yè)0.200.300產(chǎn)出投入實驗2.1線性方程組的解19（1）設有個部門，記一定時期內(nèi)第個部門的總產(chǎn)出為個部門的投入為，滿足的外部需求為,則：（2.1）記第個部門的單位產(chǎn)出需要第個部門的投入為.在每個部門的產(chǎn)出與投入

成正比的假定下，有：

（2.2）式（2.2）代入式（2.1）得

（2.3）實驗2.1線性方程組的解，其中對第20記投入系數(shù)矩陣，產(chǎn)出向量，需求向量

，則（2.3）式可寫作

（2.4）（2.5）可逆，則：（2.6）或若和外部需求給定后，即可算出各部門的總產(chǎn)出.當投入系數(shù)實驗2.1線性方程組的解21（2）將表2.2中的投入系數(shù)即對各部門的外部需求輸入MATLAB，用以下程序：formata=[0.150.10.2;0.30.050.3;0.20.30];d=[50150100]';b=eye(3)-a;x=b\dc=inv(b)xx=c*d↙實驗2.1線性方程組的解輸出結果為x=139.2801267.6056208.1377c=1.34590.25040.34430.56341.26760.49300.43820.43041.2167xx=139.2801267.6056208.1377即三個部門的總產(chǎn)出分別應為139.2801，267.6056，208.1377（億元）22（3）式（2.6）表明總產(chǎn)出對外部需求是線性的，增加一個單位時，單位時，農(nóng)業(yè)、制造業(yè)和服務業(yè)的總產(chǎn)出應分別增加1.3459，0.5634，0.4382單位，其余類似，這些數(shù)字稱為部門關聯(lián)系數(shù)．的增量由決定，其第一列數(shù)字表明，當農(nóng)業(yè)的需求增加一個實驗2.1線性方程組的解所以當23２．大型輸電網(wǎng)絡設電源電壓為一種大型輸電網(wǎng)絡可簡化為圖2.1所示電路，其中表示線路內(nèi)阻，圖2.1輸電網(wǎng)絡簡化電路(1)

列出各負載上電流的方程；記上的電流為

（2.7）實驗2.1線性方程組的解表示負載電阻，24（2.８）消去得：

（2.9）25

（2.9）記則方程（2.9）表示為：式（2.9）或式（2.10）即是電流滿足的方程．（2.10）實驗2.1線性方程組的解26（2）設，及總電流的情況下，我們可以利用以下程序求出r=1;R=6;v=18;n=10;b1=sparse(1,1,v,n,1);b=full(b1);a1=triu(r*ones(n,n));a2=R*eye(n);a3=-tril(R*ones(n,n),-1)+tril(R*ones(n,n),-2);a=a1+a2+a3;I=a\b，I0=sum(I)↙輸出結果為：I=2.00051.33440.89070.59550.39950.27020.18580.13240.10110.0867I0=5.9970實驗2.1線性方程組的解，在27（3）考慮的情況.為求出總電阻，考察圖2.2所示的第n段電路（從右向左數(shù)）和第段電路的等效電阻和，有：（2.11）圖2.2輸電網(wǎng)絡的等效電阻以代入可計算，當時記為，滿足：即：由此解得:(2.12)實驗2.1線性方程組的解28而不難得出：

（2.13）于是總電流：實驗2.1線性方程組的解29第二章方程及方程組的解實驗2.1線性方程組的解

實驗2.2非線性方程的解30實驗2.2非線性方程組的解數(shù)學實驗一、二分法二、簡單迭代法三、牛頓迭代法四、應用舉例31實驗2.2非線性方程的解一、二分法1.數(shù)學原理32離區(qū)間,重復上述做法，實驗2.2非線性方程的解332.MATLAB實現(xiàn)用二分法求方程根的近似值，一般分為兩步：第一步根的隔離用fplot函數(shù)作出y=f(x)的圖形，它與x軸的交點就是方程f(x)=0的根.確定根的存在區(qū)間[a,b](即根的隔離區(qū)間)使方程第二步根的逐次逼近若函數(shù)f(x)在根的隔離區(qū)間[a,b]上單調(diào)連續(xù),且則在(a,b)內(nèi)方程有唯一根.我們可以反復用將根的隔離區(qū)間一分為二，其中含有根的區(qū)間為新的隔離區(qū)間的方法，最終搜索到達到要求精度的根的近似值.實驗2.2非線性方程的解34實驗2.2非線性方程的解functiony=erfen(fun,a,b,esp)iffeval(fun,a).*feval(fun,b)<0n=1;c=(a+b)/2;

whileabs(b-a)>esp

if

feval(fun,a).*feval(fun,c)<0b=c;c=(a+b)/2;

elseiffeval(fun,c).*feval(fun,b)<0a=c;c=(a+b)/2;

elsey=c;esp=10000;

end打開M-edit窗口,編寫二分法程序：erfen·mn=n+1;endy=c;elseiffeval(fun,a)==0y=a;elseiffeval(fun,b)==0y=b;elsedisp('these,maynotbearootintheintercal');endn35實驗2.2非線性方程的解解fplot('x^4-8.6*x^3-35.51*x^2+464.4*x-998.46',[0,10]),gridon↙例6

求方程的隔離區(qū)間，并用二分(1)求根的隔離區(qū)間，設法求方程在這個區(qū)間內(nèi)根的近似值，使誤差不在命令窗口輸入：36實驗2.2非線性方程的解經(jīng)過15次搜索，得到方程的近似根為x=7.3485，誤差不超過functiony=f(x)y=x^4-8.6*x^3-35.51*x^2+464.4*x-998.46;erfen('f',7,8,1.0e-4)↙在命令窗口輸入:編寫函數(shù)文件f.mn=15ans=7.3485

由二分法的原理可知，當根的隔離區(qū)間較大時，收斂于根的速度很快；而當根的隔離區(qū)間較小時，靠近根的速度較緩慢，且難以滿足高精度的要求.因此常用這種方法為迭代法提供初值.37實驗2.2非線性方程的解二、簡單迭代法簡單迭代法的基本思想是將方程改寫為等價形式當方程的根難以求出時，這種求根的近似值的方法稱為簡單迭代法.稱為迭代格式.1.簡單迭代法簡介38實驗2.2非線性方程的解用迭代法求方程的根的首要問題是迭代序列是否收斂.為此我們先從幾何上觀察由簡單迭代格式得到的迭代序列的收斂性.由圖可見，根據(jù)初始值因此迭代序列是收斂的.39在下圖中，按上述同樣的方法產(chǎn)生的點列逐漸遠離交點

，可見迭代序列是發(fā)散的.通過上面的作圖可見，迭代序列是否收斂與a的取值有關，實驗2.2非線性方程的解40實驗2.2非線性方程的解2.簡單迭代法MATLAB實現(xiàn)先編寫簡單迭代程序iterate.m及關于方程的函數(shù)f.m文件存入工作窗口,再調(diào)用其中x0為初始值.編寫迭代法程序iterate.m:簡單迭代法命令iterate('f',x0),functiony=iterate(fun,x,esp)ifnargin<3esp=1e-4;endx1=feval(fun,x);n=1;while(abs(x1-x)>=esp)&(n<=1000)

x=x1;

x1=feval(fun,x);

n=n+1;endx1n其中espw為精度.41，分別用下面兩種迭代格式：例7

已知方程求方程根的近似值，使誤差不超過解在迭代法程序iterate.m中，令esp=1.0e-4.(1)編寫函數(shù)g.m文件functiony=g(x)y=(0.8+x^2)^(1/3);在命令窗口輸入：iterate('g',1.5)↙實驗2.2非線性方程的解x1=1.4052n=1042functiony=g1(x)y=(x^3-0.8)^0.5;(2)編寫函數(shù)g1.m文件在命令窗口輸入：迭代格式(2)發(fā)散.iterate(1.5)↙實驗2.2非線性方程的解經(jīng)過10次迭代,得到方程的近似根為x=1.4052,誤差不超過x1=

Infn=2043在迭代格式(1)中，設，則在迭代格式(2)中，設，則迭代格式（2）是發(fā)散的.迭代格式是收斂的.實驗2.2非線性方程的解44定理則在該鄰域內(nèi)的任意初始值收斂于α.（證明略）但如何改寫才能使迭代序列收斂？收斂的速度更快？這是用迭代法求方程的根的另一個問題.注意到迭代函也是他們的加權平均的不動點.根據(jù)上面的實驗,使得迭代序列收斂并加快收斂速度,實驗2.2非線性方程的解45得于是得修正迭代函數(shù)修正迭代格式實驗2.2非線性方程的解46在例7中，g2.mfunctiony=g2(x)y=x-[3*(0.8+x^2)-3*x*(0.8+x^2)^(2/3)]/[2*x-3*(0.8+x^2)^(2/3)]iterate('g2',1.5)↙輸入：迭代格式（1）修正后，經(jīng)過2次迭代，得到方程的近似根

x=1.4052，誤差不超過10-4，收斂速度明顯加快.利用修正迭代格式求根的近似值如下：實驗2.2非線性方程的解x1=1.4052n=247編寫函數(shù)文件g3.m:functiony=g3(x)y=x-[2*(x^3-0.8)-2*x*(x^3-0.8)^(1/2)]/[3*x^2-2*(x^3-0.8)^(1/2)]iterate('g3',1.5)↙輸入：迭代格式（2）修正后，經(jīng)過2次迭代，就可以達到所求的近似根.實驗2.2非線性方程的解x1=1.4052n=248實驗2.2非線性方程的解三、牛頓迭代法牛頓迭代法是將非線性方程f(x)=0的求解逐步化為線性方程來求解，是線性化與迭代法的結合.從幾何上看,就是用曲線上一點的切線與x軸的交點來近似代替根.1.非線性方程單根的求法如圖49重復這一過程得迭代格式這種求方程近似根的方法稱為牛頓迭代法.迭代函數(shù)為實驗2.2非線性方程的解502.非線性方程重根的求法若存在正整數(shù)m使得可以證明，在充要條件是實驗2.2非線性方程的解由此可見，則即513.牛頓迭代法的MATLAB實現(xiàn)編寫牛頓迭代法程序newton.m:functiony=newton('f',df,x0)x1=x0-feval(f,x0)/feval(df,x0);esp=1e-4;n=1;while(abs(x1-x0)>=esp）

x0=x1;x1=x0-feval(f,x0)/feval(df,x0);

n=n+1;endx1n實驗2.2非線性方程的解52例8

用牛頓迭代法求方程在[0,1]內(nèi)實根的近似值,使誤差不超過解在牛頓迭代法程序中，取編寫函數(shù)文件f.m:functiony=f(x)y=exp(x)+10*x-2;編寫f(x)的導函數(shù)文件df.m:functiony=f(x)y=exp(x)+10;實驗2.2非線性方程的解53分別取初值x0=1,0,在命令窗口輸入：newton（'f','df',1）↙newton（'f','df',0）↙當初值x0=1時，迭代4次后得到滿足要求的的近似解

x=0.0905;

當初值x0=0時，僅迭代3次后得到所求的解.由此可見,選擇適當?shù)某踔悼梢蕴岣叩蛄械氖諗克俣?實驗2.2非線性方程的解x1=0.0905n=4x1=0.0905n=354例9

求方程附近的根，使誤差不超過解由例6知，方程在x=4附近有一重根.在牛頓迭代法中,編寫函數(shù)文件f1.m:functiony=f1(x)y=x^4-8.6*x^3-35.51*x^2+464.4*x-998.46;實驗2.2非線性方程的解編寫f1(x)的導函數(shù)文件df1.m:functiony=df1(x)y=4*x^3-25.8*x^2-71.02*x+464.4;在命令窗口輸入：newton（'f1','df1',4.0）↙x1=4.2999n=

12經(jīng)過12次迭代，得到滿足要求的解

x=4.2999.55用公式編寫函數(shù)文件f2.m:functiony=f2(x)y=(x^4-8.6*x^3-35.51*x^2+464.4*x-998.46)/(4*x^3-25.8*x^2-71.02*x+464.4);實驗2.2非線性方程的解編寫f2(x)的導函數(shù)文件df2.m:functiony=df2(x)y=1-(x^4-8.6*x^3-35.51*x^2+464.4*x-998.46)*(12*x^2-51.6*x-71.02)/(4*x^3-25.8*x^2-71.02*x+464.4)^2;在命令窗口輸入：newton（'f2','df2',4.0）↙經(jīng)過3次迭代，得到滿足要求的解

x=4.3000.x1=4.3000n=356從圖中可以看出，在[3,5]內(nèi)f1(x)=0有一個重根.如果選取的初值x0不恰當，點的切線斜率的絕對值就會較大，此時用切線與x軸的交點x1作為根的近似值，與用x0作為根的近似值產(chǎn)生的誤差相差無幾,因此用f1(x)構造的牛頓迭代序列收斂性較差（迭代50次才達到要求的精度?。。?而用f2(x)構造的迭代序列收斂速度非?？欤▋H需要迭代3次就可以達到要求的進度.實驗2.2非線性方程的解

數(shù)值求解非線性方程的方法雖然比較機械，有時需要進行大量的重復計算，但其詮釋了無限逼近的過程，精確與近似的相對統(tǒng)一.57實驗2.2非線性方程的解四、應用舉例一輛汽車售價A0元，可分m個月付款，每月需交b元，計算這種分期付款的年利率.利用這種方式購車，我們需要每月把b元存入汽車銷售商的賬戶,也就相當于汽車銷售商周期的收入資金.設利率按每月收款計算為r,為使m個月后存入銀行A0元，則第一個月應該交第二個月應該交于是有58由公式（1）可以計算出月利率，月利率乘以12便可以得到年利率.若買一輛汽車15000元,可分60個月付款，每月需交3114元;也可以向銀行貸款,貸款的期限為5年,年利率為6.66%,應該選擇哪一種方式購車.例10解在公式（1）中,于是得顯然r=0是上述方程的根，但這里我們要求的是最小正根.實驗2.2非線性方程的解編寫函數(shù)文件f.m59編寫函數(shù)文件f.mfunction

y=f(r)y=150000*r*(1+r)^60-3114*(1+r)^60+3114;確定根的隔離區(qū)間fplot('f',[0.001,0.5]),gridon↙由圖知根的隔離區(qū)間［0.001,0.5］.實驗2.2非線性方程的解60調(diào)用二分法的程序,在命令窗口輸入:erfen('f',0.001,0.5,1.0e-4)↙于是得月利率年利率為所以分期付款買車要付的利率高達9%，與5年貸款的利率比較，最好選擇貸款買車.實驗2.2非線性方程的解n=

14ans=0.007561例1