數(shù)值計算方法實驗報告

1、數(shù)值計算方法實驗報告 112010XXXX XXX數(shù)值分析實驗報告實驗一、解線性方程組的直接方法梯形電阻電路問題利用追趕法求解三對角方程組的方法，解決梯形電阻電路問題：電路中的各個電流，須滿足下列線性方程組：設(shè)，運用追趕法，求各段電路的電流量。問題分析：上述方程組可用矩陣表示為：問題轉(zhuǎn)化為求解，8階方陣A滿足順序主子式，因此矩陣A存在唯一的Doolittle分解，可以采用解三對角矩陣的追趕法！追趕法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程序

2、function x= zhuiganfa( a,b,c,d )%追趕法實現(xiàn)要求：|b1|>|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);for k=n-1:-1:1x(k)=(y(k)-c(k)*x(k+1)/u(k);endendMATLAB命令窗口輸入：a=0

3、 -2 -2 -2 -2 -2 -2 -2;b=2 5 5 5 5 5 5 5;c=-2 -2 -2 -2 -2 -2 -2 0d=220/27 0 0 0 0 0 0 0;x= zhuiganfa(a,b,c,d )運行結(jié)果為：x =8.1478 4.0737 2.0365 1.0175 0.5073 0.2506 0.1194 0.0477存在問題根據(jù)電路分析中的所講到的回路電流法，可以列出8個以回路電流為獨立變量的方程，課本上給出的第八個回路電流方程存在問題，正確的應(yīng)該是；或者可以根據(jù)電路并聯(lián)分流的知識，同樣可以確定。正確的處理結(jié)果應(yīng)為：x =8.1481 4.0741 2.0370 1

4、.0185 0.5093 0.2546 0.1273 0.0637實驗二、線性方程組的迭代方法不同迭代法解線性方程組的對比分別用（1）Jacobi迭代法；（2）Gauss-Seidel迭代法；（3）共軛梯度法解線性方程組迭代初始向量??；(1)Jacobi迭代法3function x,k = jacobi(A,b,x)%x為迭代初始向量e=input('請輸入絕對誤差限:n');D=diag(diag(A);n=size(b);I=eye(n,n);B=I-DA;g=Db;for k=1:50;%最大迭代次數(shù)為50次 x=B*x+g; error=norm(b-A*x);%2-范

5、數(shù)if (error<e)break;endendsprintf('迭代次數(shù):nk= %d',k)sprintf('方程組的解為:n ')xendMatlab命令窗口輸入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);運行結(jié)果請輸入絕對誤差限:10(-3)迭代次數(shù):k= 42ans =方程組的解為:x = 1.0002 -2.0002 2.9997 -1.9999

6、0.9998(2)Gauss-Seidel迭代法function x,k = gau_sei(A,b,x)D=diag(diag(A);L=D-tril(A);U=D-triu(A);e=input('請輸入絕對誤差限:n');for k=1:50; x=(D-L)U*x+(D-L)b;error=norm(b-A*x);if (error<10(-3)break;endendsprintf('迭代次數(shù):nk= %d',k)sprintf('方程組的解為:n ')xend14Matlab命令窗口輸入A=10 1 2 3 4;1 9 -1 2

7、 -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 = gau_sei(A,b,x);運行結(jié)果請輸入絕對誤差限:10(-3)迭代次數(shù):k= 24方程組的解為:x = 1.0002 -2.0002 2.9997 -2.0000 0.9998（3）共軛梯度法function x,k = gongetidu( A,b,x ) e=input('請輸入絕對誤差限:n'); d=b-A*x; r=d; l=(d'*r)./(A*d)'*d); x=x+l*

8、d;for k=1:50 r=b-A*x; t=-(A*d)'*r)./(A*d)'*d); d=r+t*d; l=(d'*r)./(A*d)'*d); x=x+l*d;error=norm(b-A*x);if (error<e)break;endendsprintf('迭代次數(shù):nk= %d',k)sprintf('方程組的解為:n ')xendMatlab命令窗口輸入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

9、 12'x=0 0 0 0 0'x,k = gongetidu(A,b,x );運行結(jié)果請輸入絕對誤差限:10(-3)迭代次數(shù):k= 4方程組的解為:x = 1.0000 -2.0000 3.0000 -2.00001.0000結(jié)果分析在選擇相同的誤差容限時，使用以上三種不同的方法來解方程組，對比迭代次數(shù)上的不同（1）Jacboi迭代法42次（2）Gauss-Seidel24次（3）共軛梯度法4次，可以看出在收斂速度上共軛梯度法最快，Jacobi迭代法最慢，對比方程組的精確解x=1 -2 3 -2 1，可以看出在收斂效果上共軛梯度法最好！實驗三、插值法龍格現(xiàn)象在代數(shù)插值中，為了

10、提高插值多項式對函數(shù)的逼近層度，常常增加節(jié)點個數(shù)，級提高多項式的次數(shù)，但由于高次多項式插值不收斂，會產(chǎn)生Runge現(xiàn)象，本實驗使用函數(shù)為y=11+x2，通過對比Lagrange插值法和分段線性插值法可以觀察到Runge現(xiàn)場的產(chǎn)生與消失。實驗程序函數(shù)fxfunction y = fx( x ) y=1./(1+x*x);endlagrange插值方法function = lagrangechazhi( )n=input('區(qū)間等分份數(shù)：n');s=-5+10/n*0:n;%給定的定點，fx為給定的函數(shù)x=-5:0.01:5;f=0;for q=1:n+1; l=1;%求插值基函數(shù)

11、for k=1:n+1;if k=q; l=l.*(x-s(k)./(s(q)-s(k);else l=l;endend f=f+feval(fx,s(q)*l;%求插值函數(shù)endy=1./(1+x.*x);e=y-f;%誤差subplot(211);plot(x,f,'r',x,y,'k')%作出插值函數(shù)曲線legend('擬合曲線','實際曲線');gridon;subplot(212);plot(x,e,'b');legend('誤差曲線');gridonend區(qū)間等分份數(shù)為10時，如下圖所示

12、，黑色曲線為的實際圖形，紅色曲線即為采用Lagrange插值方法實際描繪出的曲線，從圖中可以觀察到明顯的Runge現(xiàn)象；Runge現(xiàn)象的是因為選擇的插值函數(shù)并不收斂的緣故。分段線性插值方法function = fenduanxianxingchazhi( )clearn=input('區(qū)間等分份數(shù)：n');s=-5+10/n*0:n;%給定的定點，Rf為給定的函數(shù)m=0;for x=-5:0.01:5;ff=0;for k=1:n+1;%求插值基函數(shù)switch kcase 1if x<=s(2); l=(x-s(2)./(s(1)-s(2);else l=0;endca

13、se n+1if x>s(n); l=(x-s(n)./(s(n+1)-s(n);else l=0;endotherwiseif x>=s(k-1)&x<=s(k); l=(x-s(k-1)./(s(k)-s(k-1);elseif x>=s(k)&x<=s(k+1); l=(x-s(k+1)./(s(k)-s(k+1);else l=0;endendend%ff=ff+1./(1+25*s(k)*s(k)*l;%求插值函數(shù)值ff=ff+feval(fx,s(k)*l;%求插值函數(shù)值end m=m+1;f(m)=ff;end%作出曲線x=-5:0.

14、01:5;y=1./(1+x.*x);e=y-f;%誤差subplot(211);plot(x,f,'r',x,y,'k')%作出插值函數(shù)曲線legend('擬合曲線','實際曲線');gridon;subplot(212);plot(x,e,'b');legend('誤差曲線');gridonend區(qū)間等分為20份時，如下圖所示，黑色曲線為的實際圖形，紅色曲線即為采用分段線性插值方法實際描繪出的曲線，從圖中可以觀察到未出現(xiàn)Runge現(xiàn)象，可以觀察誤差曲線可以看到在將區(qū)間等分為20份的時候，和實際曲

15、線的誤差限為0.5X10（-2）；繼續(xù)提高等份份數(shù)，可以繼續(xù)降低誤差！實驗四、數(shù)值積分考紐螺線考紐螺線的形狀像鐘表的發(fā)條，也稱回旋曲線，它在直角坐標(biāo)系中的參數(shù)方程為曲線關(guān)于原點對稱。取a=1，參數(shù)s的變化范圍-5，5，容許誤差限分別是10-6和10-10。選取適當(dāng)?shù)墓?jié)點個數(shù)，利用數(shù)值積分方法計算曲線上點的坐標(biāo)，并畫出曲線的圖形。實驗程序（程序中編寫的積分方法根據(jù)數(shù)值積分中常見的Romberg求積公式得到）function X,Y = kaoniuluoxian( )%命令窗口輸入X,Y = kaoniuluoxian( )formatlong; T=zeros(20,20); a=input(

16、'積分下限：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:1000 k=1; h=c(m)-a; T(1,1)=h/2*(fx1(a)+fx1(c(m);%計算積分上限為c(m)時fx1的積分值fori=1:

17、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:k T(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;end k=k+1; h=h/2;endX(m)=T(1,k+1);endfor m=1:1000 k=1; h=c(m)-a; T(1,1)=h/2*(feval(fx2,a)+feval(fx2,c(m);%feval用來求函數(shù)fx在給定點處的值fori=1:20 %計算積分上

18、限為c(m)時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:k T(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;end k=k+1; h=h/2;endY(m)=T(1,k+1);end X=-X; Y=-Y;plot(X,Y,'k-',-X,-Y,'k-');gridon;end命令窗口輸入X,Y = kaoniuluoxian( );運行結(jié)果

19、積分下限：0積分上限：5允許誤差限：10(-6)繪出圖形如下圖所示：容許誤差限為10(-2)時，所繪出圖形如下圖：結(jié)果分析對比上述兩圖，采用Romberg積分方法容許誤差限分別為10（-6）和10（-2）時，所繪出的圖形形狀像鐘表的發(fā)條，即為考紐螺線的形狀，但10(-2)對應(yīng)的曲線誤差偏大！實驗五、同軸矩形金屬波導(dǎo)中靜電場的分布問題問題：同軸矩形金屬波導(dǎo)橫截面結(jié)構(gòu)如下圖所示，其中內(nèi)金屬壁電位為100V，外金屬壁為0V，作出其橫截面內(nèi)電位線分布圖；問題背景：在一個二維矩形區(qū)域D內(nèi)，電位函數(shù)滿足二維拉普拉斯方程在邊界L上各點的電位置已給定，即。為了采用差分法求解D內(nèi)的電位函數(shù)，作平行于坐標(biāo)軸的兩組

20、平行線;將區(qū)域D劃分為許多正方形網(wǎng)絡(luò)，網(wǎng)格的邊長h稱為步長，兩組平行線的交點成為網(wǎng)格的節(jié)點。用表示節(jié)點處的電位值，利用二元函數(shù)的Taylor展開式，可將與相鄰的節(jié)點上的電位函數(shù)值表示為式中表示量級的微量，將以上四式代入拉普拉斯方程可得這樣在處的電位值滿足的拉普拉斯方程近似的就可以以一個差分方程來代替。然后以此計算每點的電位，即利用上述近似拉普拉斯方程來計算，用圍繞它的四個點的電位平均值作為它的電位值，當(dāng)所有的點計算完后，用新值代替舊值就完成了一次迭代運算，然后再進行下一次迭代運算，知道計算的新值和舊值之差小于指定的誤差容限為止。就為簡單迭代法！Matlab程序function = dianwe

21、ifenbu( )FI=100 ;hx=23 ;hy=23; %外金屬槽橫截面沿x、y方向網(wǎng)格數(shù)目，步進長度為1CX1=9;CX2=15;CY1=9;CY2=15; %內(nèi)金屬槽占據(jù)網(wǎng)格位置 : x 輸 9 -15 ，y 軸 9-15(剛好在正中間%以上數(shù)據(jù)可以按要求更改DX=1+CX2-CX1 ;DY=1+CY2-CY1;v=ones(hy,hx)*50;%每個網(wǎng)格點迭代初值mesh1=ones(hy,hx)*2;v(1,:)=zeros(1,hx);%上邊界初值mesh1(1,:)=zeros(1,hx);v(hy,:)=zeros(1,hx);%下邊界初值mesh1(hy,:)=zeros

22、(1,hx);v(:,1)=zeros(hy,1);%左邊界初值mesh1(:,1)=zeros(hy,1);v(:,hx)=zeros(hy,1);%右邊界初值mesh1(:,hx)=zeros(hy,1);v(CY1:CY2,CX1:CX2)=ones(DX,DY)*FI;%內(nèi)邊界初值mesh1(CY1:CY2,CX1:CX2)=ones(DX,DY);k=0;%記錄迭代次數(shù)difmax=1.0;while(difmax>1e-6)%誤差限 k=k+1;difmax=0.0;fori=2:hy-1for j=2:hx-1 m=mesh1(i,j);if (m>=2)oriv=v

23、(i,j);%原電勢值 v(i,j)=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 (dif>difmax) difmax=dif;%取誤差的最大值endendendendendsubplot(1,2,1);%figure(1);mesh(v);%畫三維曲面圖axis(-2,hx+3,-2,hy+3,0,100);title('等勢線三維分布圖');subplot(1,2,2);%figure(2);contour(v,13);%畫等電

24、位線圖holdon;x=1:1:hx;y=1:1:hy;xx,yy=meshgrid(x,y);%形成柵格Gx,Gy=gradient(v,0.5,0.5);%計算梯度quiver(xx,yy,Gx,Gy,'r');%根據(jù)梯度畫箭頭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('橫截面等勢線分布圖');text(CX1+0.6,CY

25、1+(CY2-CY1)/2,'U=100V','fontsize',10);text(hx/2,hy+1,'U=0V','fontsize',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;end命令窗口輸入dianweifenbu();運行結(jié)果：從上圖可以看出，在靠近內(nèi)金屬波導(dǎo)時，等勢線分布更加密集，同時電力線垂直于金屬波導(dǎo)邊緣，符合實際同軸金屬波導(dǎo)中的場強分布規(guī)律！總結(jié)實驗報告中涉及到的內(nèi)容主要包括數(shù)值計算方法課中所講