计算方法上机题答案.docx
- 文档编号:25945888
- 上传时间:2023-06-16
- 格式:DOCX
- 页数:21
- 大小:545.60KB
计算方法上机题答案.docx
《计算方法上机题答案.docx》由会员分享,可在线阅读,更多相关《计算方法上机题答案.docx(21页珍藏版)》请在冰豆网上搜索。
计算方法上机题答案
2.用下列方法求方程e^x+10x-2=0的近似根,要求误差不超过5*10的负4次方,并比较计算量
(1)二分法
(局部,大图不太看得清,故后面两小题都用局部截图)
(2)迭代法
(3)牛顿法
顺序消元法
#include
#include
#include
intmain()
{intN=4,i,j,p,q,k;
doublem;
doublea[4][5];
doublex1,x2,x3,x4;
for(i=0;i for(j=0;j scanf("%lf",&a[i][j]); for(p=0;p { for(k=p+1;k { m=a[k][p]/a[p][p]; for(q=p;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: n ifabs(A(i,k))>max1 max1=abs(A(i,k));r=i; end end ifmax1<1e-10 flag='failure';return; end ifr>k forj=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 fori=k+1: n m=A(i,k)/A(k,k); forj=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) ifabs(A(n,n))<1e-10 flag='failure';return; end x(n)=b(n)/A(n,n); fork=n-1: -1: 1 forj=k+1: n b(k)=b(k)-A(k,j)*x(j); end x(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-4 x0=y; y=B*x0+f;n=n+1; end y n 高斯赛德尔迭代法 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; end y n SOR迭代法 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; end y n 1.分段线性插值: functiony=fdxx(x0,y0,x) p=length(y0);n=length(x0);m=length(x); fori=1: m z=x(i); forj=1: n-1 ifz 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)=%f\nx1=%.3fy1=%.3f\nx2=%.3fy2=%.3f\n\n',i,y(i),x0(j),y0(j),x0(j+1),y0(j+1)); end end 结果0.394040.380070.35693 2.分段二次插值: functiony=fdec(x0,y0,x) p=length(y0);n=length(x0);m=length(x); fori=1: m z=x(i); forj=1: n-1 ifz break; end end ifj<2 j=j+1; elseif(j 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; fort=j-1: j+1 a=1.0; fork=j-1: j+1 ift~=k a=a*(z-x0(k))/(x0(t)-x0(k)); end end ans=ans+a*y0(t); end y(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)); end end 结果为0.394470.380220.35725 3.拉格朗日全区间插值 functiony=lglr(x0,y0,x) p=length(y0);n=length(x0);m=length(x); fort=1: m ans=0.0; z=x(t); fork=1: n p=1.0; forq=1: n ifq~=k p=p*(z-x0(q))/(x0(k)-x0(q)); end end ans=ans+p*y0(k); end y(t)=ans; fprintf('y(%d)=%f\n',t,y(t)); end end 结果为0.394470.380220.35722 function[p,S,mu]=polyfit(x,y,n) if~isequal(size(x),size(y)) error end x=x(: ); y=y(: ); ifnargout>2 mu=[mean(x);std(x)]; x=(x-mu (1))/mu (2); end V(: n+1)=ones(length(x),1,class(x)); forj=n: -1: 1 V(: 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) ifnargout>2 warning; else warning; end end ifnargout>1 r=y-V*p; S.R=R; S.df=max(0,length(y)-(n+1)); S.normr=norm(r); end p=p.'; functionflag=warnIfLargeConditionNumber(R) ifisa(R,'double') flag=(condest(R)>1e+10); else flag=(condest(R)>1e+05); end x=[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)); y xa min 运行截图 运行结果 function [I] = CombineTraprl(f,a,b,eps) if(nargin ==3) eps=1.0e-4; end n=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 n=n+1; h=(b-a)/n; I1=I2; I2=0; for i=0: n-1 x=a+h*i; x1=x+h; I2=I2+(h/2)*(subs(sym(f),findsym(sym(f)),x)+subs(sym(f),findsym(sym(f)),x1)); end end I=I2; 程序: >> [q]=CombineTraprl('(1-exp(-x))^0.5/x'10^-12,1) 结果: q = 1.8521 欧拉方法 function[x,y]=euler(fun,x0,xfinal,y0,n); ifnargin<5,n=50; end h=(xfinal-x0)/n; x (1)=x0;y (1)=y0; fori=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); ifnargin<5,n=50; end h=(xfinal-x0)/n; x (1)=x0; y (1)=y0; fori=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; forn=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=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: 10 y(i)=(1+x(i)^2)^(1/3); end x Y 结果:
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 计算方法 上机 答案