计算方法上机题答案

2.用下列方法求方程e^x+10x-2=0的近似根，要求误差不超过5*10的负4次方，并比较计算量（1）二分法（局部，大图不太看得清，故后面两小题都用局部截图）（2）迭代法（3）牛顿法顺序消元法#include<stdio.h>#include<stdlib.h>#include<math.h>int main(){ int N=4,i,j,p,q,k; double m;double a[4][5]; double x1,x2,x3,x4; for (i=0;i<N ;i++ )for (j=0;j<N+1; j++ )scanf("%lf",&a[i][j]);for(p=0;p<N-1;p++) {for(k=p+1;k<N;k++){m=a[k][p]/a[p][p];for(q=p;q<N+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);for k=1:n-1 max1=0;for i=k:nif abs(A(i,k))>max1max1=abs(A(i,k));r=i;endendif max1<1e-10flag='failure';return;endif r>kfor j=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; endfor i=k+1:nm=A(i,k)/A(k,k);for j=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)if abs(A(n,n))<1e-10flag='failure';return;endx(n)=b(n)/A(n,n);for k=n-1:-1：1for j=k+1:nb(k)=b(k)-A(k,j)*x(j);endx(k)=b(k)/A(k,k);end运行结果：雅可比迭代法function y=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;while norm(y-x0)>1e-4 x0=y;y=B*x0+f;n=n+1;endyn高斯赛德尔迭代法function y=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;while norm(y-x0)>10^(-4) x0=y;y=G*x0+f;n=n+1;endynSOR迭代法function y=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;while norm(y-x0)>10^(-4) x0=y;y=lw*x0+f;n=n+1;endyn1.分段线性插值：function y=fdxx(x0,y0,x)p=length(y0);n=length(x0);m=length(x);for i=1:mz=x(i);for j=1:n-1if z<x0(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.39404 0.38007 0.356932.分段二次插值：function y=fdec(x0,y0,x)p=length(y0);n=length(x0);m=length(x);for i=1:mz=x(i);for j=1:n-1if z<x0(j+1)break;endendif j<2j=j+1;elseif (j<n-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;for t=j-1:j+1a=1.0;for k=j-1:j+1if t~=ka=a*(z-x0(k))/(x0(t)-x0(k));endendans=ans+a*y0(t);endy(i)=ans;fprintf('y(%d)=%f\n x1=%.3f y1=%.3f\n x2=%.3f y2=%.3f\n x3=%.3f y3=%.3f\n\n',i,y(i),x0(j-1),y0(j-1),x0(j),y0(j),x0(j+1),y0(j+1));endend结果为0.39447 0.38022 0.357253.拉格朗日全区间插值function y=lglr(x0,y0,x)p=length(y0);n=length(x0);m=length(x);for t=1:mans=0.0;z=x(t);for k=1:np=1.0;for q=1:nif q~=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.39447 0.38022 0.35722function [p,S,mu] = polyfit(x,y,n)if ~isequal(size(x),size(y))errorendx = x(:);y = y(:);if nargout > 2mu = [mean(x); std(x)];x = (x - mu(1))/mu(2);endV(:,n+1) = ones(length(x),1,class(x));for j = n:-1:1V(:,j) = x.*V(:,j+1);end[Q,R] = qr(V,0);ws = warning;p = R\(Q'*y); warning(ws);if size(R,2) > size(R,1)warning;elseif warnIfLargeConditionNumber(R)if nargout > 2warning;elsewarning;endendif nargout > 1r = y - V*p;S.R = R;S.df = max(0,length(y) - (n+1));S.normr = norm(r);endp = p.';function flag = warnIfLargeConditionNumber(R) if isa(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; while abs(I2-I1)>epsn=n+1;h=(b-a)/n;I1=I2;I2=0;for i=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);if nargin<5,n=50;endh=(xfinal-x0)/n;x(1)=x0;y(1)=y0;for i=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); if nargin<5,n=50;endh=(xfinal-x0)/n;x(1)=x0;y(1)=y0;x(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;for n=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;for i=1:10y(i)=(1+x(i)^2)^(1/3); endxY结果：。