09级师兄数值分析上机报告.docx
- 文档编号:27491477
- 上传时间:2023-07-02
- 格式:DOCX
- 页数:23
- 大小:70.90KB
09级师兄数值分析上机报告.docx
《09级师兄数值分析上机报告.docx》由会员分享,可在线阅读,更多相关《09级师兄数值分析上机报告.docx(23页珍藏版)》请在冰豆网上搜索。
09级师兄数值分析上机报告
数值分析上机报告
1设
,
(1)由递推公式
,从
的几个近似值出发,计算
;
(2)粗糙估计
,用
,计算
;
(3)分析结果的可靠性及产生此现象的原因(重点分析原因)。
程序:
i=0.182322%i0的近似值
forn=1:
1:
20
i=-5*i+1/n
end
disp(i)%求得i20的值
%该程序称为程序一
i=0.009167%粗糙估计i20值
forn=20:
-1:
1
i=-i/5+1/5/n
end
disp(i)%求得i0的值
%该程序称为程序二
程序一结果为
i=0.182********000
.0883********
.0580********
i=0.04308333333333
i=0.03458333333334
.027*********
.0312********
i=-0.01339285714336
i=0.19196428571680
i=-0.84871031747287
i=4.34355158736433
i=-21.62684884591257
i=1.082175775628962e+002
i=-5.410109647375579e+002
i=2.705126252259218e+003
i=-1.352556459462942e+004
i=6.762788547314711e+004
i=-3.381393685422061e+005
i=1.690696898266586e+006
i=-8.453484438701350e+006
i=4.226742224350675e+007
4.226742224350675e+007
程序二结果
i=0.00916700000000
i=0.00816660000000
i=0.00889299578947
i=0.00933251195322
i=0.00989820349171
.010*********
i=0.01122926147300
i=0.01203986199111
i=0.01297664298639
i=0.01407133806939
i=0.01536755056794
i=0.01692648988641
i=0.01883692424494
.021*********
.024*********
.028*********
i=0.03430632955498
.0431********
.0580********
.0883********
i=0.182********395
0.182********395
%分析结果可知,程序一中虽然i0有相当高的精度,并且递推公式并无任何截断误差,但是
%i12~i20的计算结果却明显有很大误差。
程序二中虽然i20的精度不算很高,递推结果i2~i0
%却有相当高的精度。
即递推公式在理论上尽管完全等价,其实际计算效果却有很大不同。
%分析两种算法中误差的传播情况。
%设in有误差δn,在程序一的递推公式中,有δn=-5δn-1,这样δ20变得非常大。
而对于程序二中的递推公式,
%有δn-1=-δn/5,这样δ20在计算中逐渐减小。
故程序二结果可靠,而程序一不可靠。
2.求方程
的近似根,要求
,并比较计算量。
(1)在[0,1]上用二分法;
(2)取初值
,并用迭代
;
(3)加速迭代的结果;
(4)取初值
,并用牛顿迭代法;
(5)分析绝对误差。
(1)解:
二分法子程序【bisect.m】
function[c,err,yc]=bisect(f,a,b,delta)
%f是所要求解的函数
%a,b分别为有根区间的左右限
%delta是所允许的误差界
%c为所求近似解
%yc为函数f在c上的值
%err是c的误差估计
ifnargin<4delta=0.5*10^(-3);end
ya=feval('f',a);
yb=feval('f',b);
ifyb==0,c=b,return,end
ifya*yb>0
disp('(a,b)不是有根区间');
return,end
max1=1+round((log(b-a)-log(delta))/log
(2));
fork=1:
max1
c=(a+b)/2
yc=feval('f',c);
ifyc==0a=c;b=c;break,
elseifyb*yc>0
b=c;yb=yc;
else
a=c;ya=c;
end
if(b-a) end k,c=(a+b)/2,err=abs(b-a),yc=feval('f',c) 函数【f】 functiony=f(x) y=exp(x)+10*x-2; 主程序【homework2_1.m】 clear clc bisect('f',0,1); 运行结果: k=%次数 11 c=%近似根 0.0906 err=%误差 4.8828e-004 yc= 5.6662e-004 (2)使用迭代公式主程序: fixpt('g',0,0.5*10^(-3),200); 【fixpt.m】 function[p0,k,err,p]=fixpt(g,p0,tol,max1) %g是给定的迭代函数 %p0是给定的初值 %max1是所允许的最大迭代次数 %k是迭代次数加1 %p是不动点的近似值 %err是误差 %p(p1,p2,...pn) P (1)=p0; fork=2: max1 P(k)=feval('g',P(k-1)); k,err=abs(P(k)-P(k-1)) p=P(k); if(err break;end ifk==max1 disp('超过最大迭代次数,迭代失败! '); end end P 其中【g.m】 functiony=g(x) y=(2-exp(x))/10; 运行结果: k= 5%迭代次数 err= 1.2652e-004%误差 P= 00.10000.08950.09060.0905%最后一个为迭代最终结果。 (3)牛顿法: 主程序【homework2_3.m】 newton('f','df',0,0.5*10^(-3),100); 其中【newton.m】 function[p1,err,k,y]=newton(f,df,p0,delta,max1) %f是非线性函数 %df是f的微商 %p0是给定的初值 %max1是所允许的最大迭代次数 %k是迭代次数加1 %p是不动点的近似值 %err是误差 %p1是牛顿法求得的方程的近似值 %y=f(p1) p0,feval('f',p0) fork=1: max1 p1=p0-feval('f',p0)/feval('df',p0) err=abs(p1-p0); p0=p1; p1,err,k,y=feval('f',p1) if(err break;end end p1,err,k,y=feval('f',p1) 【df.m】 functiony=df(x) y=exp(x)+10; 迭代结果: p1= 0.0905%结果点 err= 3.8398e-004%误差 k= 2%迭代次数 3.钢水包使用次数多以后,钢包的容积增大,数据如下: x 2 3 4 5 6 7 8 9 y 6.42 8.2 9.58 9.5 9.7 10 9.93 9.99 10 11 12 13 14 15 16 10.49 10.59 10.60 10.8 10.6 10.9 10.76 试从中找出使用次数和容积之间的关系,计算均方差。 (注: 增速减少,用何种模型) 设y=f(x)具有指数形式 (a>0,b<0)。 对此式两边取对数,得 。 记A=lna,B=b,并引入新变量z=lny,t=1/x。 引入新变量后的数据表如下 x 2 3 4 5 6 7 8 9 t=1/x 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 z=lny 1.8594 2.1041 2.2597 2.2513 2.2721 2.3026 2.2956 2.3016 10 11 12 13 14 15 16 0.1000 0.0909 0.0833 0.0769 0.0714 0.0667 0.0625 2.3504 2.3599 2.3609 2.3795 2.3609 2.3888 2.3758 程序: t=[0.50000.33330.25000.20000.16670.14290.12500.11110.10000.09090.08330.07690.07140.06670.0625]; z=[1.85942.10412.25972.25132.27212.30262.29562.30162.35042.35992.36092.37952.36092.38882.3758]; polyfit(t,z,1) 运行结果: ans= -1.11072.4578 由此可得A=2.4578,B=-1.1107, ,b=B=-1.1107 方程即为 4.设 , , 分析下列迭代法的收敛性,并求 的近似解及相应的迭代次数。 (1)JACOBI迭代; (2)GAUSS-SEIDEL迭代; (3)SOR迭代( )。 程序如下: function[x,k,flag,err]=jacobi(A,b,delta,maxl) %求解线性方程组的迭代法,其中 %A为方程组的系数矩阵 %b为方程组的右端项 %delta为代数精度,默认值为1e-5 %maxl为最大迭代次数,默认值为100 %x为方程组的解; %k为迭代次数 %flag为指标变量flag='OK! '表示迭代收敛指标要求 %flag='fail'表示迭代失败 ifnargin<4 maxl=100; end ifnargin<3 delta=1e-5; end n=length(A); k=0; x=zeros(n,1); y=zeros(n,1); flag='OK! '; whilek<99 fori=1: n y(i)=b(i); forj=1: n ifj~=i y(i)=y(i)-A(i,j)*x(j); end end %ifabs(A(i,i)<1e-10)|k==maxl %flag='Fail'; %return %end y(i)=y(i)/A(i,i); end ifnorm(y-x,inf) break end x=y k=k+1 end 计算结果: x= 1.0000 1.9999 0.9999 1.9999 0.9999 1.9999 k=26 G-S法: function[x,k,flag]=gau_seid(A,b,delta,maxl) %求解线性方程组的迭代法,其中 %A为方程组的系数矩阵 %b为方程组的右端项 %delta为代数精度,默认值为1e-5 %maxl为最大迭代次数,默认值为100 %x为方程组的解; %k为迭代次数 %flag为指标变量flag='OK! '表示迭代收敛指标要求 %flag='fail'表示迭代失败 ifnargin<4 maxl=100; end ifnargin<3 delta=1e-5; end n=length(A); k=0; x=zeros(n,1); y=zeros(n,1); flag='OK! '; whilek<100 y=x; fori=1: n z=b(i); forj=1: n ifj~=i z=z-A(i,j)*x(j); end end %ifabs(A(i,i)<1e-10)|k==maxl %flag='Fail'; %return %end z=z/A(i,i); x(i)=z; end ifnorm(y-x,inf) break end k=k+1 end x= 0.9999 1.9999 0.9999 1.9999 0.9999 2.0000 k=13 SOR法: function[x,k,flag]=sor(A,b,delta,w,maxl) %求解线性方程组的迭代法,其中 %A为方程组的系数矩阵 %b为方程组的右端项 %delta为代数精度,默认值为1e-5 %w为松弛因子 %maxl为最大迭代次数,默认值为100 %x为方程组的解; %k为迭代次数 %flag为指标变量flag='OK! '表示迭代收敛指标要求 %flag='fail'表示迭代失败 ifnargin<5 maxl=100; end ifnargin<4 w=1; end ifnargin<3 delta=1e-5; end n=length(A); k=0; x=zeros(n,1); y=zeros(n,1); flag='OK! '; whilek<100 y=x; fori=1: n z=b(i); forj=1: n ifj~=i z=z-A(i,j)*x(j); end end %ifabs(A(i,i)<1e-10)|k==maxl %flag='Fail'; %return %end z=z/A(i,i); x(i)=(1-w)*x(i)+w*z; end ifnorm(y-x,inf) break end k=k+1 end 计算结果: 当w=1.334时x= 1.0000 2.0000 1.0000 2.0000 1.0000 2.0000 k=12 当w=1.95时x= 1.0071 2.0029 1.0140 1.9899 1.0095 2.0071 k=100 当w=0.95时x= 0.9999 1.9999 0.9999 1.9999 0.9999 2.0000 k=15 5用逆幂迭代法求 最接近于11的特征值和特征向量,准确到 。 先建立M文件pow_inv.m function[m,u,flag]=pow_inv(A,delta,max1) I=[100;010;001]; ifnargin<3max1=100;end ifnargin<2delta=1e-3;end n=length(A);u=ones(n,1); flag='Fail! ';k=0;m1=0; invA=inv(A-11*I); whilek<=max1 v=invA*u;[vmax,i]=max(abs(v)); m=v(i);u=v/m; ifabs(m-m1) flag='OK! ';break; end m1=m;k=k+1; end m=11+1/m; 然后在窗口输入: A=[631;321;111]; [m,u,flag]=pow_inv(A) 得到结果为 m=7.8745 u=1.0000 0.5495 0.2259 flag=OK! 6.用经典R-K方法求解初值问题 (1) , , ; (2) , , 。 和精确解 比较,分析结论。 程序: function[x,y,z]=nark4(dyfun1,dyfun2,xspan,y0,z0,h) %y'=f(x,y,z),z'=g(x,y,z),y(x0)=y0,z(x0)=z0 %dyfun为函数f(x,y),xspan为求解区间[x0,xN],y0为初值y(x0),z0为初值z(x0),h为步长 %x返回节点,y返回数值解 x=xspan (1): h: xspan (2); y (1)=y0; z (1)=z0; forn=1: length(x)-1 k1=feval(dyfun1,x(n),y(n),z(n)); l1=feval(dyfun2,x(n),y(n),z(n)); k2=feval(dyfun1,x(n)+h/2,y(n)+h/2*k1,z(n)+h/2*l1); l2=feval(dyfun2,x(n)+h/2,y(n)+h/2*k1,z(n)+h/2*l1); k3=feval(dyfun1,x(n)+h/2,y(n)+h/2*k2,z(n)+h/2*l2); l3=feval(dyfun2,x(n)+h/2,y(n)+h/2*k2,z(n)+h/2*l2); k4=feval(dyfun1,x(n+1),y(n)+h*k3,z(n)+h*l3); l4=feval(dyfun2,x(n+1),y(n)+h*k3,z(n)+h*l3); y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6; z(n+1)=z(n)+h*(l1+2*l2+2*l3+l4)/6; end (1)其中: dyfun1=inline('-2*y+z+2*sin(x)'); dyfun2=inline('y-2*z+2*cos(x)-2*sin(x)'); 结果: x00.50001.00001.50002.00002.50003.0000 y2.00001.69621.57901.44191.17470.75530.2331 z3.00002.08751.27480.5192-0.1400-0.6293-0.8824 x3.50004.00004.50005.00005.50006.00006.5000 y-0.2963-0.7229-0.9541-0.9406-0.6901-0.26650.2249 z-0.8695-0.6135-0.18890.29310.71010.95730.9726 x7.00007.50008.00008.50009.00009.500010.0000 y0.66270.93910.98620.79220.4044-0.0823-0.5487 z0.75140.3470-0.1417-0.5954-0.9031-0.9896-0.8337 和真实值比较x=10时真实值y=-0.5439,z=-0.8390.误差分别为0.48,0.47 (2)dyfun1=inline('-2*y+z+2*sin(x)'); dyfun2=inline('998*y-999*z+999*cos(x)-999*sin(x)'); 1.0e+184* xyz 00.00000.0000 0.00000.0000-0.0000 0.00000.0000-0.0000 0.00000.0000-0.0000 0.00000.0000-0.0000 0.00000.0000-0.0000 0.00000.0000-0.0000 0.00000.0000-0.0000 0.00000.0000-0.0000 0.00000.0000-0.0000 0.00000.0000-0.0000 0.00000.0000-0.0000 0.00000.0000-0.0000 0.00000.0000-0.0000 0.00000.0000-0.0000 0.00000.0000-0.0000 0.00000.0000-0.0000 0.00000.0000-0.0000 0.00000.0000-0.0000 0.00000.0000-0.0000 0.00000.0020-2.0121 7.用有限差分法求解边值问题(h=0.1): . x (1)=-1; >>h=0.1; >>N=21; >>y (1)=1; >>y(N)=1; >>fori=2: N x(i)=x (1)+(i-1)*h; q(i)=-(2+h^2*(1+x(i))^2); a(i-1)=q(i); b(i-1)=1; c(i-1)=1; end >>n=length(a); >>A=diag(a)+diag(b(1: n-1),1)+diag(c(2: n),-1); >>B (1)=-1; >>B(N-1)=-1; >>fori=2: N-2 B(i)=0; end >>y=A\B' y= 0.9524 0.9050 0.8578 0.8115 0.7664 0.7233 0.6828 0.6456 0.6126 0.5845 0.5622 0.5468 0.5392 0.5408 0.5529 0.5775 0.6169 0.6741 0.7532 0.8594
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 09 师兄 数值 分析 上机 报告
![提示](https://static.bdocx.com/images/bang_tan.gif)