计算方法上机题问题详解.doc

标准文档2.用下列方法求方程ex+10x-2=0的近似根，要求误差不超过5*10的负4次方，并比较计算量（1）二分法（局部，大图不太看得清，故后面两小题都用局部截图）（2）迭代法（3）牛顿法顺序消元法#include#include#includeint main() int N=4,i,j,p,q,k; double m; double a45; double x1,x2,x3,x4; for (i=0;iN ;i+ ) for (j=0;jN+1; j+ ) scanf(%lf,&aij); for(p=0;pN-1;p+) for(k=p+1;kN;k+) m=akp/app; for(q=p;qmax1 max1=abs(A(i,k);r=i; end end if max1k for j=k:n z=A(k,j);A(k,j)=A(r,j);A(r,j)=z; end z=b(k);b(k)=b(r);b(r)=z;det=-det; end for i=k+1:n m=A(i,k)/A(k,k); for j=k+1:n A(i,j)=A(i,j)-m*A(k,j); end b(i)=b(i)-m*b(k); end det=det*A(k,k); end det=det*A(n,n) if abs(A(n,n)1e-4x0=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:m z=x(i); for j=1:n-1 if zx0(j+1) break; end end y(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)=%fnx1=%.3fy1=%.3fnx2=%.3fy2=%.3fnn,i,y(i),x0(j),y0(j),x0(j+1),y0(j+1); end end 结果0.39404 0.38007 0.356932.分段二次插值：function y=fdec(x0,y0,x)p=length(y0);n=length(x0);m=length(x); for i=1:m z=x(i); for j=1:n-1 if zx0(j+1) break; end end if j2 j=j+1; elseif (jabs(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; end end ans=0.0; for t=j-1:j+1 a=1.0; for k=j-1:j+1 if t=k a=a*(z-x0(k)/(x0(t)-x0(k); end end ans=ans+a*y0(t); end y(i)=ans; fprintf(y(%d)=%fn x1=%.3f y1=%.3fn x2=%.3f y2=%.3fn x3=%.3f y3=%.3fnn,i,y(i),x0(j-1),y0(j-1),x0(j),y0(j),x0(j+1),y0(j+1);end end 结果为0.39447 0.38022 0.357253.拉格朗日全区间插值function y=lglr(x0,y0,x) p=length(y0);n=length(x0);m=length(x); for t=1:m ans=0.0; z=x(t); for k=1:n p=1.0; for q=1:n if q=k p=p*(z-x0(q)/(x0(k)-x0(q); end end ans=ans+p*y0(k); end y(t)=ans; fprintf(y(%d)=%fn,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 2 mu = mean(x); std(x); x = (x - mu(1)/mu(2);endV(:,n+1) = ones(length(x),1,class(x);for j = n:-1:1 V(:,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 2 warning; else warning; endendif nargout 1 r = 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);else flag = (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运行截图运行结果functionI=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/x10-12,1)结果：q=1.8521欧拉方法function x,y=euler(fun,x0,xfinal,y0,n);if nargin5,n=50;endh=(xfinal-x0)/n; x(1)=x0;y(1)=y0;for i=1:n x(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 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)-1 k1=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