北科大丁军计算方法作业

姓名：学号：班级：学院：2018年11月25日3-1试验目的：考察不动点迭代法的局部收敛性试验内容：分别构造方程和，至少采用3种迭代法，迭代100次，考察收敛性，改变初值符号，再做迭代。分析收敛与发散的原因。(1)迭代原理：若实数满足，称为函数的一个不动点，迭代称为不动点迭代，称为迭代函数。由不动点方程建立迭代法，其中称为初值，需要预先给定。(2)方程分别对应下列不同形式的不动点方程：1.2.3.取，按迭代，并分析收敛性。不动点迭代法代码1.function[p,k]=fone(p0,max,tol)k=1;whilek<=maxp=3*p0+3-exp(p0);ifabs(p-p0)<tolbreak;endk=k+1;p0=p;enddisp(p);disp(k)运行结果：2.function[p,k]=ftwo(p0,max,tol)k=1;whilek<=maxp=(exp(p0)-3)/2;ifabs(p-p0)<tolbreak;endk=k+1;p0=p;enddisp(p);disp(k)运行结果：3.function[p,k]=fthree(p0,max,tol)k=1;whilek<=maxp=log(2*p0+3);ifabs(p-p0)<tolbreak;endk=k+1;p0=p;enddisp(p);disp(k)运行结果：(3)方程分别对应下列不同形式的不动点方程：1.2.3.取，按迭代，并分析收敛性。不动点迭代法代码1.function[p,k]=fone(p0,max,tol)k=1;whilek<=maxp=3*p0^5+5*p0^2+p0-10;ifabs(p-p0)<tolbreak;endk=k+1;p0=p;enddisp(p);disp(k)运行结果： 2.function[p,k]=ftwo(p0,max,tol)k=1;whilek<=maxp=(10/(3*p0^3+5))^(1/2);ifabs(p-p0)<tolbreak;endk=k+1;p0=p;enddisp(p);disp(k)运行结果：3.function[p,k]=fthree(p0,max,tol)k=1;whilek<=maxp=p0-(3*p0^5+5*p0^2-10)/(15*p0^4+10*p0);ifabs(p-p0)<tolbreak;endk=k+1;p0=p;enddisp(p);disp(k)运行结果：(4)结果分析：一般地，一个方程的不动点迭代有多种形式，有的收敛，有的发散。在收敛的方法中，收敛的快慢也不相同。第二种和第三种方法都收敛，第一种方法发散，说明迭代法具有局部性。3-5试验目的：应用不同迭代法求解代数方程试验内容：分别采用二分法、Newton法、Newton下山法、割线法求解方程在[0.1,1]中的根；精确到10-4。1.二分法算法：为连续函数，且由题意可知[0.1，1]为含根区间，令a=0.1，b=1，取p=(a+b)/2。若f(p)=0则p是方程f(x)=0的解；若f(a)f(p)<0则根在(a,p)内，取a1=a,b1=p；否则根在区间(p,b)内，取a1=p,b1=b。重复上述过程直到达到精度要求为止。程序：functiony=f_3_5(x)y=600*x.^4-550*x.^3+200*x.^2-20*x-1;endfunction[p,k]=f_3_5_Bisection(max,a,b,e)formatlongfork=1:maxp=(b+a)/2;iff_3_5(p)==0||(b-a)/2<ebreak;elseiff_3_5(a)*f_3_5(p)>0a=p;elseb=p;endend运行结果：2.Newton法算法：建立牛顿迭代格式直到小于精度要求时迭代结束，将作为结果输出。程序：functiony=f_3_5(x)y=600*x.^4-550*x.^3+200*x.^2-20*x-1;endfunctiony=f2_3_5(x)y=600*4*x.^3-550*3*x.^2+400*x-20;endfunction[p,k]=f3_5_Newton(max,p0,e)formatlongfork=1:maxp=p0-f_3_5(p0)/f2_3_5(p0);ifabs(p-p0)<ebreak;endp0=p;end运行结果：3.Newton下山法算法：建立迭代格式 其中是下山因子，在每步迭代中依次取1，1/2，1/4，…直到满足精度要求。程序：functiony=f_3_5(x)y=600*x.^4-550*x.^3+200*x.^2-20*x-1;endfunctiony=f2_3_5(x)y=600*4*x.^3-550*3*x.^2+400*x-20;endfunction[p,k]=f3_5_Newton2(max,p0,tol,e)formatlongfork=1:maxy=1;p=p0-y*f_3_5(p0)/f2_3_5(p0);while(abs(f_3_5(p))>=abs(f_3_5(p0))&y>e)y=y/2;p=p0-y*f_3_5(p0)/f2_3_5(p0);endifabs(p-p0)<tolbreak;endp0=p;end运行结果：4.割线法算法：建立迭代格式 直到小于精度要求时迭代结束，将作为结果输出。程序：function[p,k]=f3_5_Secant(max,a,b,e)formatlongfork=1:maxp=b-f_3_5(b)*(b-a)/(f_3_5(b)-f_3_5(a));ifabs(p-b)<ebreak;enda=b;b=p;end运行结果：4-2试验目的：熟悉LU分解法求解方程。试验内容：用MATLAB编程做Doolittle分解，求解方程组。算法：对进行Doolittle分解使得，则线性方程组的求解问题就转化为对三角形方程组的求解。Doolittle分解过程为：然后依次求解和：程序：function[L,U,x]=f4_2(a,b)n=length(a);L=eye(n);U(1,:)=a(1,:);fork=2:nifU(k-1,k-1)==0disp('分解失败');return;endL(k:n,k-1)=a(k:n,k-1)/U(k-1,k-1);U(k,k:n)=a(k,k:n)-L(k,1:k-1)*U(1:k-1,k:n);ifk<na(k+1:n,k)=a(k+1:n,k)-L(k+1:n,1:k-1)*U(1:k-1,k);endendy=zeros(n,1);y(1)=b(1);fori=2:ny(i)=b(i)-L(i,1:i-1)*y(1:i-1);endx=zeros(n,1);x(n)=y(n)/U(n,n);fori=n-1:-1:1x(i)=(y(i)-U(i,i+1:n)*x(i+1:n))/U(i,i);end运行结果：4-3试验目的：考察Hilbert矩阵的性态，了解病态方程组求解。试验内容：生成Hilbert矩阵：Hn=hij,hij=1i+j-1,n=5,…,10。计算H（1）计算Hn算法实现：formatratforn=5:10fori=1:nforj=1:nH(i,j)=1/(i+j-1);endenddisp(H)cond(H,1)end运行结果：cond1H5=943656cond1Hcond1H7cond1H8cond1H9cond1H10矩阵的条件数越来越大，而且远远大于1，这是一个病态问题。（2）LU分解法算法实现：n=10fori=1:nforj=1:nH(i,j)=1/(i+j-1);endenddisp(H)cond(H,1)formatshort;n=size(H);L=eye(n);U(1,:)=H(1,:);fork=2:nifU(k-1,k-1)==0disp('分解失败');return;endL(k:n,k-1)=H(k:n,k-1)/U(k-1,k-1);U(k,k:n)=H(k,k:n)-L(k,1:k-1)*U(1:k-1,k:n);ifk<nH(k+1:n,k)=H(k+1:n,k)-L(k+1:n,1:k-1)*U(1:k-1,k);endendB([1:10],1)=1;Y=L\B;LUX=U\Y运行结果：（2）LU分解迭代求精法算法实现：n=10;fori=1:nforj=1:nH(i,j)=1/(i+j-1);endendformatshort;n=size(H);L=eye(n);U(1,:)=H(1,:);tol=10^(-3);fork=2:nifU(k-1,k-1)==0disp('分解失败');returnendL(k:n,k-1)=H(k:n,k-1)/U(k-1,k-1);U(k,k:n)=H(k,k:n)-L(k,1:k-1)*U(1:k-1,k:n);ifk<nH(k+1:n,k)=H(k+1:n,k)-L(k+1:n,1:k-1)*U(1:k-1,k);endendB([1:10],1)=1;Y=L\B;X=U\Y;r=B-H*X;Y=L\r;w=U\Y;formatlongwhilenorm(w)<tolr=B-H*X;Y=L\r;w=U\Y;X=X+w;endX运算结果：结果分析：对于同一复杂程度的方程组，LU分解迭代求精法所得的解和LU分解法所得到的解在精确度上相一致。5-1实验目的：熟悉Jacobi、Seidel、Sor迭代法，了解松弛因子对收敛速度的影响。实验内容：分别用Jacobi、Seidel、Sor迭代法求解下面的方程组，并作结果分析。初值，精度要求：（1）（2）Jacobi:A=[12.3,-2,-1,3.4,-3.7;1.4,9,-3,2.4,2.7;2.1,1,8,2.6,5.8;3.5,-2.1,1,13,4.6;2.5,-1,-2,5.3,14.8];b=[4.8;2.3;2.5;3.6;2.2];X0=[0;0;0;0;0];X=X0;K=1;whileK<=30fori=1:5X(i)=(b(i)-A(i,:)*X0)/A(i,i)+X0(i);endifnorm(X-X0)/X<0.00001disp(X);disp(K);return;endK=K+1;X0=X;end运行结果：Seidel:A=[12.3,-2,-1,3.4,-3.7;1.4,9,-3,2.4,2.7;2.1,1,8,2.6,5.8;3.5,-2.1,1,13,4.6;2.5,-1,-2,5.3,14.8];b=[4.8;2.3;2.5;3.6;2.2];X0=[0;0;0;0;0];X=X0;K=1;whileK<=20fori=1:5X(i)=(b(i)-A(i,:)*X)/A(i,i)+X(i);endifnorm(X-X0)/X<0.00001disp(X);disp(K);return;endK=K+1;X0=X;end运行结果：SOR:A=[12.3,-2,-1,3.4,-3.7;1.4,9,-3,2.4,2.7;2.1,1,8,2.6,5.8;3.5,-2.1,1,13,4.6;2.5,-1,-2,5.3,14.8];b=[4.8;2.3;2.5;3.6;2.2];X0=[0;0;0;0;0];X=X0;K=1;w=[0.8,1.1,1.2,1.3,1.4,1.5];forj=1:6m=w(j);whileK<=20fori=1:5X(i)=m*(b(i)-A(i,:)*X)/A(i,i)+X(i);endifnorm(X-X0,inf)/norm(X,inf)<0.00001disp(X);disp(K);break;endK=K+1;X0=X;endend运行结果：结果分析：1.Seidel迭代法收敛速度最快，迭代次数最少。2.当松弛因子越大，其收敛速度越慢，迭代次数越多。3.当松弛因子选取合适时（在1.0附近），Sor迭代法的收敛速度比Jacobi迭代法快。5-2实验目的：掌握Newton法与最速下降法求解非线性方程组，观察各自的优势。实验内容：（1）分别用Newton法和最速下降法求解下面非线性方程组，初值，精度要求。（2）改变初值，再用如上两种方法求解，得到什么结果？（3）采用初值，先用最速下降法求解3步，再用Newton迭代法求解，得到什么结果？对以上运算结果作分析。Newton法解非线性方程组程序：symsx1x2x3f1=3*x1-cos(x2*x3)-(1/2);f2=(x1)^2-81*(x2+0.1)^2+sin(x3)+1.06;f3=exp(-x1*x2)+20*x3+((10*pi-3)/3);g1=diff(f1,x1);g2=diff(f1,x2);g3=diff(f1,x3);g4=diff(f2,x1);g5=diff(f2,x2);g6=diff(f2,x3);g7=diff(f3,x1);g8=diff(f3,x2);g9=diff(f3,x3);x0=[0.1;0.1;-0.1];%x0=[20;20;20];Tol=10^(-5);F2=[subs(f1,{x1,x2,x3},{x0(1),x0(2),x0(3)});subs(f2,{x1,x2,x3},{x0(1),x0(2),x0(3)});subs(f3,{x1,x2,x3},{x0(1),x0(2),x0(3)})];F3=[subs(g1,{x1,x2,x3},{x0(1),x0(2),x0(3)})subs(g2,{x1,x2,x3},{x0(1),x0(2),x0(3)})subs(g3,{x1,x2,x3},{x0(1),x0(2),x0(3)});subs(g4,{x1,x2,x3},{x0(1),x0(2),x0(3)})subs(g5,{x1,x2,x3},{x0(1),x0(2),x0(3)})subs(g6,{x1,x2,x3},{x0(1),x0(2),x0(3)});subs(g7,{x1,x2,x3},{x0(1),x0(2),x0(3)})subs(g8,{x1,x2,x3},{x0(1),x0(2),x0(3)})subs(g9,{x1,x2,x3},{x0(1),x0(2),x0(3)})];k=1;z=-inv(F3);y=z*F2;whilek<=20x=x0+y;p=norm(x,inf);q=norm(x0,inf);if((p-q)/p<Tol)breakelsex0=x;F2=[subs(f1,{x1,x2,x3},{x0(1),x0(2),x0(3)});subs(f2,{x1,x2,x3},{x0(1),x0(2),x0(3)});subs(f3,{x1,x2,x3},{x0(1),x0(2),x0(3)})];F3=[subs(g1,{x1,x2,x3},{x0(1),x0(2),x0(3)})subs(g2,{x1,x2,x3},{x0(1),x0(2),x0(3)})subs(g3,{x1,x2,x3},{x0(1),x0(2),x0(3)});subs(g4,{x1,x2,x3},{x0(1),x0(2),x0(3)})subs(g5,{x1,x2,x3},{x0(1),x0(2),x0(3)})subs(g6,{x1,x2,x3},{x0(1),x0(2),x0(3)});subs(g7,{x1,x2,x3},{x0(1),x0(2),x0(3)})subs(g8,{x1,x2,x3},{x0(1),x0(2),x0(3)})subs(g9,{x1,x2,x3},{x0(1),x0(2),x0(3)})];z=-inv(F3);y=z*F2;x=x0+y;k=k+1;endenddisp(x);disp(k)最速下降法解非线性方程程序：symsx1x2x3tf1=3*x1-cos(x2*x3)-(1/2);f2=(x1)^2-81*(x2+0.1)^2+sin(x3)+1.06;f3=exp(-x1*x2)+20*x3+((10*pi-3)/3);F=(f1)^2+(f2)^2+(f3)^2;g1=diff(F,x1);g2=diff(F,x2);g3=diff(F,x3);G11=[g1;g2;g3];G1=-G11;x0=[0.1;0.1;-0.1];%x0=[20;20;20];G0=subs(G1,{x1,x2,x3},{x0(1),x0(2),x0(3)});z0=x0+(t*G0);F1=subs(F,{x1,x2,x3},{z0(1),z0(2),z0(3)});FD=diff(F1,t);xs0=solve(FD,t);FD_xs0=vpa(subs(FD,t,xs0),6);F1_xs0=vpa(subs(F1,t,xs0),6);t=xs0;k=1;Tol=10^(-5);whilek<=20x=x0+(t*G0);p=norm(x,inf);q=norm(x0,inf);if((p-q)/p<Tol)breakelsex0=x;G0=subs(G1,{x1,x2,x3},{x0(1),x0(2),x0(3)});x=x0+t*G0;k=k+1;endenddisp(x);disp(k)程序运行结果：（1）将初值定为分别运行上述两个程序，所得结果分别为：牛顿法：最速下降法：（2）改变初值为，分别运行上述两个程序，所得结果分别为：牛顿法：最速下降法：（3）采用初值，先用最速下降法求解3步，再用牛顿迭代法求解。首先将第一个程序运行3次后，将结果作为第二个程序的初值再做迭代，运行结果如下：结果分析：Newton法的收敛速度较快，并且用相邻两次迭代的无穷范数相对误差来控制迭代次数，由程序运行的结果可知，改变初值后Newton法迭代的结果发散，由此可见Newton法对初值的要求较高。而最速下降法，对初值的要求不高，它对任意的初值都是收敛的，因此，为了提高迭代的效率，应先用最速下降法迭代几次来确定初值，再用牛顿迭代法进行迭代。7-3试验目的：熟悉三次样条插值试验内容：利用固支三次样条函数画狗的轮廓图（上半部）。运行程序：x1=[125678101317];y1=[3.03.73.94.25.76.67.16.74.5];pp1=csape(x1,y1,'complete',[1,-2/3]);xx1=linspace(1,17,100);plot(xx1,ppval(pp1,xx1));holdonx2=[17202324252727.7];y2=[4.57.06.15.65.85.24.1];pp2=csape(x2,y2,'complete',[3,-4]);xx2=linspace(17,27.7,100);plot(xx2,ppval(pp2,xx2));holdonx3=[27.7282930];y3=[4.14.34.13.0];pp3=csape(x3,y3,'complete',[1/3,-1.5]);xx3=linspace(27.7,30,100);plot(xx3,ppval(pp3,xx3));gridon运行结果：8-1试验目的：熟悉最小二乘法拟合多项式。（1）程序：formatlongx=[0.40.550.650.800.901.05];y=[0.410750.578150.696750.888111.026521.25386];A=zeros(4,4);form=1:4forn=1:4fori=1:6A(m,n)=A(m,n)+x(i)^(m+n-2);endendendd=zeros(4,1);forn=1:4fori=1:6d(n,1)=d(n,1)+y(i)*x(i)^(n-1);endenda=A\d;R=0;fori=1:6R=R+(y(i)-(a(1)+a(2)*x(i)+a(3)*x(i)^2+a(4)*x(i)^3))^2;enddisp('系数a=');disp(a);disp('平方误差R=');disp(R);运行结果：（2）程序：functionb=Schmidt_orthogonalization(a)[m,n]=size(a);b=zeros(m,n);b(:,1)=a(:,1);fori=2:nforj=1:i-1b(:,i)=b(:,i)-dot(a(:,i),b(:,j))/dot(b(:,j),b(:,j))*b(:,j);endb(:,i)=b(:,i)+a(:,i);endfork=1:nb(:,k)=b(:,k)/norm(b(:,k));endx1=[0.40.550.650.800.901.05];y1=[0.410750.578150.696750.888111.026521.25386];x=Schmidt_orthogonalization(x1')y=Schmidt_orthogonalization(y1')formatlongA=zeros(4,4);form=1:4forn=1:4fori=1:6A(m,n)=A(m,n)+x(i)^(m+n-2);endendendd=zeros(4,1);forn=1:4fori=1:6d(n,1)=d(n,1)+y(i)*x(i)^(n-1);endenda=A\d;R=0;fori=1:6R=R+(y(i)-(a(1)+a(2)*x(i)+a(3)*x(i)^2+a(4)*x(i)^3))^2;enddisp('系数a=');disp(a);disp('平方误差R=');disp(R);运行结果：结果分析：比较(1)和（2）可以发现正交之后得到的最小二乘法的误差约为直接多项式拟合的五分之一，通过正交的pj(x)的函数系解决最小二乘法的病态结果问题，提供了求最小二乘法拟合多项式的快捷方法。9-1试验目的：熟悉数值积分公式，掌握数值计算定积分的方法。试验内容：采用不同方法数值计算积分。编写复化梯形公式和复化Simpson公式通用子程序，分别采用4,8，16,32,64等分区间计算。使用Romberg求积公式。（1）复化梯形公式：functionI_n=ComTrapezoidal(a,b,n)h=(b-a)/n;for(k=0:n)x(k+1)=a+k*h;if(x(k+1)==0)x(k+1)=10^(-10);endendI_1=h/2*(f(x(1))+f(x(n+1)));for(i=2:n)F(i)=h*f(x(i));e