數(shù)值計(jì)算方法實(shí)驗(yàn)報(bào)告

1、數(shù)值計(jì)算方法實(shí)驗(yàn)報(bào)告112010XXXXXXX數(shù)值分析實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)一、解線性方程組的直接方法梯形電阻電路問(wèn)題利用追趕法求解三對(duì)角方程組的方法，解決梯形電阻電路問(wèn)題：電路中的各個(gè)電流 i1 ， i 2 ， ， i8 須滿足下列線性方程組：2i12i2V R2i 15i22i302i25i32i 402i35i 42i 502i45i 52i 602i 55i 62i 702i65i 72i 802i75i 80設(shè) V220V ， R 27，運(yùn)用追趕法，求各段電路的電流量。問(wèn)題分析：上述方程組可用矩陣表示為：22000000i18.148125200000i 2002520000i30002520

2、00i 4000025200i5000002520i6000000252i7000000025i801數(shù)值計(jì)算方法實(shí)驗(yàn)報(bào)告112010XXXXXXX問(wèn)題轉(zhuǎn)化為求解Axb ， 8 階方陣 A 滿足順序主子式Ai (i1,2.7)0 ，因此矩陣A存在唯一的Doolittle 分解，可以采用解三對(duì)角矩陣的追趕法！追趕法a=0 -2 -2 -2 -2 -2 -2 -2;b=2 5 5 5 5 5 5 5；c=-2 -2 -2 -2 -2 -2 -2 0;d=220/27 0 0 0 0 0 0 0;Matlab 程序functionx= zhuiganfa( a,b,c,d )%追趕法實(shí)現(xiàn)要求：|b1

3、|C1|0,|bi|=|ai|+|ci|n=length(b);u=ones(1,n);L=ones(1,n);y=ones(1,n);u(1)=b(1);y(1)=d(1);fori=2:nL(i)=a(i)/u(i-1);u(i)=b(i)-c(i-1)*L(i);y(i)=d(i)-y(i-1)*L(i);endx(n)=y(n)/u(n);fork=n-1:-1:1x(k)=(y(k)-c(k)*x(k+1)/u(k);endendMATLAB命令窗口輸入：a=0 -2 -2 -2 -2 -2 -2 -2;b=2 5 5 5 5 5 5 5;c=-2 -2 -2 -2 -2 -2 -2

4、 0d=220/27 0 0 0 0 0 0 0;x= zhuiganfa(a,b,c,d )運(yùn)行結(jié)果為：x =8.14784.07372.03651.01750.50730.25060.11940.0477存在問(wèn)題根據(jù)電路分析中的所講到的回路電流法，可以列出8 個(gè)以回路電流為獨(dú)立變量的方程，課本上給出的第八個(gè)回路電流方程存在問(wèn)題，正確的應(yīng)該是2i74i80 ；或者可以根據(jù)電路并聯(lián)分流的知識(shí)，同樣可以確定2i 74i80 。正確的處理結(jié)果應(yīng)為：x =8.14814.07412.03701.01850.50930.25460.12730.06372數(shù)值計(jì)算方法實(shí)驗(yàn)報(bào)告112010XXXXXXX實(shí)

5、驗(yàn)二、線性方程組的迭代方法不同迭代法解線性方程組的對(duì)比分別用（ 1）Jacobi迭代法；（ 2）Gauss-Seidel迭代法；（ 3）共軛梯度法解線性方程組101234x11219123x22721735x314323121x417435115x512迭代初始向量取 x(0)T00000；(1)Jacobi迭代法function x,k = jacobi(A,b,x)error=norm(b-A*x); %2-范數(shù)%x為迭代初始向量if (errore)e=input( 請(qǐng)輸入絕對(duì)誤差限:n );break;D=diag(diag(A);endn=size(b);endI=eye(n,n);

6、sprintf( 迭代次數(shù) :nk= %d ,k)B=I-DA;sprintf( 方程組的解為 :n )g=Db;xfor k=1:50;%最大迭代次數(shù)為50次endx=B*x+g;Matlab 命令窗口輸入A=10 1 2 3 4;1 9 -1 2 -3;2 -1 7 3 -5;3 2 3 12 -1;4 -3 -5 -1 15;b=12 -27 14 -17 12;x=0 0 0 0 0;x,k = jacobi(A,b,x);運(yùn)行結(jié)果請(qǐng)輸入絕對(duì)誤差限:10(-3)迭代次數(shù) :k= 42ans =方程組的解為 :x =1.0002-2.00022.9997-1.99990.9998(2)G

7、auss-Seidel迭代法function x,k = gau_sei(A,b,x)U=D-triu(A);D=diag(diag(A);e=input( 請(qǐng)輸入絕對(duì)誤差限:n );L=D-tril(A);for k=1:50;3數(shù)值計(jì)算方法實(shí)驗(yàn)報(bào)告112010XXXXXXXx=(D-L)U*x+(D-L)b;enderror=norm(b-A*x);sprintf( 迭代次數(shù) :nk= %d ,k)if (error10(-3)sprintf( 方程組的解為 :n )break;xendendMatlab 命令窗口輸入A=10 1 2 3 4;1 9 -1 2 -3;2 -1 7 3 -5

8、;3 2 3 12 -1;4 -3 -5 -1 15;b=12 -27 14 -17 12;x=0 0 0 0 0;x,k = gau_sei(A,b,x);運(yùn)行結(jié)果請(qǐng)輸入絕對(duì)誤差限:10(-3)迭代次數(shù) :k= 24方程組的解為 :x =1.0002-2.00022.9997-2.00000.9998（ 3）共軛梯度法function x,k = gongetidu( A,b,x )x=x+l*d;e=input( 請(qǐng)輸入絕對(duì)誤差限:n );error=norm(b-A*x);d=b-A*x;if (errore)r=d;break;l=(d*r)./(A*d)*d);endx=x+l*d;

9、endfor k=1:50sprintf( 迭代次數(shù) :nk= %d ,k)r=b-A*x;sprintf( 方程組的解為 :n )t=-(A*d)*r)./(A*d)*d);xd=r+t*d;endl=(d*r)./(A*d)*d);Matlab 命令窗口輸入A=10 1 2 3 4;1 9 -1 2 -3;2 -1 7 3 -5;3 2 3 12 -1;4 -3 -5 -1 15;b=12 -27 14 -17 12;x=0 0 0 0 0;x,k = gongetidu(A,b,x );運(yùn)行結(jié)果請(qǐng)輸入絕對(duì)誤差限:10(-3)迭代次數(shù) :k= 44數(shù)值計(jì)算方法實(shí)驗(yàn)報(bào)告112010XXXXX

10、XX方程組的解為 :x =1.0000-2.00003.0000-2.00001.0000結(jié)果分析在選擇相同的誤差容限時(shí)，使用以上三種不同的方法來(lái)解方程組，對(duì)比迭代次數(shù)上的不同（ 1） Jacboi 迭代法 42 次（ 2） Gauss-Seidel24 次（ 3）共軛梯度法4 次，可以看出在收斂速度上共軛梯度法最快，Jacobi 迭代法最慢，對(duì)比方程組的精確解x=1 -2 3 -2 1，可以看出在收斂效果上共軛梯度法最好！5數(shù)值計(jì)算方法實(shí)驗(yàn)報(bào)告112010XXXXXXX實(shí)驗(yàn)三、 插值法龍格現(xiàn)象在代數(shù)插值中， 為了提高插值多項(xiàng)式對(duì)函數(shù)的逼近層度，常常增加節(jié)點(diǎn)個(gè)數(shù)，級(jí)提高多項(xiàng)式的次數(shù)，但由于高次

11、多項(xiàng)式插值不收斂，會(huì)產(chǎn)生Runge 現(xiàn)象，本實(shí)驗(yàn)使用函數(shù)為1y = 1+x2 ，通過(guò)對(duì)比Lagrange 插值法和分段線性插值法可以觀察到Runge 現(xiàn)場(chǎng)的產(chǎn)生與消失。實(shí)驗(yàn)程序函數(shù) fxfunction y = fx( x )y=1./(1+x*x);endlagrange插值方法function = lagrangechazhi( )f=f+feval(fx,s(q)*l; %求插值函數(shù)n=input( 區(qū)間等分份數(shù)：n );ends=-5+10/n*0:n; %給定的定點(diǎn)， fx為給y=1./(1+x.*x);定的函數(shù)e=y-f; %誤差x=-5:0.01:5;subplot(211);f

12、=0;plot(x,f, r,x,y,k)%作出插值函數(shù)曲線for q=1:n+1;legend(擬合曲線 ,實(shí)際曲線 );l=1;%求插值基函數(shù)gridon;for k=1:n+1;subplot(212);if k=q;plot(x,e,b);l=l.*(x-s(k)./(s(q)-s(k);legend(誤差曲線 );elsegridonl=l;endendend1區(qū)間等分份數(shù)為10 時(shí)，如下圖所示，黑色曲線為的實(shí)際圖形，紅色曲線即為采2用 Lagrange插值方法實(shí)際描繪出的曲線， 從圖中可以觀察到明顯的 Runge現(xiàn)象； Runge 現(xiàn)象的是因?yàn)檫x擇的插值函數(shù)并不收斂的緣故。6數(shù)值計(jì)

13、算方法實(shí)驗(yàn)報(bào)告112010XXXXXXX分段線性插值方法function = fenduanxianxingchazhi()elseclearl=0;n=input( 區(qū)間等分份數(shù)：n );ends=-5+10/n*0:n; %給定的定點(diǎn)， Rf 為給end定的函數(shù)endm=0;%ff=ff+1./(1+25*s(k)*s(k)*l;%求插值函數(shù)for x=-5:0.01:5;值ff=0;ff=ff+feval(fx,s(k)*l;%求插值函數(shù)值for k=1:n+1; % 求插值基函數(shù)endswitch km=m+1;case1f(m)=ff;if xs(n);plot(x,f, r,x,y

14、,k)%作出插值函數(shù)曲線l=(x-s(n)./(s(n+1)-s(n);legend(擬合曲線 ,實(shí)際曲線 );elsegridon;l=0;subplot(212);endplot(x,e,b);otherwiselegend(誤差曲線 );if x=s(k-1)&x=s(k)&x=s(k+1);l=(x-s(k+1)./(s(k)-s(k+1);7數(shù)值計(jì)算方法實(shí)驗(yàn)報(bào)告112010XXXXXXX區(qū)間等分為 20 份時(shí)，如下圖所示，黑色曲線為1的實(shí)際圖形，紅色曲線即為采用1+x2分段線性插值方法實(shí)際描繪出的曲線，從圖中可以觀察到未出現(xiàn)Runge 現(xiàn)象，可以觀察誤差曲線可以看到在將區(qū)間等分為20

15、 份的時(shí)候，和實(shí)際曲線的誤差限為0.5X10（ -2）；繼續(xù)提高等份份數(shù)，可以繼續(xù)降低誤差！8數(shù)值計(jì)算方法實(shí)驗(yàn)報(bào)告112010XXXXXXX實(shí)驗(yàn)四、數(shù)值積分考紐螺線考紐螺線的形狀像鐘表的發(fā)條，也稱回旋曲線，它在直角坐標(biāo)系中的參數(shù)方程為x(s)s1 at2dtcos0 2y(s)s1at 2 dtsin0 2曲線關(guān)于原點(diǎn)對(duì)稱。取 a=1，參數(shù) s 的變化范圍 -5，5，容許誤差限分別是10-6 和 10-10。選取適當(dāng)?shù)墓?jié)點(diǎn)個(gè)數(shù)，利用數(shù)值積分方法計(jì)算曲線上點(diǎn)的坐標(biāo)，并畫(huà)出曲線的圖形。實(shí)驗(yàn)程序 （程序中編寫(xiě)的積分方法根據(jù)數(shù)值積分中常見(jiàn)的Romberg求積公式得到）function X,Y = ka

16、oniuluoxian( )%命令窗口輸入X,Y = kaoniuluoxian()formatlong;T=zeros(20,20);a=input( 積分下限： n );b=input( 積分上限： n );e=input( 允許誤差限： n );c=linspace(a,b,1000);X=zeros(1,1000);Y=zeros(1,1000);fx1=inline( cos(x.2/2) );%考紐螺線方程 fx2=inline( sin(x.2/2) );%inline函數(shù)建立方程for m=1:1000k=1;h=c(m)-a;T(1,1)=h/2*(fx1(a)+fx1(c(

17、m); %計(jì)算積分上限為 c(m)時(shí)fx1 的積分值for i=1:20%最大迭代次數(shù)20次x=a+h/2:h:c(m)-h/2;T(k+1,1)=T(k,1)/2+h*sum(fx1(x)/2;for j=1:kT(k-j+1,j+1)=(4j*T(k-j+2,j)-T(k-j+1,j)/(4j-1);endif abs(T(1,k+1)-T(1,k)ebreak;endk=k+1;h=h/2;命令窗口輸入X,Y = kaoniuluoxian( );endX(m)=T(1,k+1);endfor m=1:1000k=1;h=c(m)-a;T(1,1)=h/2*(feval(fx2,a)+f

18、eval(fx2,c(m); %fev al用來(lái)求函數(shù) fx在給定點(diǎn)處的值fori=1:20%計(jì)算積分上限為c(m)時(shí) fx2的積分值x=a+h/2:h:c(m)-h/2;T(k+1,1)=T(k,1)/2+h*sum(feval(fx2,x)/2;for j=1:kT(k-j+1,j+1)=(4j*T(k-j+2,j)-T(k-j+1,j)/(4j-1);endif abs(T(1,k+1)-T(1,k)1e-6)%誤差限k=k+1;difmax=0.0;for i=2:hy-1for j=2:hx-1m=mesh1(i,j);if (m=2)oriv=v(i,j); %原電勢(shì)值v(i,j)

19、=1/4*(v(i-1,j)+v(i,j-1)+v(i+1,j)+v(i,j+1); % 拉普拉斯方程差分式dif=v(i,j)-oriv; %前后兩次迭代之差 dif=abs(dif);if (difdifmax)difmax=dif; %取誤差的最大值endendendendendsubplot(1,2,1);%figure(1);mesh(v);%畫(huà)三維曲面圖axis(-2,hx+3,-2,hy+3,0,100);title(等勢(shì)線三維分布圖);subplot(1,2,2);%figure(2);contour(v,13);%畫(huà)等電位線圖holdon;x=1:1:hx;y=1:1:hy;

20、xx,yy=meshgrid(x,y);%形成柵格Gx,Gy=gradient(v,0.5,0.5);%計(jì)算梯度quiver(xx,yy,Gx,Gy,r);%根據(jù)梯度畫(huà)箭頭axis(-2,hx+3,-2,hy+3);plot(1,1,hx,hx,1,1,hy,hy,1,1,k);%外框邊界線plot(CX1,CX1,CX2,CX2,CX1,CY1,CY2,CY2,CY1,CY1,k);%內(nèi)框邊界線title(橫截面等勢(shì)線分布圖 ); text(CX1+0.6,CY1+(CY2-CY1)/2, U=100V ,fo ntsize,10);text(hx/2,hy+1, U=0V ,fontsiz

21、e,10); text(hx/2,0, U=0V ,fontsize,10); text(-1.7,hy/2, U=0V ,fontsize,10); text(hx+0.4,hy/2, U=0V ,fontsize,10); gridon;holdoff ;end12數(shù)值計(jì)算方法實(shí)驗(yàn)報(bào)告112010XXXXXXX命令窗口輸入dianweifenbu();運(yùn)行結(jié)果：從上圖可以看出， 在靠近內(nèi)金屬波導(dǎo)時(shí)， 等勢(shì)線分布更加密集， 同時(shí)電力線垂直于金屬波導(dǎo)邊緣，符合實(shí)際同軸金屬波導(dǎo)中的場(chǎng)強(qiáng)分布規(guī)律！13數(shù)值計(jì)算方法實(shí)驗(yàn)報(bào)告112010XXXXXXX總結(jié)實(shí)驗(yàn)報(bào)告中涉及到的內(nèi)容主要包括數(shù)值計(jì)算方法課中所講解過(guò)的（1）特殊矩陣的三角分解使用追趕法解決電路中常見(jiàn)的電阻網(wǎng)絡(luò)問(wèn)題，通過(guò)列出回路電流方程可以發(fā)現(xiàn)方程組的系數(shù)矩陣恰好為三對(duì)角矩陣，既可以使用追趕法解決；（2）使用迭代法 Jacobi迭代法、 Gauss-Seidel 迭代法和共軛梯度法來(lái)解線性方程組，通過(guò)對(duì)同一線性方程組設(shè)置相同的迭代初始向量和迭代誤