丁軍計算方法作業(yè)北京科技大學(xué)

第一次作業(yè)3-1(1)2x-ex+3=0采用了3種形式g(x),m文件寫法相同：m文件formatlongp0=-1;N=100;k=1;whilek<=NP=g(p0);k=k+1;p0=p;enddisp(p);——函數(shù)functiony=g(x)y=(exp(x)-3)/2;end結(jié)果(p0=1) p=-1.373374545351944(p0=-1)p=-1.373374545351944函數(shù)functiony=g(x)y=exp(x)-3-x;end結(jié)果(p0=1)p=-1.373374545351969(p0=-1)p=-1.373374545351871函數(shù)functiony=g(x)y=log(2*x+3);end結(jié)果(p0=1) p=1.923938750345880(p0=-1)p=1.923938750345880(2)3x5+5x2-10=0采用了3種形式g(x),m文件寫法相同：M文件formatlongp0=1;N=100;k=1;whilek<=Np=g(p0)k=k+1;p0=p;enddisp(p);g(x)=3x5+5x2+x-10函數(shù)functiony=f(x)y=3*x.A5+5*x.A2+x-10;end結(jié)果不收斂 函數(shù)functiony=f2(x)y=(10-5*x.A2)/(3*x.A4);end結(jié)果不收斂 函數(shù)functiony=f3(x)y=((10-5*x.A2)/3).A0.2;end結(jié)果p=1.0722276979308363-5(1)二分法函數(shù)functiony=f(x)y=600*x.A4-550*x.A3+200*x.A2-20*x-1;endM文件a=0.1;b=1;N=100;tol=0.5*10A(-4);k=0;fa=f(a);fork=1:Np=(a+b)/2;fp=f(p);if(fp==0|(b-a)/2<tol)breakendiffa*fp<0b=p;elsea=p;endenddisp(p);disp(k)結(jié)果P=0.232357788085938K=15(2)牛頓迭代法函數(shù)functiony=f(x)y=600*x.A4-550*x.A3+200*x.A2-20*x-1;endfunctiony=df(x)y=2400*x.A3-1650*x.A2+400*x-20;endM文件p0=0.1;N=100;tol=0.5*10A(-4);fork=1:Np1=p0-f(p0)/df(p0);ifabs(p1-p0)<tolbreakendp0=p1;enddisp(p);disp(k)結(jié)果P=0.232352964749985K=5(3)牛頓下山法M文件p0=0.1；N=100;tol=0.5*10A(-4);e=10A-8;fork=1:Nw=1;p1=p0-w*f(p0)/df(p0);while(abs(f(p1))>=abs(f(p0))&w>e)w=w/2;p1=p0-w*f(p0)/df(p0);endifabs(p1-p0)<tolbreakendp0=p1;enddisp(p1);disp(k)結(jié)果P=0.232352965423057K=3(4)割線法M文件p0=0.1;p1=1;N=100;tol=0.5*10A(-4);fork=1:Np=p1-f(p1)*(p1-p0)/(f(p1)-f(p0));ifabs(p-p1)<tolbreakendp0=p1;p1=p;enddisp(p);disp(k)結(jié)果P=0.232352956733991K=7第二次作業(yè)4-2代碼A=[1.3,-1,-1,2.4,-3.4;2.4,-1,-1,1.4,3.7;2,1,-2,3.6,6.8;2.5,-1,4,3,6.6;1.5,-1,-1,5.3,2.8];b=[4.2;6.3;5.5;3.6;6.2];L=eye(5);U(1,:)=A(1,:);y=[0;0;0];x=[0;0;0];fork=2:5ifU(k-1,k-1)==0disp('分解失敗');ReturnendL(k:5,k-1)=A(k:5,k-1)/U(k-1,k-1);U(k,k:5)=A(k,k:5)-L(k,1:k-1)*U(1:k-1,k:5);ifk<5A(k+1:5,k)=A(k+1:5,k)-L(k+1:5,1:k-1)*U(1:k-1,k);endendL1=inv(L);y=L1*b;U1=inv(U);x=U1*y;LUX結(jié)果1.000000001.8462L00000001.53853,00001,0000001.9231L0909-L5667L000001.15380,18180Q1846L00001.3000-L0000-L000G工4000一風(fēng)4QQQ0Q,S462d3462-3.03089.97690QT.000090000-LL90000QQ6909-27.57880Q00LQ00131.5128-1.0257-D.7460D.3490D.L1064-3代碼formatratforn=5:10fori=1:nforj=1:nH(i,j)=1/(i+j-1);endendcond(H,1)end結(jié)果=943656=29070279=985194887=33872788559=1099649476470=35352333500164矩陣的條件數(shù)越來越大，而且遠遠大于1,這是一個病態(tài)問題(1)應(yīng)用LU分解法求解代碼formatlong;n=10fori=1:nforj=1:nH(i,j)=1/(i+j-1);endendn=size(H);L=eye(n);U(1,:)=H(1,:);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;L1=inv(L);y=L1*b;U1=inv(U);x=U1*y;LUx結(jié)果',Oe+06*-0.0000099964952120.000989693993574-0.023753430904197a240179914949146260972005136S533.782965146456109-6.725411S69085a887.000012023760501-3.93^527372606900.923632957130733(2)應(yīng)用LU分解迭代求精法求解代碼clearformatlong;n=10;fori=1:nforj=1:nH(i,j)=1/(i+j-1);endendn=size(H);L=eye(n);U(1,:)=H(1,:);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;L1=inv(L);y=L1*b;U1=inv(U);x=U1*y;r=b-H*x;y=L1*r;w=U1*y;tol=10A(12)whilenorm(w)<tolr=b-H*x;y=L1*r;w=U1*y;x=x+w;endx5-1(1)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<=300fori=1:5X(i)=(b(i)-A(i,:)*X0)/A(i,i)+X0(i);endifnorm(X-X0,inf)/norm(X,inf)<1*10A(-5)XKreturn;endK=K+1;X0=X;End結(jié)果0.3M60.16780.09960.17540.044717Seidel:clearclcA=[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<=300fori=1:5X⑴二(b(i)-A(i,:)*X)/A(i,i)+X⑴；endifnorm(X-X0,inf)/norm(X,inf)<1*10A(-5)XKreturn;endK=K+1;X0=X;enddisp(發(fā)散')；結(jié)果0.39060.16730.09960.17540.0447K=SOR:代碼clc;clearA=[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=zeros(5,1);w=0.8;%1.1;1.2;1.3;1.4;1.5;k=1;x=x0;whilek<=300fori=1:5x(i)=w*(b(i)-A(i,:)*x)/A(i,i)+x(i);endifnorm(x-x0,inf)/norm(x,inf)<10A(-5)wxkreturnendk=k+1;x0=x;enddisp(發(fā)散')；結(jié)果w=w二w=w=W=0.8000L1000L2000L.3000L4000X=x-K—K一x-0.39060.39060.39060.39060.39060.16780.16780.1&7S0.16780.16780.09960.0996CL09950.09960.0996W=0.17540.17j40.17540.1754ClL7540.04470.94470.04470.04470.04471.5000k-k-k=k-K-發(fā)散二12142565(2)Jacobi:代碼clearclcA=[13.3,-4,-1,3.5,-3.8;3.4,9,-3,4.4,2.3;4.1,1,7,2.7,5.9;2.5,-2.4,1,13,5.6;1.5,-1,-3,4.3,14.9];b=[5.8;4.3;2.6;3.8;4.2];X0=[0;0;0;0;0];X=X0;K=1;whileK<=300fori=1:5X(i)=(b⑴-A(i,:)*X0)/A(i,i)+X0⑴；endifnorm(X-X0,inf)/norm(X,inf)<1*10A(-5)XKreturn;endK=K+1;X0=X;End結(jié)果Q.47200.-.330-0.13020.16280.1701K=25Seidel:clearclcA=[13.3,-4,-1,3.5,-3.8;3.4,9,-3,4.4,2.3;4.1,1,7,2.7,5.9;2.5,-2.4,1,13,5.6;1.5,-1,-3,4.3,14.9];b=[5.8;4.3;2.6;3.8;4.2];X0=[0;0;0;0;0];X=X0;K=1;whileK<=300fori=1:5X⑴二(b(i)-A(i,:)*X)/A(i,i)+X⑴；endifnorm(X-X0,inf)/norm(X,inf)<1*10A(-5)XKreturn;endK=K+1;X0=X;enddisp(發(fā)散')；結(jié)果0.47200.1330-0.l：020.16280.1701K-12SOR:代碼clearclcA=[13.3,-4,-1,3.5,-3.8;3.4,9,-3,4.4,2.3;4.1,1,7,2.7,5.9;2.5,-2.4,1,13,5.6;1.5,-1,-3,4.3,14.9];b=[5.8;4.3;2.6;3.8;4.2];x0=zeros(5,1);w=0.8;%1.1;1.2;1.3;1.4;1.5;k=1;x=x0;whilek<=300fori=1:5x(i)=w*(b(i)-A(i,:)*x)/A(i,i)+x(i);endifnorm(x-x0,inf)/norm(x,inf)<10A(-5)wxkreturnendk=k+1;x0=x;enddisp（發(fā)散'）；結(jié)果叫=w-w-w一0.SOOO1.10001.2000[.3000003X=0.4720K-0.4720x=0.47200.13300.13300.13300.1330-Q.1902-0.1302-aisos-0.1302llj=0.16230.16230.16280.L62Snrw=0.17010.17010.17010.17011.40001.5000k=k=k=k=10182S64發(fā)散發(fā)散兩道題的結(jié)果表明，三種迭代方法都能正確求得方程組，Jacobi迭代法的迭代次數(shù)較多，Seidel的迭代次數(shù)較少。對于Sor迭代法，松弛因子w較小時，迭代次數(shù)較少。當松弛因子w過大時，結(jié)果發(fā)散，無法求解方程組。5-2(1)Newton法代碼functionf=F(x)f(1)=3*x(1)-cos(x(2)*x(3))-0.5;f(2)=x(1)A2-81*(x(2)+0.1)A2+sin(x(3))+1.06;f(3)=exp(-x(1)*x(2))+20*x(3)+(10*pi-3)/3;f=[f(1);f(2);f(3)];functiondf=DF(x)df=[3,x(3)*sin(x(2)),x(2)*sin(x(3));2*x(1),-162*(x(2)+0.1),cos(x(3));-x(2)*exp(-x(1)*x(2)),-x(1)*exp(-x(1)*x(2)),20];clc;clearx0=[0.1;0.1;-0.1];fork=1:300y=-DF(x0)\F(x0);x=x0+y;ifnorm(x-x0,inf)/norm(x,inf)<1*10A(-5)xkbreak;endx0=x;end結(jié)果0.500000000000000CL000000009000000-0,523508775598299=最速下降法代碼clc;clearsymsx1x2x3tf1=3*x1-cos(x2*x3)-(1/2);f2=(x1)A2-81*(x2+0.1)A2+sin(x3)+1.06;f3=exp(-x1*x2)+20*x3+((10*pi-3)/3);F=(f1)A2+(f2)A2+(f3)A2;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];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=so