1、數值計算方法上機實習題1 設，（1） 由遞推公式，從, 出發(fā)，計算；（2） ，, 用，計算；（3） 分析結果的可靠性及產生此現象的原因(重點分析原因)。解：（1）程序如下:clear allclcI=0.1822; %題中的已知數據for n=1:20; I=(-5)*I+1/n; %由遞推公式所得endfprintf('I20=%fn',I)M=0.1823; %與I的計算結果形成對比for i=1:20; M=(-5)*M+1/i; %由遞推公式所得endfprintf('M20=%fn',M)輸出結果為：I20=-11592559237.912731M20

2、=-2055816073.851284（2） 程序如下:clear allclcI=0; %賦予I20的初始值for n=0:19; I=(-1/5)*I+1/(5*(20-n); %有遞推公式得endfprintf('I0=%fn',I)M=10000; for i=0:19; M=(-1/5)*M+1/(5*(20-i);%有遞推公式得end fprintf('M0=%fn',M)輸出結果為：I0=0.182322M0=0.182322（3）由輸出結果可看出第一種算法為不穩(wěn)定算法，第二中算法為穩(wěn)定算法。由于誤差第一種算法為正向迭代算法，每計算一步誤差增長5倍

3、，雖然所給的初始值很接近，隨著n的增大，誤差也越來越大。第二種算法為倒向迭代算法，每計算一步誤差縮小5倍，雖然所給的初始值之間差很多，隨著n的增大，誤差也越來越小。2 求方程的近似根，要求，并比較計算量。（1） 在0，1上用二分法；（2） 取初值，并用迭代；（3） 加速迭代的結果；（4） 取初值，并用牛頓迭代法；（5） 分析絕對誤差。解：（1）程序如下（二分法）：clear allclca=0;b=1;f=(x)(exp(x)+10*x-2);%是定義函數句柄的運算符c=(a+b)/2;%取區(qū)間中點i=0;%分割次數while abs(f(c)>5*10(-4) %判斷f(x)的精度是否

4、滿足要求 if f(a)*f(c)<0 b=c;c=(a+b)/2; elseif f(b)*f(c)<0 a=c;c=(b+a)/2; end i=i+1;endfprintf('二分法運算次數為%in',i)fprintf('二分法計算結果為%fn',c)運行結果如下：二分法運算次數為13二分法計算結果為0.090515（2） 程序如下（不動點迭代）clear allclcx0=0;x=x0; for k=1:10000 %規(guī)定迭代次數上限 y=(2-exp(x)/10; %迭代結果存到y(tǒng)中 if abs(x-y)<5*10(-4) fpr

5、intf('初始值x0為%in迭代次數為%in',x0,k); break end x=y; if k=10000; fprintf('迭代次數超出上限%in',k); endendfprintf('迭代法計算結果為%fn',y);運行結果為：初始值x0為0迭代次數為4迭代法計算結果為0.090513（3） 程序如下（艾肯特迭代加速）clear allclcx0=0;x=x0; %給定迭代初值p(1000)=0; %先定義p矩陣的長度p(1)=x;for k=2:1000 ; %規(guī)定迭代次數上限 y=(2-exp(x)/10; z=(2-exp(

6、y)/10;x=x-(y-x)2)/(z-2*y+x); （4） 程序如下（牛頓法迭代）clear allclcx0=0;x=x0; %給定迭代初值for k=1:1000 ; %規(guī)定迭代次數上限y=x-(exp(x)+10*x-2)/(exp(x)+10); %牛頓計算公式 if abs(y-x)<5*10(-4) fprintf('初始值x0為%fn迭代次數為%in',x0,k);%艾特肯加速公式 p(k)=x; %p是用來存儲每一步迭代結果的矩陣 if abs(p(k-1)-p(k)<5*10(-4) fprintf('初始值x0為%fn迭代次數為%i

7、n',x0,k-1); break end if k=1000; fprintf('迭代次數超出上限%in',k);endendfprintf('艾特肯加速迭代法計算結果為%fn',x);運行結果為：初始值x0為0.000000迭代次數為2艾特肯加速迭代法計算結果為0.090525 break end x=y; if k=1000; fprintf('迭代次數超出上限%in',k); endendfprintf('牛頓迭代法計算結果為%fn',y);運行結果為：初始值x0為0.000000迭代次數為2牛頓迭代法計算結果為0

8、.090525(5) 分析絕對誤差通過指令求得方程精確解的近似解：>> solve('exp(x)+10*x-2=0') ans =1/5 - lambertw(0, exp(1/5)/10)>> 1/5 - lambertw(0, exp(1/5)/10)ans = 0.090525101307255各種計算方法的絕對誤差為:二分法的絕對誤差：1.0e-04 * 0.508986927450078不動點迭代方法的絕對誤差： 1.0e-04 * 0.121013072549997艾特肯加速迭代的絕對誤差： 1.0e-04 * 0.001013072550

9、016牛頓法的絕對誤差: 1.0e-04 * 0.001013072550016由上述各種算法計算出方程的值，二分法迭代了11次，收斂速度最慢，不動點迭代法迭代了4次，艾特肯迭代法迭代了2次，牛頓法迭代了2次，后兩種方法的收斂速度都比較快，但計算量大。由絕對誤差可以看出二分法的誤差最大，而艾特肯加速迭代和牛頓法誤差最小。3鋼水包使用次數多以后，鋼包的容積增大，數據如下：x23456789y6.428.29.589.59.7109.939.991011121314151610.4910.5910.6010.810.610.910.76試從中找出使用次數和容積之間的關系，計算均方差。（用擬合）解：

10、用已知函數模型來擬合，程序如下：先編輯函數function err=f31(f,x,y)a=f(1);b=f(2);c=f(3);err=0;for k=1:15 err=err+(c+x(k)*y(k)-(a*x(k)+b)2;end主函數如下：x=2 3 4 5 6 7 8 9 10 11 12 13 14 15 16;y=6.42 8.2 9.58 9.5 9.7 10 9.93 9.99 10.49 10.59 10.6 10.8 10.6 10.9 10.76;f,fval,exitflag=fminsearch(f31,0;0;0,optimset,x,y)%fminsearch函

11、數求解多元函數的最小值%f返回多元函數y=f(x)在初始值x0附近的局部極小值點%fval返回局部極小值%exitflag返回函數輸出條件值，exitflag=1說明函數收斂到一個解x；exitflag=0說明函數達到最大迭代次數；exitflag=-1說明輸出函數終止算法。fprintf('n 擬合出來的方程式為t (%7.4f+x)y=%7.4f+%7.4fxnn',f(3),f(2),f(1);z=linspace(x(1),x(end),15);Y=(f(1)*x+f(2)./(f(3)+x);plot(x,y,'r*-',z,Y,'b-'

12、;)legend('y','Y');title('非線性的最小二乘擬合');%均方差for i=1:15yt(i)=(f(1)*x(i)+f(2)./(f(3)+x(i);endc=0;for i=1:15 c=c+(y(i)-yt(i)2;endjfc=sqrt(c/15);fprintf('均方差為%fn',jfc)運行結果如下：擬合出來的方程式為 (-0.7110+x)y=-15.5024+11.2657x均方差為0.3316514設，分析下列迭代法的收斂性，并求的近似解及相應的迭代次數。（1） JACOBI迭代；（2） G

13、AUSS-SEIDEL迭代；（3） SOR迭代（取，找到迭代步數最少的）。解：(1)JACOBI迭代(程序如下)：先編輯JACOBI函數：function tx=jacobi(A,b,imax,x0,tol) %初始值x0,次數imax,精度toldel=10-10;tx=x0;n=length(x0);for i=1:n dg=A(i,i); if abs(dg)<del disp('對角元太小');%防止現溢出現象 return endendfor k=1:imax %jacobi循環(huán) for i=1:n sm=b(i); for j=1:n if j=i sm=sm

14、-A(i,j)*x0(j); end end x(i)=sm/A(i,i);%x(1)到x(n)即完成一次迭代 end tx=tx;x;%矩陣中又加一行。 if norm(x-x0)<tol return else x0=x; end end主函數程序如下：clear allclcA=4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1 0 -1 4 -1 0;0 -1 0 -1 4 -1;0 0 -1 0 -1 4;b=0 5 -2 5 -2 6;x0=0 0 0 0 0 0;imax=100;%迭代次數tol=10-4;%精度tx=jacobi

15、(A,b,imax,x0,tol);for j=1:length(tx)fprintf('%d %f %f %f %f %f %fn',j,tx(j,1),tx(j,2),tx(j,3),tx(j,4),tx(j,5),tx(j,6);end運行結果如下:1 0.000000 0.000000 0.000000 0.000000 0.000000 0.0000002 0.000000 1.250000 -0.500000 1.250000 -0.500000 1.50000028 0.999947 1.999964 0.999927 1.999964 0.999927 1.99

16、997329 0.999982 1.999950 0.999975 1.999950 0.999975 1.999964(2)GAUSS-SEIDEL迭代(程序如下)：先編輯GAUSS-SEIDEL函數：function tx=gseidel(A,b,imax,x0,tol) %初始值x0,次數imax,精度toldel=10-10;tx=x0;n=length(x0);%tx是個二維矩陣，存儲著每一步迭代的結果for i=1:n dg=A(i,i); if abs(dg)<del disp('對角元太小'); return endend for k=1:imax %G-

17、S循環(huán) x=x0; for i=1:n sm=b(i); for j=1:n if j=i sm=sm-A(i,j)*x(j); end end x(i)=sm/A(i,i); end tx=tx;x; if norm(x-x0)<tol return else x0=x; end end主程序如下:clear allclcA=4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1 0 -1 4 -1 0;0 -1 0 -1 4 -1;0 0 -1 0 -1 4;b=0 5 -2 5 -2 6;x0=0 0 0 0 0 0;imax=100;tol=

18、10-4;tx=gseidel(A,b,imax,x0,tol);for j=1:length(tx)fprintf('%d %f %f %f %f %f %fn',j,tx(j,1),tx(j,2),tx(j,3),tx(j,4),tx(j,5),tx(j,6)end運行結果如下：1 0.000000 0.000000 0.000000 0.000000 0.000000 0.0000002 0.000000 1.250000 -0.187500 1.203125 0.113281 1.48144515 0.999917 1.999910 0.999924 1.999931

19、0.999943 1.99996716 0.999960 1.999957 0.999964 1.999967 0.999973 1.999984(3)SOR迭代(程序如下)：先編輯SOR函數：function tx=sor(A,b,imax,x0,tol,w) %初始值x0,次數imax,精度toldel=10-10;tx=x0;n=length(x0);for i=1:n dg=A(i,i); if abs(dg)<del disp('對角元太小'); return endend for k=1:imax %SOR循環(huán) x=x0; for i=1:n sm=b(i);

20、 for j=1:n if j=i sm=sm-A(i,j)*x(j); %先高斯迭代 end end x(i)=sm/A(i,i); x(i)=w*x(i)+(1-w)*x0(i);%比較高斯與超松弛迭代公式，補上省缺的項 end tx=tx;x; if norm(x-x0)<tol return else x0=x; endend 主程序如下：clear allclcA=4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1 0 -1 4 -1 0;0 -1 0 -1 4 -1;0 0 -1 0 -1 4;b=0 5 -2 5 -2 6;x0=0

21、 0 0 0 0 0;imax=100;tol=10-4;w=1.2;%給定不同的w，w=0.1:0.1:1.9，找出使SOR迭代法收斂速度最快tx=sor(A,b,imax,x0,tol,w);for j=1:size(tx,1)fprintf('%d %f %f %f %f %f %fn',j,tx(j,1),tx(j,2),tx(j,3),tx(j,4),tx(j,5),tx(j,6)end 運行結果如下:1 0.000000 0.000000 0.000000 0.000000 0.000000 0.0000002 0.000000 1.500000 -0.150000

22、 1.455000 0.286500 1.84095010 0.999957 2.000022 0.999979 1.999982 1.000018 1.99999211 1.000010 1.999998 0.999996 2.000011 0.999997 1.999999w值得范圍在02之間，SOR迭代才會收斂，令w=0.1:0.1:1.9找出w*，使得SOR迭代的收斂速度最快。w=0.1時需要迭代101次；w=1.9時需要迭代101次；w=1時需要迭代16次；w=0.9時需要迭代20次；w=1.1時需要迭代13次；w=1.3時需要迭代13次；w=1.2時需要迭代11次。有這幾次試代可得

23、出w=1.2時迭代次數最少，SOR迭代法收斂速度最快。總結:由于A為嚴格對角占優(yōu)矩陣，根據定理可知雅可比迭代和GS迭代都收斂，SOR迭代收斂的必要條件是0<w<2，即SOR迭代也收斂。三種迭代算法的結果分析:(1) JACOBI迭代是根據A分解成上三角、下三角和對角陣，將線性方程組的求解歸結為對角方程組求解過程的重復，思路簡單易于編程。迭代次數比較多，收斂速度慢。(2) GAUSS-SEIDEL迭代，是在JACOBI收斂的基礎上將計算出來的未知量馬上投入下一個迭代方程中使用，使得迭代的數度加快，效果更好。但GS迭代與JACOBI迭代的收斂域并不能相互包含，所以不能相互代替。但如果兩

24、者都收斂時，一般說GS迭代比JACOBI迭代的收斂速度快。(3)SOR迭代為了加快收斂速度，在GS的基礎上加入一修正參數即松弛因子w，使得收斂速度更快。但選著不同的w，收斂速度不一樣，為了達到最好的效果，要選著最合適的松弛因子。5用逆冪迭代法求最接近于11的特征值和特征向量，準確到。解：程序如下先編輯函數文件：function l,v=fanpower(A,v0,p,ep,max1)%p為A最接近的特征值，ep為精度，max1為最大迭代次數n=length(A);B=A-p*eye(n);k=0;err=1;v0=v0'%求v0的轉秩while (k<=max1)&(er

25、r>ep) m j=max(abs(v0);%找出v0中模最大的元素，j為此元素的序號 y=v0/(v0(j); x0=0 1 0;%解B*v1=y'求逆矩陣占用空間較大,用高斯迭代。 tx=gseidel(B,y',100,x0,10-7); r q=size(tx);%得出的tx是一個二維矩陣 v1=tx(r,1) tx(r,2) tx(r,3)' %提取出最后一行的三個元素 err=abs(1/norm(v1,2)-1/norm(v0,2); %norm(v,2)求矩陣v的2范數，即：最大特征值 v0=v1; k=k+1;end m j=max(abs(v0

26、);l=p+(1/v0(j);v=y;fprintf('最接近給定數的特征值是ln');fprintf('對應的特征向量為v'); 主程序如下：clear allclcA=6 3 1;3 2 1;1 1 1;v0=1 1 1;p=11;ep=0.001;max1=1000;l,v=fanpower(A,v0,p,ep,max1)運行結果如下:最接近給定數的特征值是l對應的特征向量為vl = 7.8731v = 1.0000 0.5493 0.2256總結:此方法為結合原點平移的反冪法，可以根據矩陣A的已知的特征值的粗略近似值p，計算出與數p最接近的特征值及特征向

27、量。6用經典R-K方法求解初值問題（1）， ；（2）， 。4解：（1）程序如下：clc;clear;f=(x,y1,y2)(-2*y1+y2+2*sin(x); g=(x,y1,y2) (y1-2*y2+2*cos(x)-2*sin(x);h=0.1;y1(1)=2;y2(1)=3;x(1)=0;for i=1:100;K1=f(x(i),y1(i),y2(i);L1=g(x(i),y1(i),y2(i);K2=f(x(i)+0.5*h,y1(i)+0.5*h*K1,y2(i)+0.5*h*L1);L2=g(x(i)+0.5*h,y1(i)+0.5*h*K1,y2(i)+0.5*h*L1);K

28、3=f(x(i)+0.5*h,y1(i)+0.5*h*K2,y2(i)+0.5*h*L2);L3=g(x(i)+0.5*h,y1(i)+0.5*h*K2,y2(i)+0.5*h*L2);K4=f(x(i)+h,y1(i)+h*K3,y2(i)+h*L3);L4=g(x(i)+h,y1(i)+h*K3,y2(i)+h*L3);x(i+1)=x(i)+h;y1(i+1)=y1(i)+h*(1/6)*(K1+2*K2+2*K3+K4);y2(i+1)=y2(i)+h*(1/6)*(L1+2*L2+2*L3+L4);endx=0:0.1:10;for i=1:101;Y1(i)=2*exp(-x(i)

29、+sin(x(i);Y2(i)=2*exp(-x(i)+cos(x(i);endc=y1-Y1;d=y2-Y2;subplot(2,1,1)plot(x,y1,'r-',x,y2,'b-','LineWidth',4)legend('y1','y2');title('R-K四階龍格庫塔算法下方程組的解');ylabel('y1曲線 y2曲線')subplot(2,1,2)plot(x,c,'r-',x,d,'b-')legend('c'

30、,'d');title('誤差值');ylabel('c曲線 d曲線')x;y1 ;c;y2;d'（2） 程序如下：clc;clear;f=(x,y1,y2)(-2*y1+y2+2*sin(x); g=(x,y1,y2)(998*y1-999*y2+999*cos(x)-999*sin(x);h=0.001;y1(1)=2;y2(1)=3;x(1)=0;for i=1:10000;K1=f(x(i),y1(i),y2(i);L1=g(x(i),y1(i),y2(i);K2=f(x(i)+0.5*h,y1(i)+0.5*h*K1,y2(i)

31、+0.5*h*L1);L2=g(x(i)+0.5*h,y1(i)+0.5*h*K1,y2(i)+0.5*h*L1);K3=f(x(i)+0.5*h,y1(i)+0.5*h*K2,y2(i)+0.5*h*L2);L3=g(x(i)+0.5*h,y1(i)+0.5*h*K2,y2(i)+0.5*h*L2);K4=f(x(i)+h,y1(i)+h*K3,y2(i)+h*L3);L4=g(x(i)+h,y1(i)+h*K3,y2(i)+h*L3);x(i+1)=x(i)+h;y1(i+1)=y1(i)+h*(1/6)*(K1+2*K2+2*K3+K4);y2(i+1)=y2(i)+h*(1/6)*(L1+2*L2+2*L3+L4);endx=0:0.001:10;Y1=2*exp(-x)+sin(x);Y2=2*exp(-x)+cos(x);%精確解x;y1 ;Y1;y2;Y2'subplot(121)plot(x,y1,