浙大計算方法上機報告

1、學號：3100300038 姓名： 專業(yè)：作業(yè)1 :用列主元高斯消去法和列主元三角分解法解p227頁第3題1.列主元高斯消去法目的：用高斯消去法解ax=b時，其中設a為非奇異矩陣，可能出現(xiàn)akk =0情況，這時 必須進行帶行交換的高斯消去法。但在實際計算中即使akk) 0但其絕對值很小時，用ak?作除數(shù)，會導致中間結(jié)果 矩陣a(k)元素數(shù)量級嚴重增長和舍入誤差的擴散，使得最后的計算結(jié)果不可靠。列主元高斯消去法可以難過一般高斯法的這些缺點。、列主元高斯消去法解方程的 matlab程序如下:function a=columneli(a)%對矩陣a進行列主元消去m n=size(a);%求取a的行數(shù)

2、m和列數(shù)nfor i=1:m-1maxele,pos=max(abs(a(i:end,i);maxrow=pos+i-1;%在每次變換前尋找絕對值最大主所在列maxrowif a(maxrow,i)=0 disp(矩陣為奇異矩陣)return%對于非奇異矩陣，在程序中給予提示，結(jié)束程序end if(maxrow=i) temp=a(maxrow,:);a(maxrow,:)=a(i,:);a(i,:)=temp; end%與列主元絕對值最大的行進行行交換for j=i+1:m a(j,i)=a(j,i)/a(i,i);%求取第 j 歹u主元for k=i+1:n a(j,k)=a(j,k)-a

3、(j,i)*a(i,k); %對第j列主元進行行變換end end end function x=elisolve(a) %利用列主元消去的結(jié)果求方程的解，a為方程組的增廣矩陣a=columneli(a);im,n=size(a); x=zeros(m,1); for i=m:-1:1x(i)=(a(i,n)-a(i,i:m)*x(i:m)/a(i,i);end二、列主元高斯消去法解p227頁第3題:034pd -11-11 x2 = 2-212 j?3 j 口答案：x=7/6 -1/3 1/2 t三、程序流程圖如下:條件圖框里面沒有條件式，因為i是從m到1所以i=1后運行下面的命令，直接輸出

4、結(jié)果。這里也是一樣的，直接運行下面的命令。四、matlab中的命令截圖如下：command vindow? x elisolve(0 3 4 l;l 7 1 2;2 1 2 3)1. 1667-0, 33330.5000aiis7/5-1/31/2運行結(jié)果：x=7/6 -1/3 1/2上機結(jié)果和正確答案完全相同。2.列主元三角分解法一，列主元三角分解法解方程的matlab程序如下：function l,u,p=tridecom(a) %對矩陣 a 進行三角分解m n=size(a); % 求取a的行數(shù) m和列數(shù)nif(m=n) |disp(矩陣的行數(shù)與列數(shù)不相等，無法進行三角分解)return

5、endl=zeros(n);u=zeros(n);p=eye(n);%p*a=l*umaxele,goodrow=max(abs(a(:,1);temp=a(1,:);a(1,:)=a(goodrow,:);a(goodrow,:尸temp;temp=p(1,:);p(1,:)=p(goodrow,:);p(goodrow,:尸temp; %行交換u(1,:)=a(1,:);l(1:end,1)=a(1:end,1)/u(1,1);for r=2:ntempurr=zeros(n,1); for k=r:ntempurr(k尸a(k,r)-l(k,1:r-1)*u(1:r-1,r); endm

6、axele,goodrow=max(abs(tempurr);|temp=a(r,:);a(r,:)=a(goodrow,:);a(goodrow,:)=temp;temp=l(r,:);l(r,:)=l(goodrow,:);l(goodrow,:尸temp;temp=p(r,:);p(r,:尸p(goodrow,:);p(goodrow,:尸temp; for i=r:nu(r,i尸a(r,i)-l(r,1:r-1)*u(1:r-1,i);endfor i=r:nl(i,r)=(a(i,r)-l(i,1:r-1)*u(1:r-1,r)/u(r,r); endendfunction x=lu

7、solve(a) %用三角分解法解線性方程組ax=b , a為方程組的增廣矩陣a|bm,n=size(a);l,u,p=tridecom(a(:,1:m); %對 a進行三角分解 a=p*l*ux=p*a(:,n);%對a的最后一列進行與方程組系數(shù)三角分解時相同的行變換for i=2:mx(i)=x(i)-l(i,1:i-1)*x(1:i-1); %求出(la-1)*p*bend x(m)=x(m)/u(m,m); for i=m-1:-1:1 x(i)=(x(i)-u(i,i+1:m)*x(i+1:m)/u(i,i);%求出方程組的解 x=(ua-1)*(la-1)*p*bend、程序流程圖

8、如下:ihritelt?. yuajb門“i -hum (dkm (博(i llj, u.p p =t r i dee om ( 0lu 00000001.dood0q.5。0。=q_ 5 uqq i, qqqu0011 2 lusolve(0ans =command vindow3 i 1;1 -1 1 2;2 1 2 3)commcind vindow p*0 3 4;1 1 2 ans -1. 1667-0, 33330. 50。 sym ans =7/6-1/31/2可見，p*a=l*u作業(yè)2;編寫matlab程序，用最小二乘法擬合多項式曲線(解p89頁第1題)目的和意義；根據(jù)觀測或?qū)?/p>

9、驗得到一系列的數(shù)據(jù)，確定了與自變量的某些點相應的函數(shù)值。當函數(shù)比較復雜或根本無法寫出解析式時，往往根據(jù)觀測數(shù)據(jù)構(gòu)造一個 適當?shù)暮唵蔚暮瘮?shù)近似地代替要尋求的函數(shù)。利用最小二乘法擬合多項式來描述其函數(shù)。一、用最小二乘法擬合多項式曲線的matlab程序如下：function a=polyercheng(x,y,n)x=x;y=y;m=size(x,1);c=ones(m,1+n);for i=1:nc(:,i+1)=x.ai;enda=c*c;b=c*y;a=(ab);a=fliplr(a);c=polyval(a,x)-y;d=max(abs(c);fprintf(最大偏差：%fn,d);d=sq

10、rt(c*c);fprintf(均方誤差：f,d);b=linspace(min(x),max(x),500);c=polyval(a,b);plot(x,y,g*,b,c,r);end、程序的流程圖如下:輸出a 結(jié)束程序a=polyercheng(x,y,n)開始輸入一組數(shù)據(jù) (x,y)和多項式次數(shù)n結(jié)果：一次擬合多項式最大偏差：2.700000均方誤差：3.271085擬合多項式曲線：y=6.5500x-12.5000二次擬合多項式最大偏差：1.950000均方誤差：2.906888擬合多項式曲線：y=0.1875x2+4.6750x-8.7500三、最小二乘法擬合多項式曲線解 p89頁第

11、1題: 已知實驗數(shù)據(jù)如下：x2468y2112840用最小二乘法求一次和二次擬合多項式，分別算出均方誤差和最大偏差。四、解p89頁第1題在matlab中的截圖如下: a=poly?rcheng x,yy 1) 最大偏差:2. 700000均方誤差；3.2710b5a -6.5500 -12.5000y=6.5500x-12.5000 a=p iyercheng 2 最大偏差：1.950000均方誤差；2,9068380.18754.6750-8. 7500y=0.1875x2+4.6750x-8.75008其中a為擬合多項式的系數(shù)陣列。通過最小二乘法擬合多項式得到的多項式接近與實驗值(代入后)

12、作業(yè)3:用龍貝格算法計算，1但幽dx使截斷誤差不超過0.5e-6 o0 x目的和意義；在一元函數(shù)的積分學中，我們已經(jīng)熟知，若函數(shù)f(x)在區(qū)間a, b 上連b續(xù)且其原函數(shù)為 f(x),則可用牛頓一萊布尼茲公f f(x)dx= f(a) - f(b) ,a來求定積分。牛頓萊布尼茲公式雖然在理論上或在解決實際問題中都起了很大的作用，但它并不能完全解決定積分的計算問題。當被積函數(shù)f(x)在區(qū)間a,b 上連續(xù)時，要使得復合梯形公式比較精確地代替定積分b f (x)dx = f(a)-f(b)可將分點(即基點)加密，也就是將區(qū)間a,b 細分，然后利用復合梯形公式求積將收斂緩慢的梯形值序列加工成收斂迅速的

13、sn,cn,rn ,這種加速方法稱為龍貝格算法。它可以把復雜的積分比較精確的求出來。，龍貝格算法matlab程序如下：function z=romberg(ff,a,b,epsi,k) %用龍貝格算法計算積分if nargin=4k=16;end%如果只有4個傳入?yún)?shù)，則將二分次數(shù)k設置為默認值16t=zeros(1,k+1);format long;function y=f(x) %嵌套函數(shù)，用來判斷積分區(qū)域內(nèi)的奇點y=ff(x); %若 ff 參數(shù)為 (x)sin(x)/x, 則有 y=sin(x)/xif isnan(y)=true|isinf(y尸truedisp(請確保所給函數(shù)在積分

14、區(qū)域內(nèi)無奇點)；endendt(1)=(b-a)/2*(f(a)+f(b);for l=1:ksum=0;for i=1:(2a(l-1)sum=sum+f(a+(2*i-1)*(b-a)/(2n);endt(l+1)=t(l)/2+(b-a)/(2al)*sum;end%計算梯形序列for m=1:ktemp=t;for i=m+1:k+1%將新序列放入 t(m+1)-t(k+1) 中t(i)=(4am)*temp(i)-temp(i-1)/(4am-1)；endif abs(t(m+1)-t(m)epsi|break;endendz=t(m+1);fprintf(截斷誤差為 dn,t(m+

15、1)-t(m);end、用龍貝格算法計算p120頁第5題：1 cin()計算sn3dx,使截斷誤差不超過0.5e-60 x結(jié)果：截斷誤差為6.632355e-008=0.946083070387222、程序流程圖如下嵌套函數(shù)y=f(x)nargin=4no識別積分區(qū)域內(nèi)奇點for l=1:k1for i=1:(2a(l-1)disp(請確保所給函數(shù) 在積分區(qū)域內(nèi)無奇點yesdisp(t()=(4am)*temp(i)-temp(i-1)/(4am-l)；if abs(t(m+1)-t(m) romberg(k sin x /a, eps, 1，0. 5e-6j 截斷誤差為5. 532355e-008arts -0. 946083070337222指定二分次數(shù)為10次時：cornmand 1 vindow， 口 x romberg (x) sincxj/xj eps, i, 0. be