數(shù)值計算技巧Matlab實題訓練(內(nèi)附程序,模型)

1、數(shù)值計算方法訓練實習報告題 目：6 A組院 系：上海電力學院數(shù)理學院專業(yè)年級：信息與計算科學專業(yè) 2009級學生姓名：XX遠學號： 200924262011年7月8日第1題：含炭量與時間的關(guān)系在某冶煉過程中，鋼的含炭量y與時間t的統(tǒng)計數(shù)據(jù)如下t 0510152025303540455055y 01.22.12.83.43.84.14.34.54.54.04.676647571824（1）畫出原始數(shù)據(jù)分布趨勢圖；（2） 用最小二乘法求鋼的含炭量y與時間t的擬合曲線y二at bt2 ct3 ；（3）打印出擬合曲線；（4）另外選用y=atb進行擬合，比較二種擬合的效果。解：分析：使用到曲線擬合的最小

2、二乘法，對于擬合函數(shù)，盡量轉(zhuǎn)化為可以方便 提煉出基函數(shù)的方程。在明確基函數(shù)的基礎(chǔ)上，通過計算，得到各個系數(shù)，得到 法方程組（1），程序：function yuan（y）t=0:5:55； plot（t,y,'*'） legend（'原始數(shù)據(jù)分布趨勢圖'）運行結(jié)果：yua n(0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.584.02 4.64)*原貽數(shù)據(jù)分布趨勢圖2.51.5104050匕LI圖1原始數(shù)據(jù)分布趨勢：取基函數(shù)為：-0 =t1 =t2 ：由基函數(shù)和y求法方程組的系數(shù):- 12(0,0)八0(Xi)0(Xi)i

3、 412(1,1)八 i(Xi) i(Xi)12(0,1)八o(Xi)1(Xi)i 二12(1,2)八 l(Xi)2(Xi)12（0, ：2）八 0（X）2（Xi）7（,%）化,®0)(®2,®0)lA =（唧1）(®1,®1)（咒浮1）L（叩2）化,® 2）（鳴咒）一(f，%)1B =(f，)一（仁篤）一12 12(f,%)=2； w ®°(x)(仁鳴)=送 ®1(xJLi£i£：由這些系數(shù)，確定法方程組:A X 二B：解這個法方程組：AX二Bsum(k3.*k1);sum(k1.*

4、k2)sum(k2.*k2)i二i£12(2,2)八 2 (Xi)2 (Xi)i 412(f, 2)八 yi :2(x)idX = b，得到擬合函數(shù)：y = at + bt2+ct32程序：function a,b,c=xian(y0) t0=0:5:55;k1=tO;k2=t0*t0;k3=tO*tO.*tO;A=sum(k1.*k1)sum(k2.*k1) sum(k3.*k2);sum(k1.*k3) sum(k2.*k3) sum(k3.*k3);B=sum(k1*yO);sum(k2.*yO);sum(k3.*yO); x=p in v(A)*B;a=x(1,1);b=X(

5、2,1);c=x(3,1);t=0:55;y=a.*t+b.*t.A2+c.*t.A3; plot(t,y,'-') hold on plot(t0,y0,'*')legend('y=a*t+b*tA2+c*tA3擬合效果','真實值') 運行結(jié)果：a,b,c=xia n(0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64)a =0.2657b =-0.0053c =3.5168e-005(3)擬合的圖形，即上一題顯示的圖像45y=a*t+b*t2+c+t3 擬合效果

6、+真實值1.510203040506023圖2擬合函數(shù)y = at bt ct效果(4)，用于這種非線性模型的擬合：把其化作線性：y= tb f兩邊同時取以e為底的對數(shù)fln( y) = ln 二"b ln t：重復(fù)上面第二題的步驟進行，其中需要強調(diào)的是(0,0 )的點需要另外輸入，因為ln (0)不存在，為了在同圖出現(xiàn)，故對第二條擬合函數(shù)，取m=a n=b程序：function m,n,a,b,c=fei(y2,y0)%y2=y0除了 0 以外的數(shù)y1=log(y2);t1=5:5:55;n=le ngth(t1);k1= on es(1, n);k2=log(t1);A=sum(

7、k1.*k1) sum(k2.*k1);sum(k1.*k2) sum(k2.*k2);B=sum(k1.*y1);sum(k2.*y1);x=p in v(A)*B;m=exp(x(1,1);n =x(2,1);t=0:55;y=m*t.A n;plot(t,y,'-')hold ona,b,c=xia n(y0)plot(t,y,'-')hold onplot(t1,y2,'*',0,0,'*')legend('y=m*t.An擬合效果','y=a*t+b*tA2+c*tA3 擬合效果',

8、9;真實值')得到的擬合圖像：圖3兩種擬合函數(shù)擬合效果對比結(jié)論：在實際生活當中，不免需要對一組數(shù)據(jù)進行擬合，通過采用最佳的擬 合，找到一個近似的函數(shù)來研究數(shù)據(jù)的共性。 通過這一道題目，發(fā)現(xiàn)不同的函數(shù), 擬合效果差別也是蠻大的。第2題：特征值與特征向量用幕法求下列矩陣的主特征值與相應(yīng)的特征向量246一3-43|(1)A1 =3915(2)A?=-46341636331解：利用幕法求矩陣A的主特征值與相應(yīng)的特征向量，首先要給一個初始向量：1 ：定義一個和A行數(shù)一致的1列全一矩陣，即V。= 1Mjn A ' Oi=A“，為了方便計算，減少計算量，需要求出-1中按模最大的那個分量的值P

9、1,同時得到向量U1 =-1 ,由此可知Uo = VoR ：重復(fù)第二步，得到v2 P2 ：計算P2 -R，看其是否大于給定的誤差，如果大于，則令V1=V2, R = P2并從第三步開始重復(fù)；如果小于，則 P2為所求的按模最大的特征值，即主特征值，V2為其對應(yīng)的特征向量 流程圖:表1幕法求主特征值的流程圖程序：function v2,p2=maxtr(A,err)n=size(A,2);vO=on es( n,1);v1=A*v0;p1=max(v1);u1=v1/p1;v2=A*u1;p2=max(v2);i=1;while i>0if abs(p2-p1)v=errbreakendv1

10、=v2;P仁 p2;u1= v1/p1;v2=A*u1;p2=max(v2);i=i+1;endfprintf('求得矩陣A的按模最大的特征值 p2=:%10.8fn對應(yīng)的特征向量為v2=n',p2)disp(v2)246對問題(1), A, = 3915，誤差 err 1 =0.00014 16 36一運行結(jié)果：v2,p2=maxtr(2 4 6;3 9 15;4 16 36,0.0001);求得矩陣A的按模最大的特征值p2=:43.87998819對應(yīng)的特征向量為、8.155919.571943.8800v2=13_431對問題（2），A2 =-463，誤差 err 2 =

11、0.00011331 _運行結(jié)果：v2,p2=maxtr(3 -4 3;-4 6 3;3 3 1,0.0001); 求得矩陣A的按模最大的特征值p2=:8.86979621 對應(yīng)的特征向量為v2=-5.36028.86981.3380結(jié)論：對一個向量進行單位化，可以為下面的計算節(jié)省很多的計算步驟。 我 用到的程序，就是循環(huán)的賦予1，卩!的值，并通過1， P1求出V2 P2，最后達到 誤差要求，就求出了矩陣的主特征值與相應(yīng)的特征向量。第3題：微分方程求解用四階龍格-庫塔法求解初值問題* 2yxy (0汰蘭5) .y(0) = 2取步長h=0.25，并畫出0乞x乞5數(shù)值解的曲線。解：對于四階龍格-

12、庫塔法，我們在求第k+1個點的值時，必須知道第k個點的 值，在通過中間的一些變量，使結(jié)果更加精確。(1)：由題意可以確定初始值，即x0 = 0 y0 = 2(2):分別計算出：k1 =h f (xn, yn)k2 =h f (Xn + 陸 yn +%) k3 = h，f (Xn +%, yn +k%) N =h f (Xn +h,yn +k3)(3)：ynyn6(k12k22k3-k4)Xn 1=Xn h,由 n = 0開始，直至 X =5為方程的初值程序：function x,y=rk(x0,y0,a,b,h)%x0,y0 x(1)=x0;y(1)=y0;n=(b-a)/h;for j=1:

13、 nk1=h*f(x(j),y(j); k2=h*f(x(j)+h/2,y(j)+k1/2); k3=h*f(x(j)+h/2,y(j)+k2/2); k4=h*f(x(j)+h,y(j)+k3);y(j+1)=y(j)+1/6*(k1+2*k2+2*k3+k4); x(j+1)=x(j)+h;endplot(x,y,'r:')hold on xx=0:0.2:5;yy=2./(xx.A2+1);plot(xx,yy,'g-')legend('h=0.25的數(shù)值解','解析解')function z=f(x,y)z=-x*yA2;

14、運行結(jié)果： x,y=rk(0,2,0,5,0.25)x =Columns 1 through 60 0.2500 0.5000 0.7500 1.0000 1.2500Columns 7 through 121.5000 1.7500 2.0000Columns 13 through 183.0000 3.2500 3.5000Columns 19 through 214.5000 4.7500 5.0000y =Columns 1 through 62.0000 1.8823 1.5999Columns 7 through 120.6155 0.4924 0.4001Columns 13 t

15、hrough 180.2000 0.1730 0.1510Columns 19 through 210.0941 0.0849 0.07692.2500 2.5000 2.75003.7500 4.0000 4.25001.2799 1.0000 0.78060.3299 0.2759 0.23360.1328 0.1177 0.104921.8h=0.25的數(shù)值解 解析解1.61.41.210.80.60.40.2000.511.522.533.544.55圖4四階龍格-庫塔法擬合圖結(jié)論：四階龍格-庫塔的精度為四階，相對而言，擬合的效果非常好。但在計 算中，計算量比較大，借助計算機才可能實現(xiàn)

16、對數(shù)據(jù)的處理。第4題：求方程的根分別應(yīng)用二分法和牛頓法求解下面方程的根xsin x = 1要求滿足精度C -Xk <丄2解：(1)二分法： ：取上下限a、b,使f (a) f(b) ：0 ：取a、b的中點，即c二乎 ：看b-a >er是否成立，成立則繼續(xù)進行，不成立，則直接進行第四步計算f(c)，并判斷是否有f(a) f(c) ： 0 ,如果存在該關(guān)系式，則取b二c，不然 取a二c,重復(fù)第三步 ：輸出數(shù)值解b 流程圖:程序：function erfen(err)a=1;b=2;n=0;while abs(b-a)>=errn=n+1;c=(a+b)/2;if (a*si n(

17、a)-1)*(c*s in (c)-1)<0b=c;elsea=c;endfprintf('第%2d次迭代值為:%10.8fn',n,a)endfprintf('第 %2d次求得的近似值 x*=:%10.8fn',n,a)運行結(jié)果：erfen (0.00005)第1次迭代值為:1.00000000第2次迭代值為:1.00000000第3次迭代值為:1.00000000第4次迭代值為:1.06250000第5次迭代值為：1.09375000第6次迭代值為:1.10937500第7次迭代值為:1.10937500第8次迭代值為:1.11328125第9次迭代值

18、為:1.11328125第10次迭代值為:1.11328125第11次迭代值為:1.11376953第12次迭代值為:1.11401367第13次迭代值為:1.11413574第14次迭代值為:1.11413574第15次迭代值為:1.11413574第15次求得的近似值x*=:1.11413574(2),牛頓法:f (Xk)f'(Xk)k 二 0,1,2Xo：由k =0開始，不斷計算Xk1，直到XH - Xk <err，最后的數(shù)值解為xk出程序：function x1,n=mynewton(x0,err) x1=x0-f(x0)/df(x0);n=1;fprintf('第%2d次迭代值為：10.8fn',n,x1)while abs(x1-x0)>=errx0=x1;x1=x0-f(x0)/df(x0);n=n+1;fprintf('第%2d次迭代值為:%10.8fn',n,x1)endfprintf('第%2d次求得的近似值 x*=:%10.8fn',n,x1)function z=f(x)z=x*sin(x)-1;function z=df(x)z=sin(x)+x*cos(x);運行