计算方法上机题答案

1、2.用下列方法求方程ex+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 a45; double x1,x2,x3,x4; for (i=0;i<N ;i+ ) for (j=0;j<N+1; j+ ) scanf("%lf&q

2、uot;,&aij); for(p=0;p<N-1;p+) for(k=p+1;k<N;k+) m=akp/app; for(q=p;q<N+1;q+) akq=akq-m*apq; x4=a34/a33; x3=(a24-x4*a23)/a22; x2=(a14-x4*a13-x3*a12)/a11; x1=(a04-x4*a03-x3*a02-x2*a01)/a00; printf("%f,%f,%f,%f",x1,x2,x3,x4); scanf("%lf",&aij); （这一步只是为了看到运行的结果） 运行结果

3、列主元消元法functionx,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:n if abs(A(i,k)>max1 max1=abs(A(i,k);r=i; end end if max1<1e-10 flag='failure'return;end if r>k 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)

4、;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-10 flag='failure'return;end x(n)=b(n)/A(n,n);for k=n-1:-1：1 for j=k+1:n b(k)=b(k)-A(k,j)*x(j); end x(k)=b(k)/A(k,k);end 运行结果

5、：雅可比迭代法function y=jacobi(a,b,x0)D=diag(diag(a);U=-triu(a,1);L=-tril(a,-1);B=D(L+U);f=Db;y=B*x0+f;n=1;while norm(y-x0)>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+

6、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 z<x0(j+1)

7、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(

8、i); for j=1:n-1 if z<x0(j+1) break; end end if j<2 j=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; 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

9、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)

10、; 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

11、 - 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*

12、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)

13、/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)>eps

14、    n=n+1;      h=(b-a)/n;    I1=I2;    I2=0;    for i=0:n-1            x=a+h*i;       

15、60;   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,xf

16、inal,y0,n);if nargin<5,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<5,n=50; endh=(xfinal-x0)/n; x(1)=x0;y(1)=y0; for i=1:n 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)-1 k1=feval(dyfun,x(n),y(n); k2=feval(