计算方法上机题答案

2.用下列方法求方程e^x+10x-2=0的近似根，要求误差不超过5*10的负4次方，并比较计算量（1）二分法（局部，大图不太看得清，故后面两小题都用局部截图）（2）迭代法（3）牛顿法顺序消元法#includestdio.h#includestdlib.h#includemath.hintmain(){intN=4,i,j,p,q,k;doublem;doublea[4][5];doublex1,x2,x3,x4;for(i=0;iN;i++)for(j=0;jN+1;j++)scanf(%lf,&a[i][j]);for(p=0;pN-1;p++){for(k=p+1;kN;k++){m=a[k][p]/a[p][p];for(q=p;qN+1;q++)a[k][q]=a[k][q]-m*a[p][q];}}x4=a[3][4]/a[3][3];x3=(a[2][4]-x4*a[2][3])/a[2][2];x2=(a[1][4]-x4*a[1][3]-x3*a[1][2])/a[1][1];x1=(a[0][4]-x4*a[0][3]-x3*a[0][2]-x2*a[0][1])/a[0][0];printf(%f,%f,%f,%f,x1,x2,x3,x4);scanf(%lf,&a[i][j]);（这一步只是为了看到运行的结果）}运行结果列主元消元法function[x,det,flag]=Gauss(A,b)[n,m]=size(A);nb=length(b);flag='OK';det=1;x=zeros(n,1);fork=1:n-1max1=0;fori=k:nifabs(A(i,k))max1max1=abs(A(i,k));r=i;endendifmax11e-10flag='failure';return;endifrkforj=k:nz=A(k,j);A(k,j)=A(r,j);A(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;det=-det;endfori=k+1:nm=A(i,k)/A(k,k);forj=k+1:nA(i,j)=A(i,j)-m*A(k,j);endb(i)=b(i)-m*b(k);enddet=det*A(k,k);enddet=det*A(n,n)ifabs(A(n,n))1e-10flag='failure';return;endx(n)=b(n)/A(n,n);fork=n-1:-1：1forj=k+1:nb(k)=b(k)-A(k,j)*x(j);endx(k)=b(k)/A(k,k);end运行结果：雅可比迭代法functiony=jacobi(a,b,x0)D=diag(diag(a));U=-triu(a,1);L=-tril(a,-1);B=D\(L+U);f=D\b;y=B*x0+f;n=1;whilenorm(y-x0)1e-4x0=y;y=B*x0+f;n=n+1;endyn高斯赛德尔迭代法functiony=seidel(a,b,x0)D=diag(diag(a));U=-triu(a,1);L=-tril(a,-1);G=(D-L)\U;f=(D-L)\b;y=G*x0+f;n=1;whilenorm(y-x0)10^(-4)x0=y;y=G*x0+f;n=n+1;endynSOR迭代法functiony=sor(a,b,w,x0)D=diag(diag(a));U=-triu(a,1);L=-tril(a,-1);lw=(D-w*L)\((1-w)*D+w*U);f=(D-w*L)\b*w;y=lw*x0+f;n=1;whilenorm(y-x0)10^(-4)x0=y;y=lw*x0+f;n=n+1;endyn1.分段线性插值：functiony=fdxx(x0,y0,x)p=length(y0);n=length(x0);m=length(x);fori=1:mz=x(i);forj=1:n-1ifzx0(j+1)break;endendy(i)=y0(j)*(z-x0(j+1))/(x0(j)-x0(j+1))+y0(j+1)*(z-x0(j))/(x0(j+1)-x0(j));fprintf('y(%d)=%f\nx1=%.3fy1=%.3f\nx2=%.3fy2=%.3f\n\n',i,y(i),x0(j),y0(j),x0(j+1),y0(j+1));endend结果0.394040.380070.356932.分段二次插值：functiony=fdec(x0,y0,x)p=length(y0);n=length(x0);m=length(x);fori=1:mz=x(i);forj=1:n-1ifzx0(j+1)break;endendifj2j=j+1;elseif(jn-1)if(abs(x0(j)-z)abs(x0(j+1)-z))j=j+1;elseif((abs(x0(j)-z)==abs(x0(j+1)-z))&&(abs(x0(j-1)-z)abs(x0(j+2)-z)))j=j+1;endendans=0.0;fort=j-1:j+1a=1.0;fork=j-1:j+1ift~=ka=a*(z-x0(k))/(x0(t)-x0(k));endendans=ans+a*y0(t);endy(i)=ans;fprintf('y(%d)=%f\nx1=%.3fy1=%.3f\nx2=%.3fy2=%.3f\nx3=%.3fy3=%.3f\n\n',i,y(i),x0(j-1),y0(j-1),x0(j),y0(j),x0(j+1),y0(j+1));endend结果为0.394470.380220.357253.拉格朗日全区间插值functiony=lglr(x0,y0,x)p=length(y0);n=length(x0);m=length(x);fort=1:mans=0.0;z=x(t);fork=1:np=1.0;forq=1:nifq~=kp=p*(z-x0(q))/(x0(k)-x0(q));endendans=ans+p*y0(k);endy(t)=ans;fprintf('y(%d)=%f\n',t,y(t));endend结果为0.394470.380220.35722function[p,S,mu]=polyfit(x,y,n)if~isequal(size(x),size(y))errorendx=x(:);y=y(:);ifnargout2mu=[mean(x);std(x)];x=(x-mu(1))/mu(2);endV(:,n+1)=ones(length(x),1,class(x));forj=n:-1:1V(:,j)=x.*V(:,j+1);end[Q,R]=qr(V,0);ws=warning;p=R\(Q'*y);warning(ws);ifsize(R,2)size(R,1)warning;elseifwarnIfLargeConditionNumber(R)ifnargout2warning;elsewarning;endendifnargout1r=y-V*p;S.R=R;S.df=max(0,length(y)-(n+1));S.normr=norm(r);endp=p.';functionflag=warnIfLargeConditionNumber(R)ifisa(R,'double')flag=(condest(R)1e+10);elseflag=(condest(R)1e+05);endx=[1,3,4,5,6,7,8,9,10];y=[10,5,4,2,1,1,2,3,4];p=polyfit(x,y,2);y=poly2sym(p);a=-p(1)/p(2)*0.5;xa=-p(2)/(2*p(1));min=(4*p(1)*p(3)-p(2)^2)/(4*p(1));yxamin运行截图运行结果function[I]=CombineTraprl(f,a,b,eps)if(nargin==3)eps=1.0e-4;endn=1;h=(b-a)/2;I1=0;I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h;whileabs(I2-I1)epsn=n+1;h=(b-a)/n;I1=I2;I2=0;fori=0:n-1x=a+h*i;x1=x+h;I2=I2+(h/2)*(subs(sym(f),findsym(sym(f)),x)+subs(sym(f),findsym(sym(f)),x1));endendI=I2;程序：[q]=CombineTraprl('(1-exp(-x))^0.5/x'10^-12,1)结果：q=1.8521欧拉方法function[x,y]=euler(fun,x0,xfinal,y0,n);ifnargin5,n=50;endh=(xfinal-x0)/n;x(1)=x0;y(1)=y0;fori=1:nx(i+1)=x(i)+h;y(i+1)=y(i)+h*feval(fun,x(i),y(i));end程序：[x,y]=euler('doty',0,1,1,10)结果：改进欧拉方法function[x,y]=eulerpro(fun,x0,xfinal,y0,n);ifnargin5,n=50;endh=(xfinal-x0)/n;x(1)=x0;y(1)=y0;fori=1:nx(i+1)=x(i)+h;y1=y(i)+h*feval(fun,x(i),y(i));y2=y(i)+h*feval(fun,x(i+1),y1);y(i+1)=(y1+y2)/2;end程序：[x,y]=eulerpro('doty',0,1,1,10)结果：经典RK法function[x,y]=RungKutta4(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;forn=1:length(x)-1k1=feval(dyfun,x(n),y(n));k2=feval(dyfun,x(n)+h/2,y(n)+h/2*k1);k3=feval(dyfun,x(n)+h/2,y(n)+h/2*k2);k4=feval(dyfun,x(n+1),y(n)+h*k3);y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6;end程序：dyfun=inline('(2*x*y^-2)/3');[x,y]=RungKutta4(dyfun,[0,1],1,0.1)结果：标准函数x(1)=0;h=0.1;fori=1:10y(i)=(1+x(i)^2)^(1/3);endxY结果：