數(shù)值分析報(bào)告資料

1、也為fl乂扎熔JianpXi University Of Scinc? And Technology數(shù)值分析課程名稱數(shù)值分析實(shí)驗(yàn)報(bào)告專業(yè)班級應(yīng)用數(shù)學(xué)111班 姓名陳偉學(xué)號8_教學(xué)教師岳雪芝實(shí)驗(yàn)二 函數(shù)逼近與曲線擬合、問題提出從隨機(jī)的數(shù)據(jù)中找出其規(guī)律性，給出其近似表達(dá)式的問題，在生產(chǎn)實(shí)踐和科學(xué)實(shí)驗(yàn) 中大量存在，通常利用數(shù)據(jù)的最小二乘法求得擬合曲線。在某冶煉過程中，根據(jù)統(tǒng)計(jì)數(shù)據(jù)的含碳量與時(shí)間關(guān)系，試求含碳量與時(shí)間t的擬合曲線。t（分）0 510152025303540455055y(x10)0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.

2、64二、要求1用最小二乘法進(jìn)行曲線擬合；2、 近似解析表達(dá)式為：:：(t)= a)t a2t a3t3 ；3、 打印出擬合函數(shù) (t)，并打印出(tj)與y(tj)的誤差，j =1,2,|1,12 ；4、另外選取一個(gè)近似表達(dá)式，嘗試擬合效果的比較；5、*繪制出曲線擬合圖。三、目的和意義1掌握曲線擬合的最小二乘法；2、最小二乘法亦可用于解超定線代數(shù)方程組；3、 探索擬合函數(shù)的選擇與擬合精度間的關(guān)系。四、實(shí)驗(yàn)步驟：輸入代碼：t=0 5 10 15 20 25 30 35 40 45 50 55;y=0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.0

3、2 4.64;for i=1:12A(i,1)=t(i);A(i,2)=t(iF2;A(i,3)=t(i).A3;endp_star=Ay：y_star=p_star(1)*t+p_star(2)*t.A2+p_star(3)*t.A3;s_star=0;for i=1:12s_star=s_star+(y_star(i)-y(i)A2;endplot(t,y,*,t,y_star,g);比較擬和結(jié)果，分別取二次，三次和四次多項(xiàng)式，再做最小二乘法，代碼如下：t=0 5 10 15 20 25 30 35 40 45 50 55;y=0 1.27 2.16 2.86 3.44 3.87 4.15

4、 4.37 4.51 4.58 4.02 4.64;%3次擬合p3=polyfit(t,y,3)y_n ew=polyval(p3,t);s_3=0;for i=1:12s_3=s_3+(y_ new(i)-y(i)A2;end%2次擬合p2=polyfit(t,y,2);y_n ew2=polyval(p2,t);s_2=0;for i=1:12s_2=s_2+(y_ new2(i)-y(i)A2;end%4次擬合p4=polyfit(t,y,4);y_n ew4=polyval(p4,t);s_4=0;for i=1:12s_4=s_4+(y_ new4(i)-y(i)A2;end比較4種

5、擬合函數(shù)，結(jié)論：常數(shù)項(xiàng)為0的三次擬合函數(shù)在保證擬合精度的同時(shí)，保證0點(diǎn)的函數(shù)值仍為0,是四種擬合方法中最好的一種。0102030405060-4原圖與曲線擬合圖如下所示：6543210-1分析：從上面的擬合中也可以知道多項(xiàng)式擬合誤差平方和隨著擬合多項(xiàng)式次數(shù)的增加而 逐漸減少，擬合曲線更靠近實(shí)際數(shù)據(jù)，擬合更準(zhǔn)確。實(shí)驗(yàn)三數(shù)值積分與數(shù)值微分基本題Romberg算法，計(jì)算(1)10選用復(fù)合梯形公式，復(fù)合Simpson公式,sin xdx( f (0) =1,10.9460831)x、要求1、編制數(shù)值積分算法的程序；2、分別用兩種算法計(jì)算同一個(gè)積分，并比較其結(jié)果；3、分別取不同步長h=(b-a)/n，試

6、比較計(jì)算結(jié)果(如 n = 10, 20等)；4、給定精度要求,試用變步長算法，確定最佳步長。三、目的和意義1、深刻認(rèn)識數(shù)值積分法的意義；2、明確數(shù)值積分精度與步長的關(guān)系；3、根據(jù)定積分的計(jì)算方法，結(jié)合專業(yè)考慮給出一個(gè)二重積分的計(jì)算問題。四、實(shí)驗(yàn)步驟：代碼 1，復(fù)合梯形公式function result=trapezoid_integration(f,count_node,a,b)%composite trapezoid integration,f,函數(shù)表達(dá)式，count_node 插入節(jié)點(diǎn)數(shù)量，a,開始點(diǎn),b 結(jié)束 點(diǎn)h=(b-a)/count_node;%h 步長np1=count_node

7、+1;% 包括端點(diǎn) ,總共的點(diǎn)的數(shù)量x=a:h:b;%x 自變量c=ones(np1,1);c(1,1)=0.5;c(np1,1)=0.5;% 兩個(gè)端點(diǎn)都取一半syms symbol_x;fx=subs(f,symbol_x,x);result=vpa(h*fx*c,9);代碼 2，復(fù)合 simpson 公式function result=simpson_integration(f,count_node,a,b)%composite simps on in tegrati on,f,函數(shù)表達(dá)式 ,co unt_node 插入節(jié)點(diǎn)數(shù)量，a,開始點(diǎn),b結(jié)束點(diǎn) h=(b-a)/count_node;

8、%h 步長np1=2*count_node+1;% 包括端點(diǎn) ,總共的 2n+1 點(diǎn)x=a:h/2:b;%x 自變量c=ones(np1,1);for i=1:np1if(mod(i,2)=0)% 偶數(shù)c(i,1)=2/3;else%奇數(shù)c(i,1)=1/3;endendc(1,1)=1/6;c(np1,1)=1/6;% 兩個(gè)端點(diǎn)都取 1/6syms symbol_x;fx=subs(f,symbol_x,x);result=vpa(h*fx*c,9);代碼 3，檢驗(yàn)復(fù)合梯形公式和復(fù)合simpson 公式function result1=test_integration()%測試復(fù)合梯形公式，

9、復(fù)合simpson 公式 ,用到 trapezoid_integration,simpson_integrationsyms symbol_x;f1=sin(symbol_x)/symbol_x;syms A1;% 用于輸出 f1 積分結(jié)果for i=1:9A1(i,1)=trapezoid_integration(f1,i,1e-9,1);A1(i,2)=simpson_integration(f1,i,1e-9,1);endk=10;for i=10:10:100A1(k,1)=trapezoid_integration(f1,i,1e-9,1);A1(k,2)=simpson_i nte

10、gratio n(f1,i,1e-9,1); k=k+1;endIS(2)=i nt(f1,symbol_x,0,1);vpa(IS(2)-result仁 A1;結(jié)果分析從上面的計(jì)算結(jié)果可以看出復(fù)合梯形公式和復(fù)合simps on公式的穩(wěn)定性，并且步長越小精度越高，當(dāng)n趨于正無窮，即步長趨于0時(shí)，上述兩公式的計(jì)算值都收斂到積分值。從收斂的速度來看，復(fù)合 simps on公式優(yōu)于復(fù)合梯形公式。實(shí)驗(yàn)四 線方程組的直接解法一、問題提出：1、三對角形線性方程組4-1-140-10000000000000X1 17 150X20-14-1000000X3-1300-14-100000X42000-14-1

11、0000X60000-14-1000X-1200000-14-100X714000000-14-10X-40000000-14-1X95.00000000-14.x10-1-5 一*x= (2, 13, 01,2,3：0,1,1)二、要求1、 對上述三個(gè)方程組分別利用Gauss順序消去法與 Gauss列主元消去法；平方根 法與改進(jìn)平方根法；追趕法求解（選擇其一） ；2、應(yīng)用結(jié)構(gòu)程序設(shè)計(jì)編出通用程序；3、比較計(jì)算結(jié)果，分析數(shù)值解誤差的原因；4、盡可能利用相應(yīng)模塊輸出系數(shù)矩陣的三角分解式。三、目的和意義1、通過該課題的實(shí)驗(yàn)，體會模塊化結(jié)構(gòu)程序設(shè)計(jì)方法的優(yōu)點(diǎn)；2、運(yùn)用所學(xué)的計(jì)算方法，解決各類線性方程

12、組的直接算法；3、提高分析和解決問題的能力，做到學(xué)以致用；通過三對角形線性方程組的解法，體會稀疏線性方程組解法的特點(diǎn)。四、實(shí)驗(yàn)步驟:高斯順序消去法function result=Gauss(A,d)%高斯順序消元法n,m=size(A) k=ones(1,n);% 用于記錄 -A(j,i)/A(i,i) for i=1:nif(A(i,i)=0)主元為 0，不能使用高斯順序消元法 return;endfor j=i+1:nk(j)=-A(j,i)/A(i,i)A(j,:)=A(i,:).*k(j)+A(j,:);% 加到第 j 行 d(j)=d(i).*k(j)+d(j);endendx=ze

13、ros(n,1)x(n)=d(n)/A(n,n);for i=n-1:-1:1x(i)=(d(i)-A(i,:)*x)/A(i,i);endresult=xGauss 列主元消去法function result=Gauss2(A,d)%高斯列主元消元法n,m=size(A); d_length,d_width=size(d);if(d_length=n | d_width=1) 方程右端向量輸入有誤 return;endif(n=m) 系數(shù)矩陣不是方陣 return;endAandD=A d;% 增廣陣k=ones(1,n);% 用于記錄 -A(j,i)/A(i,i) record=1:n;%

14、 用于記錄行的交換情況 for i=1:n%選取第 i 列的最大元進(jìn)行換位 l,index=max(abs(AandD(i:n,i); if(index=i)temp_record=record(index); record(index)=record(i+index-1); record(i+index-1)=temp_record; temp=AandD(i,:);AandD(i,:)=AandD(i+index-1,:);加到第 j 行AandD(i+index-1,:)=temp; end %順序消元 if(AandD(i,i)=0) 主元為 0 ，不能使用高斯列主元順序消元法 ret

15、urn; end for j=i+1:n k(j)=-AandD(j,i)/AandD(i,i); AandD(j,:)=AandD(i,:).*k(j)+AandD(j,:);% end endA=AandD(1:n,1:n);d=AandD(1:n,n+1);x=zeros(n,1);x(n)=d(n)/A(n,n);for i=n-1:-1:1 x(i)=(d(i)-A(i,:)*x)/A(i,i); end %由于不是順序消元，所以還要再換回來 for i=1:n result(i)=x(record(i);end追趕法function result=chase_after(A,d)n

16、,m=size(A);d_length,d_width=size(d); if(d_length=n | d_width=1) 方程右端向量輸入有誤 return;endif(n=m)系數(shù)矩陣不是方陣 return;endfor i=2:ntoCheck1=diag(A,i);toCheck2=diag(A,-i); if(norm(toCheck1)=0|norm(toCheck2)=0) 系數(shù)矩陣不是三對角陣 return;endendb=diag(A);a=zeros(n,1);a(2:n)=diag(A,-1);c=diag(A,1);if(b(1)=0) 因順序主子式為 0, 不能進(jìn)

17、行 Crout 分解 return;endl(1)=b(1);y(1)=d(1)/l(1);u(1)=c(1)/l(1);for i=2:nl(i)=b(i)-a(i)*u(i-1);if(l(i)=0) 因順序主子式為 0, 不能進(jìn)行 Crout 分解 return;endy(i)=(d(i)-y(i-1)*a(i)/l(i);if(i=n)u(i)=c(i)/l(i);endendx(n)=y(n);for i=n-1:-1:1x(i)=y(i)-u(i)*x(i+1);endy=y;result=x;檢驗(yàn)三種直接解線性方程組的方法function test_chase_afterA= 4

18、 -1 0 0 0 0 0 0 0 0;-1 4 -1 0 0 0 0 0 0 0;0 -1 4 -1 0 0 0 0 0 0;0 0 -1 4 -1 0 0 0 0 0;0 0 0 -1 4 -1 0 0 0 0;0 0 0 0 -1 4 -1 0 0 0;0 0 0 0 0 -1 4 -1 0 0;0 0 0 0 0 0 -1 4 -1 0;0 0 0 0 0 0 0 -1 4 -1;0 0 0 0 0 0 0 0 -1 4;d=7 5 -13 2 6 -12 14 -4 5 -5;% 按列存儲L,U = lu(A)%檢驗(yàn)高斯順序消去法 result_Gauss=Gauss(A,d) %檢

19、驗(yàn)高斯列主元消去法 result_Gauss2=Gauss2(A,d) %檢驗(yàn)追趕法 result_chase=chase_after(A,d)運(yùn)行結(jié)果：高斯順序消去法的計(jì)算結(jié)果：result_Gauss =21-31.7849e-0161-23-1.4874e-0161-1高斯列主元消去法的計(jì)算結(jié)果：result_Gauss =21-31.7849e-0161-23-1.4874e-0161-1追趕法的計(jì)算結(jié)果：result_Gauss =21-31.7849e-0161-23-1.4874e-0161-1結(jié)果分析：高斯順序消去法與精確解x* =(2,1,-3,0,1,-2,3,0,1,-1

20、)丁，比較其誤差是 1.7849*10(16),該誤差可以忽略，事實(shí)上這是由計(jì)算機(jī)的除法引起的。而在這個(gè)問題中，因?yàn)轫?序消元的除法中用到的數(shù) k的絕對值都小于1,這就避免使用絕對值較大的數(shù)，使得總體舍入誤差得到有效控制，所以高斯順序消元法的效果也很好。同樣高斯列主元消去法與精確解的誤差都是(1 )式的結(jié)果，說明在這個(gè)問題中誤差都控制得很好。實(shí)驗(yàn)五解線性方程組的迭代法Jacobi 迭代法，Gauss-Seidel 迭代、問題提出對實(shí)驗(yàn)四所列目的和意義的線性方程組，試分別選用 法和SOR方法計(jì)算其解。二、要求1、體會迭代法求解線性方程組，并能與消去法做以比較；2、 分別對不同精度要求，如；=10

21、,10二10,由迭代次數(shù)體會該迭代法的收斂快慢；3、對方程組2, 3使用SOR方法時(shí)，選取松弛因子 3 =0.8 , 0.9 , 1, 1.1 , 1.2等，試看對算法 收斂性的影響，并能找出你所選用的松弛因子的最佳者；4、給出各種算法的設(shè)計(jì)程序和計(jì)算結(jié)果。三、目的和意義1、通過上機(jī)計(jì)算體會迭代法求解線性方程組的特點(diǎn)，并能和消去法比較；2、運(yùn)用所學(xué)的迭代法算法，解決各類線性方程組，編出算法程序；3、 體會上機(jī)計(jì)算時(shí)，終止步驟|x(k xk：；或k (給予的迭代次數(shù))，對迭代法斂散性的意義；4、體會初始解x0，松弛因子的選取，對計(jì)算結(jié)果的影響。四、實(shí)驗(yàn)步驟：Jacobi迭代法function r

22、esult=Jacobi(A,b,x0,yupsil ong)%Jacobi迭代方法,x0是迭代用的初始向量，yupsilong是解的誤差上限n ,m=size(A);b_le ngth,b_width=size(b);if(b_le ngth=n | b_width=1)方程右端向量輸入有誤return;endif(n=m)系數(shù)矩陣不是方陣return;endif(det(A)=0)系數(shù)矩陣是奇異陣returnend%分解A=D-L-U;D是對角元，L是對角元為0的下三角陣，U是對角元為0的上三角陣L=tril(A,-1);U=triu(A,1);D=A-L-U;L=-L;U=-U;%計(jì)算迭

23、代矩陣B0=DA(-1);B=B0*(L+U);g=B0*b;x=x0;delta=1;coun t_iterati on=0;while delta=yupsil ongif(cou nt_iterati on 1000) 超過迭代次數(shù)上限，認(rèn)為算法此時(shí)不收斂 return;end x1=B*x+g;delta=norm(x1-x);x=x1; count_iteration=count_iteration+1;end count_iteration result=x;Gauss_Seidol 迭代法function result=Gauss_Seidol(A,b,x0,yupsilong)

24、 %Gauss_Seidol 方法 ,x0 是迭代用的初始向量 ,yupsilong 是解的誤差上限 n,m=size(A);b_length,b_width=size(b); if(b_length=n | b_width=1) 方程右端向量輸入有誤 return;end if(n=m) 系數(shù)矩陣不是方陣 return;endif(det(A)=0) 系數(shù)矩陣是奇異陣 returnend%分解A=D-L-U;D是對角元，L是對角元為0的下三角陣，U是對角元為0的上三角陣 L=tril(A,-1);U=triu(A,1);D=A-L-U;L=-L;U=-U;%計(jì)算迭代矩陣B0=(D-LF(-1

25、); B=B0*U; g=B0*b;x=x0;delta=1; count_iteration=0;while delta=yupsilong if(count_iteration1000) 超過迭代次數(shù)上限，認(rèn)為算法此時(shí)不收斂 return;end x1=B*x+g;delta=norm(x1-x);x=x1; count_iteration=count_iteration+1;endcount_iterationresult=x;SOR 迭代法function result=SOR(A,b,x0,yupsilong,omiga)%S0迭代方法，xO是迭代用的初始向量，yupsilong是解

26、的誤差上限，omiga是松弛因子n ,m=size(A);b_length，b_width=size(b); if(b_length=n | b_width=1) 方程右端向量輸入有誤 return;endif(n=m) 系數(shù)矩陣不是方陣 return;endif(det(A)=0) 系數(shù)矩陣是奇異陣 returnend%分解A=D-L-U;D是對角元，L是對角元為0的下三角陣，U是對角元為0的上三角陣L=tril(A，-1);U=triu(A，1);D=A-L-U;L=-L;U=-U;%計(jì)算迭代矩陣B0=(D-omiga*L)A(-1);B=B0*(1-omiga)*D+omiga*U);

27、g=omiga*B0*b;x=x0;delta=1;count_iteration=0;while delta=yupsilong if(count_iteration1000) 超過迭代次數(shù)上限，認(rèn)為算法此時(shí)不收斂 return;endx1=B*x+g;delta=norm(x1-x);x=x1;count_iteration=count_iteration+1;endcount_iterationresult=x;檢驗(yàn) Jocobi,GS,SOR 迭代方法function test_Gauss_Seidol clc;A= 4 -1 0 0 0 0 0 0 0 0; -1 4 -1 0 0

28、0 0 0 0 0; 0 -1 4 -1 0 0 0 0 0 0; 0 0 -1 4 -1 0 0 0 0 0; 0 0 0 -1 4 -1 0 0 0 0; 0 0 0 0 -1 4 -1 0 0 0; 0 0 0 0 0 -1 4 -1 0 0; 0 0 0 0 0 0 -1 4 -1 0; 0 0 0 0 0 0 0 -1 4 -1; 0 0 0 0 0 0 0 0 -1 4; b=7 5 -13 2 6 -12 14 -4 5 -5; n,m=size(A);yupsilong=1e-3 1e-4 1e-5 1e-6 1e-7 1e-8 1e-9 1e-10; x_real=2 1 -

29、3 0 1 -2 3 0 1 -1;% 精確解 測試初始點(diǎn)對迭代的影響for i=1:8 xO=o nes( n,1)*10F-1);x0 result_Jacobi=(Jacobi(A,b,x0,yupsilong(3); if(norm(result_Jacobi-x_real)yupsilong(3) 誤差過大 end result_GS=(Gauss_Seidol(A,b,x0,yupsilong(3); if(norm(result_GS-x_real)yupsilong(3) 誤差過大 end result_SOR=(SOR(A,b,x0,yupsilong(3),1.1); if

30、(norm(result_GS-x_real)yupsilong(3) 誤差過大 end end 測試 yupsilong 對迭代的影響 x0=ones(n,1);for i=1:length(yupsilong) (yupsilong(i) result_Jacobi=(Jacobi(A,b,x0,yupsilong(i); result_GS=(Gauss_Seidol(A,b,x0,yupsilong(i); result_SOR=(SOR(A,b,x0,yupsilong(i),1.1); end測試omiga松弛因子對迭代的影響 omiga=0.1:0.1:2 for i=1:len

31、gth(omiga) omiga(i) result_SOR=(SOR(A,b,x0,yupsilong(3),omiga(i); end運(yùn)行結(jié)果 初步計(jì)算結(jié)果1) Jocobi 迭代方法，初始點(diǎn)(1,1,1,1, 1,1, 1,1,1,1), , f10a(-6),迭代 22 次：Result_Jocobi=(2 1 -3 1.6191e-007 1 -2 3 1.4127e-007 1 -1)2) Gauss-Seidol 迭代方法，初始點(diǎn)(1, 1, 1, 1, 1, 1, 1, 1, 1, 1), 日0八(-6),迭代 13 次： Result_GS=(2 1 -3 8.0741e-008 1 -2 3 4.3362e-0