數(shù)值計(jì)算方法上機(jī)實(shí)習(xí)題

1、數(shù)值計(jì)算方法上機(jī)實(shí)習(xí)題1 設(shè)，（1） 由遞推公式，從I0=0.1824, 出發(fā)，計(jì)算；（2） ，, 用，計(jì)算；(1)由I0計(jì)算I20遞推子程序：function f=fib(n,i)if n>=1 f=fib(n-1,i)*(-5)+(1/(n);elseif n=0 f=i;end計(jì)算和顯示程序：I=0.1824;I1=0.1823;fib1=fib(20,I);fib2=fib(20,I1);fprintf('I_0=0.1824時，I_20=%dn',fib1);fprintf('I_0=0.1823時，I_20=%dn',fib2);計(jì)算結(jié)果顯示：

2、I_0=0.1824時，I_20=7.480927e+09I_0=0.1823時，I_20=-2.055816e+09(2)由I20計(jì)算I0程序：n=21;i1=0;i2=10000;f1=i1;f2=i2;while n=0; f1=f1*(-1/5)+(1/(5*n); f2=f2*(-1/5)+(1/(5*n); n=n-1;endfprintf('I_20=0 時，I_0=%4.8fn',f1);fprintf('I_20=10000時，I_0=%4.8fn',f2);計(jì)算結(jié)果顯示：I_20=0 時，I_0=0.18232156I_20=10000時，I

3、_0=0.18232156（3） 分析結(jié)果的可靠性及產(chǎn)生此現(xiàn)象的原因(重點(diǎn)分析原因)。答：第一個算法可得出e0=I0-I0*en=In-In*=5ne0易知第一個算法每一步計(jì)算都把誤差放大了5倍，n次計(jì)算后更是放大了5n倍，可靠性低。 第二個算法可得出en=In-In*e0=15nen可以看出第二個算法每一步計(jì)算就把誤差縮小5倍，n次后縮小了5n倍，可靠性高。2 求方程的近似根，要求，并比較計(jì)算量。（1） 在0，1上用二分法；(1) 0，1上的二分法二分法子程序：function root,n=EFF3(f,x1,x2)%第二題(1)二分法f1=subs(f,symvar(f),x1);%函數(shù)

4、在x=x1的值f2=subs(f,symvar(f),x2);%x=x2n=0;%步數(shù)er=5*10-4;%誤差if(f1=0) root=x1; return;elseif(f2=0) root=x2; return;elseif(f1*f2>0) disp('兩端點(diǎn)函數(shù)值乘積大于0！'); return;else while(abs(x1-x2)>er)%循環(huán) x3=(x1+x2)/2; f3=subs(f,symvar(f),x3); n=n+1; if(f3=0) root=x3; break; elseif(f1*f3>0) x1=x3; else

5、x2=x3; end end root=(x1+x2)/2;%while循環(huán)少一步需加上end計(jì)算根與步數(shù)程序：fplot(x) exp(x)+10*x-2,0,1);grid on;syms x;f=exp(x)+10*x-2;root,n=EFF3(f,0,1);fprintf('root=%6.8f ,n=%d n',root,n);計(jì)算結(jié)果顯示：root=0.09057617 ,n=11 （2） 取初值，并用迭代；(2) 初值x0=0迭代迭代法子程序：function root,n=DDF(g,x0,err,max) （接下頁）%root根，n+1步數(shù)，g函數(shù)，x0初值

6、，err誤差，max最大迭代次數(shù)X(1)=x0;for n=2:max X(n)=subs(g,symvar(g),X(n-1); c=abs(X(n)-X(n-1); root=X(n); if(c<err) break; end if n=max disp('超出預(yù)設(shè)迭代次數(shù)'); endendn=n-1;%步數(shù)減1計(jì)算根與步數(shù)程序：syms x;f=(2-exp(x)/10; （接下頁）x0=0; err=5*10(-4);max=100;root,n=DDF(f,x0,err,max);fprintf('root=%6.8f ,n=%d n',ro

7、ot,n);計(jì)算結(jié)果顯示：root=0.09051262 ,n=4（3） 加速迭代的結(jié)果;(3) 加速迭代加速迭代計(jì)算程序：x0=0;err=5*10(-4);max=100;syms x;g=x-(f(x)-x)2/(f(f(x)-2*f(x)-x);root,n=DDF(g,x0,err,max);fprintf('root=%6.8f, n=%d',root,n);程序函數(shù)設(shè)置：function g=f(x)g=(2-exp(x)/10;計(jì)算結(jié)果顯示：root=0.09048375, n=2（4） 取初值，并用牛頓迭代法；(4) 牛頓迭代法牛頓迭代法子程序：functio

8、n root,n=NDDDFx(g,x0,err,max)%root根，n+1步數(shù)，g函數(shù)，x0初值，err誤差，max最大迭代次數(shù)X(1)=x0;for n=2:max X(n)=subs(g,symvar(g),X(n-1); c=abs(X(n)-X(n-1); root=X(n); if(c<err)|(root=0) break; end if n=max disp('超出預(yù)設(shè)迭代次數(shù)'); （接下頁） endendn=n-1;%步數(shù)減1牛頓迭代法計(jì)算程序：x0=0;err=5*10(-4);max=100;syms x;g=exp(x)+10*x-2;g=x-

9、(g/diff(g);root,n=NDDDFx(g,x0,err,max);fprintf('root=%6.8f, n=%d n',root,n)；計(jì)算結(jié)果顯示：root=0.09052511, n=2（5） 分析絕對誤差。答：可以看到，在同一精度下，二分法運(yùn)算了11次，題設(shè)迭代算式下運(yùn)算了4次，加速迭代下運(yùn)算了2次，牛頓迭代下運(yùn)算了2次。因不動點(diǎn)迭代法和二分法都是線性收斂的，但二分法壓縮系數(shù)比題設(shè)迭代方法大，收斂速度較慢。加速迭代速度是超線性收斂，牛頓法是二階，收斂速度快。3鋼水包使用次數(shù)多以后，鋼包的容積增大，數(shù)據(jù)如下：x23456789y6.428.29.589.59

10、.7109.939.991011121314151610.4910.5910.6010.810.610.910.76試從中找出使用次數(shù)和容積之間的關(guān)系，計(jì)算均方差。（用擬合）擬合曲線程序：x=2:16;y=6.42 8.2 9.58 9.5 9.7 10 9.93 9.99 10.49 10.59 10.60 10.8 10.6 10.9 10.76;f,fval=fminsearch(delta1,0,0,0,optimset,x,y);%fminsearch為求解多元函數(shù)的最小值函數(shù)%f為多元函數(shù)初值x0附近的極小值點(diǎn)%fval為極小值k=2:100;y1=(f(1).*k+f(2)./(

11、f(3)+k);figure1=figure('Color',1 1 1);gxt=plot(x,y,'xr',k,y1,'k-');legend('原數(shù)據(jù)','擬合曲線','Location','northwest'); %y為數(shù)據(jù)點(diǎn)連接曲線，y1為擬合曲線title('函數(shù)y=(ax+b)/(x+c)的擬合','FontSize',14,'FontWeight','Bold');xlabel('次數(shù)'

12、,'FontSize',14,'FontWeight','Bold');ylabel('容積','FontSize',14,'FontWeight','Bold');set(gxt,'LineWidth',1.5);grid on;%計(jì)算均方差for i=1:15y2(i)=(f(1).*x(i)+f(2)./(f(3)+x(i); （接下頁）endj=0;for i=1:15j=j+(y(i)-y2(i)2;endjfc=sqrt(j/15);fprintf(

13、9;擬合出的方程為:(x+%4.4f)y=%4.4fx+%4.4f n均方差為:%4.8f n',f(3),f(1),f(2),jfc);構(gòu)造函數(shù)子程序：function delta=delta1(f,x,y)a=f(1);b=f(2);c=f(3);delta=0;for k=1:15 delta=delta+(x(k)+c)*y(k)-(a*x(k)+b)2;end計(jì)算結(jié)果顯示：擬合出的方程為:(x+-0.7110)y=11.2657x+-15.5024 均方差為:0.33165089總結(jié)：指標(biāo)選擇，因題設(shè)方程為非線性的，要轉(zhuǎn)化為線性方程故需提指標(biāo)為：其駐點(diǎn)方程為:計(jì)算結(jié)果顯示：4

14、設(shè)，分析下列迭代法的收斂性，并求的近似解及相應(yīng)的迭代次數(shù)。（1） JACOBI迭代；(1) Jacobi迭代Jacobi迭代子程序：function x,k=Jacobifl2(A,b,x0,eps,max1)%jacobi按矩陣的分量算法%x0初值，eps誤差，max1最大迭代次數(shù)n=length(x0);for k=1:max1 x(1)=(b(1)-A(1,2:n)*x0(2:n)/A(1,1); for i=2:n x(i)=(b(i)-A(i,1:i-1)*x0(1:i-1)-A(i,i+1:n)*x0(i+1:n)/A(i,i); end if norm(x'-x0,2)&

15、lt;eps return else x0=x' endendJacobi迭代計(jì)算程序：A=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'eps=10(-4);max1=200;X,k=Jacobifl2(A,b,x0,eps,max1);fprintf('近似解x=n%4.8fn %4.8f n %4.8f n %4.8fn %4.8f n %4.8fn迭代次數(shù)n=%

16、dn',X(1),X(2),X(3),X(4),X(5),X(6),k);計(jì)算結(jié)果顯示：近似解x=0.99998177 1.99995018 0.99997509 1.99995018 0.99997509 1.99996353迭代次數(shù)n=28（2） GAUSS-SEIDEL迭代；(2) GAUSS-SEIDEL迭代GAUSS-SEIDEL迭代子程序：function x,k=GSDD(A,b,x0,eps,max1)n=length(x0);for k=1:max1 x(1)=(b(1)-A(1,2:n)*x0(2:n)')/A(1,1); for i=2:n x(i)=(b

17、(i)-A(i,1:i-1)*x(1:i-1)'-A(i,i+1:n)*x0(i+1:n)')/A(i,i); end if norm(x-x0,2)<eps return else x0=x; endendGAUSS-SEIDEL迭代計(jì)算程序：A=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;eps=10(-4);max1=200;X,k=GSDD(A,b,x0,eps

18、,max1);fprintf('近似解x=n%4.8fn %4.8f n %4.8f n %4.8fn %4.8f n %4.8fn迭代次數(shù)n=%dn',X(1),X(2),X(3),X(4),X(5),X(6),k);計(jì)算結(jié)果顯示：近似解x=0.99996014 1.99995676 0.99996351 1.99996662 0.99997253 1.99998401迭代次數(shù)n=15（3） SOR迭代（取，找到迭代步數(shù)最少的）。(3) SOR迭代SOR迭代子程序：function x,k,m=SOR1(A,b,x0,eps,max1)n=length(x0);a=x0;fo

19、r s=1:19 x0=a; w(s)=s/10; x(1)=x0(1)+w(s)*(b(1)-A(1,1:n)*x0(1:n)')/A(1,1); for i=2:n x(i)=x0(i)+w(s)*(b(i)-A(i,1:i-1)*x(1:i-1)'-A(i,i:n)*x0(i:n)')/A(i,i); end k(s)=0; while norm(x-x0,2)>=eps x0=x; x(1)=x0(1)+w(s)*(b(1)-A(1,1:n)*x0(1:n)')/A(1,1); for i=2:n x(i)=x0(i)+w(s)*(b(i)-A(i

20、,1:i-1)*x(1:i-1)'-A(i,i:n)*x0(i:n)')/A(i,i); end k(s)=k(s)+1; if (k(s)>=max1) disp('超出迭代次數(shù)上限！n'); return; end end Datas=x;endmin1,index=min(k);x=Dataindex;k=min1;m=w(index);fprintf('n解x=n%4.8fn%4.8fn%4.8fn%4.8fn%4.8fn%4.8f',x);fprintf('n步數(shù)最小omega=%2.2f, 步數(shù)n=%d 。n',

21、m,k);SOR迭代計(jì)算程序：A=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;eps=10(-4);max1=1000;X,k,m=SOR1(A,b,x0,eps,max1);fprintf('n近似解x=n%4.8fn%4.8fn%4.8fn%4.8fn%4.8fn%4.8f',x);fprintf('n步數(shù)最小omega=%2.1f, 步數(shù)n=%dn',m

22、,k);計(jì)算結(jié)果顯示：近似解x=1.000009771.999997820.999995672.000010710.999996621.99999923步數(shù)最小omega=1.20, 步數(shù)n=9 ?？偨Y(jié)：Jacobi迭代次數(shù)為28次；GAUSS-SEIDEL迭代次數(shù)為15次；SOR迭代次數(shù)在w=1.2時達(dá)到最小次數(shù)9次。系數(shù)矩陣嚴(yán)格對角占優(yōu)，則Jacobi迭代法與Gauss-Seidel迭代法均收斂；又因系數(shù)矩陣對稱正定，且0w2，則SOR迭代法也收斂。三種迭代法的結(jié)果分析及比較：（1）Jacobi法，其設(shè)計(jì)思想是將線性方程組的求解歸結(jié)為對角方程組求解過程的重復(fù)，計(jì)算規(guī)則簡單而易于編寫計(jì)算程序

23、。但收斂速度慢，譜半徑相對偏大（2）Gauss-Seidel充分利用新信息進(jìn)行計(jì)算，其迭代的效果比Jacobi迭代好，譜半徑比Jacobi法小。（3）SOR迭代法，計(jì)算公式簡單、編制程序容易，是求解大型稀疏方程組的一種有效方法，且若松弛因子選擇得當(dāng)超松弛法收斂速度比Gauss-Seidel迭代和Jacobi迭代都要快。 5用逆冪迭代法求最接近于11的特征值和特征向量，準(zhǔn)確到。逆冪迭代法子程序：function lamda,V,k=NMSF2(A,v,eps,p)n,n=size(A);A=A-p*eye(n);v0=v;tmax,tindex=max(abs(v0);lamd0=v0(tind

24、ex);u0=v0/lamd0;flag=0;k=0;while(flag=0) V=Au0; tmax,tindex=max(abs(V); lamd1=V(tindex); u0=V/lamd1; if (abs(lamd0)(-1)-(lamd1)(-1)<=eps flag=1; end lamd0=lamd1; k=k+1;end （接下頁）lamda=(lamd1)(-1)+p;V=u0'逆冪迭代計(jì)算程序：A=6 3 1;3 2 1;1 1 1;v=1 1 1'eps=0.001;p=11;lamda,U,k=NMSF2(A,v,eps,p);fprintf(

25、'最接近11的特征值為:%4.8fn特征向量：n%4.8fn %4.8fn%4.8fn',lamda,U(1),U(2),U(3);計(jì)算結(jié)果顯示：最接近于11的特征值為：7.87313712特征向量： 1.00000000 0.54922698 0.225456186用經(jīng)典R-K方法求解初值問題（1）， ；（2）， 。和精確解比較，進(jìn)行誤差分析得到結(jié)論，圖形顯示精確解和數(shù)值解。經(jīng)典四階R-K方法R-K程序：f=(x,y) -2*y(1)+y(2)+2*sin(x);y(1)-2*y(2)+2*cos(x)-2*sin(x);g=(x,y) -2*y(1)+y(2)+2*sin(

26、x);998*y(1)-999*y(2)+999*cos(x)-999*sin(x); （接下頁）x=0:0.1:10;Y1=2*exp(-x)+sin(x);Y2=2*exp(-x)+cos(x);x1=0:10/10000:10;Y3=2*exp(-x1)+sin(x1);Y4=2*exp(-x1)+cos(x1);N=100;N1=10000;a1,b1=ode23(f,0:0.1:10,2 3);a2,b2=ode23(g,0:1/1000:10,2 3);a3,b3=ode45(f,0:0.1:10,2 3);a4,b4=ode45(g,0:1/1000:10,2 3);n=leng

27、th(b1);figure1=figure('Color',1 1 1);dbt=plot(x,b1(:,1),'r-',x,b1(:,2),'k-',x,Y1,'g-',x,Y2,'y-');legend('y1','y2','精確解Y1','精確解Y2'); title('四階RK算法下的方程組的數(shù)值解&方程的精確解','FontSize',14,'FontWeight','Bold&

28、#39;);ylabel('Y','FontSize',14,'FontWeight','Bold');xlabel('X','FontSize',14,'FontWeight','Bold');set(dbt,'LineWidth',1.5);grid on;Q=b1(:,1)-Y1'W=b1(:,2)-Y2'Q1=b2(:,1)-Y3'W1=b2(:,2)-Y4'Q2=b3(:,1)-Y1'W2=b3(:,2

29、)-Y2'Q3=b4(:,1)-Y3'W3=b4(:,2)-Y4' ER1=max(abs(Q);ER2=max(abs(W);ER3=max(abs(Q1);ER4=max(abs(W1);ER5=max(abs(Q2);ER6=max(abs(W2);ER7=max(abs(Q3);ER8=max(abs(W3);fprintf('n（1）中y1在ode23下求出的最大誤差為：%4.8fn（1）中y2在ode23下求出的最大誤差為：%4.8fn',ER1,ER2);fprintf('n（2）中y1在ode23下求出的最大誤差為：%4.8fn（

30、2）中y2在ode23下求出的最大誤差為：%4.8fn',ER3,ER4);fprintf('n（1）中y1在ode45下求出的最大誤差為：%4.8fn（1）中y2在ode45下求出的最大誤差為：%4.8fn',ER5,ER6);fprintf('n（2）中y1在ode45下求出的最大誤差為：%4.8fn（2）中y2在ode45下求出的最大誤差為：%4.8fn',ER7,ER8);計(jì)算與顯示程序：（1）中y1在ode23下求出的最大誤差為：0.00169850（1）中y2在ode23下求出的最大誤差為：0.00207422（2）中y1在ode23下求出的最大誤差為：0.00000583（2）中y2在ode23下求出的最大誤差為：0.00581329（1）中y1在ode45下求出的最大誤差為：0.00052155（1）中y2在ode45下求出的最大誤差為：0.00045793（2）中y1在ode45下求出的最大誤差為：0.00000276（2）中y2在ode45下求出的最大誤差為：0.00275331總結(jié)：對比2種方程在ode23，ode45命令下，與精確值的誤差可以發(fā)現(xiàn)，在此題中，ode45的誤差相對于ode23的誤