我的数值计算上机报告上海电力学院.docx
- 文档编号:4686534
- 上传时间:2022-12-07
- 格式:DOCX
- 页数:25
- 大小:138.86KB
我的数值计算上机报告上海电力学院.docx
《我的数值计算上机报告上海电力学院.docx》由会员分享,可在线阅读,更多相关《我的数值计算上机报告上海电力学院.docx(25页珍藏版)》请在冰豆网上搜索。
我的数值计算上机报告上海电力学院
数值计算方法上机
实习报告
专业年级:
学生姓名:
学号:
2012年12月25日
1.设
,
(1)由递推公式
,从
的几个近似值出发,计算
;
(2)给出
的近似值,用
,计算
;
(3)分析结果的可靠性及产生此现象的原因(重点分析原因)。
解:
(1)
ln6-ln5=0.18232,
将
代入递推公式
中,MATLAB程序如下:
I1=0.18232;
I2=0.1823;
I3=0.18
forn=1:
1:
20
I1=-5*I1+1/n;
I2=-5*I2+1/n;
I3=-5*I3+1/n;
fprintf('%.4f\n',I1,I2,I3)
end
I1
I2
I3
结果为:
I1=-1.4847e+008;I2=-2.0558e+009;I3=-2.2140e+011
分析:
精度不同时计算出来的结果是不同的
(2)
>>i=0.009167%粗糙估计i20值
i=
0.0092
>>forn=20:
-1:
1
i=-i/5+1/5/n;
end
>>disp(i)
0.1823
(3)对两种递推公式的误差进行分析:
设第一递推式中开始时的误差为
,递推过程的舍入误差不计。
并记
,则有
。
由于
,因此递推式不可靠。
而在第二种递推式中
,误差在缩小,所以此递推式是可靠的。
2.求方程
的近似根,要求
,并比较计算量。
(1)在[0,1]上用二分法;
(2)取初值
,并用迭代
;
(3)加速迭代的结果;
(4)取初值
,并用牛顿迭代法;
(5)分析绝对误差。
解:
(1)用二分法,编写程序代码为下:
r=0.0005;
a=0;
b=1;
y1=exp(a)+10*a-2;
y2=exp(b)+10*b-2;
k=1;
whileabs(b-a)>r
k=k+1;
x=(a+b)/2;
y=exp(x)+10*x-2;
ify1*y<0
b=x
else
a=x
end
end
disp(x)
disp(k)
其运算结果为:
b=
0.5000
b=
0.2500
b=
0.1250
a=
0.0625
b=
0.0938
a=
0.0781
a=
0.0859
a=
0.0898
b=
0.0918
b=
0.0908
a=
0.0903
0.0903
12
(2)用题设迭代法,取初值x
(1)=0,并用迭代x(i)=(2-exp(x(i-1)))/10,程序代码如下:
x
(1)=0;
r=0.0005;
fori=2:
20
x(i)=(2-exp(x(i-1)))/10;
ifabs(x(i)-x(i-1)) return; end end disp(x(i)) disp(i) 其运算结果为: 0.0905 5 (3)加速迭代: 构造序列{ },使 = ,其中 是x的一个近似值,经计算: =0; =0.1; =0.0894829; =0.090639; =0.0905126; =0.090484; =0.090524449; =0.090525536。 编程如下: x=0;A=0;B=1;n=0 whileabs(B-A)>5*1e-4 A=x; y=exp(x)+10*x-2; z=exp(y)+10*y-2; x=x-(y-x)^2/(z-2*y+x); B=x; n=n+1 end x n 结果为: x=0.0995;n=3 (4)用牛顿法,取初值x=0,编写程序代码为下: r=0.00005; x=0; fori=1: 20 x=x-(exp(x)+10*x-2)/(exp(x)+10); ifabs((exp(x)+10*x-2)/(exp(x)+10)) return; end end disp(x) disp(i) 运算结果为: 0.0905 2 (5) 不同的计算方式是有收敛速度上的不同的。 由本题可以看到,在要求同一运算精度的情况下,采用二分法运算了12次,采用题设的迭代方法运算了5次,采用加速迭代法运算了3次,采用牛顿迭代法则只需运算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 试从中找出使用次数和容积之间的关系,计算均方差。 (注: 增速减少,用何种模型) 解: (1)由数据点分布图可知,拟合曲线y=f(x)随着x的增加而上升,但上升速度由快到慢,当x趋于无穷大时,y趋于某个常数,故曲线有一水平渐进线。 根据上述特征很容易想到用Logistic模型来拟合该曲线。 数据点分布图 设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 方程即为 (2)校核方差 x=[2345678910111213141516]; y=[6.428.209.589.509.7010.009.939.9910.4910.5910.6010.8010.6010.9010.76]; f(x)=11.6791*exp(-1.1107./x); plot(x,y,'*') holdon; plot(x,f(x),'-') 可以得到散点与拟合曲线图如下所示: sum=(11.6791*exp(-1.1107/x (1))-y (1))^2; fori=1: 14 sum=sum+(11.6791*exp(-1.1107/x(i+1))-y(i+1))^2; fc=sum^0.5; end fc 运行结果: fc= 0.9441 所以均方差是0.9441 4.设 , , 分析下列迭代法的收敛性,并求 的近似解及相应的迭代次数。 (1)JACOBI迭代; (2)GAUSS-SEIDEL迭代; (3)SOR迭代( )。 解: (1) A=[4,-1,0,-1,0,0;-1,4,-1,0,-1,0;0,-1,4,-1,0,-1;-1,0,-1,4,-1,0;0,-1,0,-1,4,-1;0,0,-1,0,-1,4]; b=[0;5;-2;5;-2;6]; esp=0.0001; x0=[0;0;0;0;0;0]; L=-tril(A,-1); U=-triu(A,1); D=diag(diag(A)); B=inv(D)*(L+U); f=inv(D)*b; J=B*x0+f; n=1; whilenorm(J-x0)>=esp x0=J; J=B*x0+f; n=n+1; end disp(J); disp(n) 运行的结果如下: 1.0000 2.0000 1.0000 2.0000 1.0000 2.0000 28 (2)用Gauss-Seidel迭代,程序如下: esp=0.0001; A=[4,-1,0,-1,0,0;-1,4,-1,0,-1,0;0,-1,4,-1,0,-1;-1,0,-1,4,-1,0;0,-1,0,-1,4,-1;0,0,-1,0,-1,4]; b=[0;5;-2;5;-2;6]; x0=[0;0;0;0;0;0] I=eye(6,6) L=-tril(A,-1); U=-triu(A,1); D=diag(diag(A)); L1=inv(D)*L U1=inv(D)*U B=inv(D)*(L+U); f=inv(D)*b; J=inv(I-L1)*U1*x0+inv(I-L1)*f; n=1; whilenorm(J-x0)>=esp x0=J; J=inv(I-L1)*U1*x0+inv(I-L1)*f; n=n+1; end ifn>400 disp('该方程不收敛') else disp(J); disp(n) end 运行结果如下: x0= 0 0 0 0 0 0 I= 100000 010000 001000 000100 000010 000001 L1= 000000 0.250000000 00.25000000 0.250000.2500000 00.250000.250000 000.250000.25000 U1= 00.250000.250000 000.250000.25000 0000.250000.2500 00000.25000 000000.2500 000000 1.0000 2.0000 1.0000 2.0000 1.0000 2.0000 15 (3)SOR迭代法(ω依次取1.334,1.95,0.95),程序如下: 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 以文件名sor.m保存。 程序: a=[4-10-100;-14-10-10;0-14-10-1;-10-14-10;0-10-14-1;00-10-14]; b=[05-25-26]'; x0=[000000]'; c=[1.3341.950.95]; fori=1: 3 w=c(i); sor(a,b,w,x0); end 运行结果分别为: y= 1.0000 2.0000 1.0000 2.0000 1.0000 2.0000 n= 13 y= 1.0000 2.0000 1.0000 2.0000 1.0000 2.0000 n= 241 y= 1.0000 2.0000 1.0000 2.0000 1.0000 2.0000 n= 17 可见w取1.334时,n最小,迭代次数最少。 上机心得: 从题中可以看出,系数矩阵严格对角占优,则Jacobi迭代法与Gauss-Seidel迭代法均收敛;又因系数矩阵对称正定,且0<w<2,则SOR迭代法也收敛。 三种迭代法的结果各有特色: (1)Jacobi法,其设计思想是将线性方程组的求解归结为对角方程组求解过程的重复,计算规则简单而易于编写汁算程序。 但,我们通常要设置最大迭代次数N,以防止迭代过程不收敛或者收敛速度过于缓慢。 (2)Gauss-Seidel充分利用新信息进行计算,其迭代的效果比Jacobi迭代好。 (3)SOR迭代法, 计算公式简单、编制程序容易,是求解大型稀疏方程组的一种有效方法,且超松弛法收敛速度比Gauss-Seidel迭代和Jacobi迭代都要快。 从对松弛因子三种不同选择的分析发现,若松弛因子选择得当,能显著地提高收敛速度。 5.用逆幂迭代法求 最接近于11的特征值和特征向量,准确到 。 解: A=[631;321;111]; [x,lamda]=eig(A), v=[111]'; q=11, B=inv(A-q*eye(3)); fork=1: 20 z=B*v; m=max(z); v=z/m; fprintf('%6d%f%f%f%f\n',k,v (1),v (2),v(3),m^-1+q), ifabs(m^-1+q-lamda(3,3))<0.001 break end end 运行结果: x= -0.3065-0.4082-0.8599 0.78120.4082-0.4722 -0.54380.8165-0.1938 lamda= 0.127000 01.00000 007.8730 q= 11 12.3571431.5714291.0000004.928571 23.5195172.0566911.0000006.576208 34.1081822.3010681.0000007.409250 44.3291432.3923811.0000007.721524 54.4024932.4225711.0000007.825064 64.4258342.4321421.0000007.857975 74.4331612.4351361.0000007.868297 84.4354522.4360701.0000007.871521 94.4361672.4363601.0000007.872527 6.用经典R-K方法求解初值问题 (1) , , ; (2) , , 。 和精确解 比较,分析结论。 解: (1)程序: functionydot=lorenzeq(x,y) ydot=[-2*y (1)+y (2)+2*sin(x);y (1)-2*y (2)+2*cos(x)-2*sin(x)] 以文件名lorenzeq.m保存。 主窗口输入: [x,y]=ode45('lorenzeq',[0: 10],[2;3]) 运行结果为: x= 0 1 2 3 4 5 6 7 8 9 10 y= 2.00003.0000 1.57751.2758 1.1802-0.1457 0.2406-0.8903 -0.7202-0.6170 -0.94540.2971 -0.27450.9652 0.65890.7557 0.9901-0.1449 0.4124-0.9109 -0.5440-0.8389 (2)程序: functionydot=lorenzeq1(x,y) ydot=[-2*y (1)+y (2)+2*sin(x);998*y (1)-999*y (2)+999*cos(x)-999*sin(x)]; 以文件名lorenzeq1.m保存。 程序: x=0: 10; y1=2*exp(-x)+sin(x); y2=2*exp(-x)+cos(x); [x,y]=ode45('lorenzeq1',[0: 10],[2;3]); fprintf('xy (1)y1y (2)y2\n') forj=1: length(y) fprintf('%4d%.4f%.4f%.4f%.4f\n',x(j),y(j,1),y1(j),y(j,2),y2(j)) end 运行结果为: xy (1)y1y (2)y2 02.00002.00003.00003.0000 11.57721.57721.27591.2761 21.18001.1800-0.1455-0.1455 30.24070.2407-0.8904-0.8904 4-0.7202-0.7202-0.6169-0.6170 5-0.9454-0.94540.29720.2971 6-0.2745-0.27450.96480.9651 70.65880.65880.75540.7557 80.99000.9900-0.1448-0.1448 90.41240.4124-0.9106-0.9109 10-0.5439-0.5439-0.8389-0.8390 结论: R-K方法求解的结果精度较高。 7.用有限差分法求解边值问题(h=0.1): 解: clear 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.95242759082864 0.90495042441637 0.85783523817387 0.81149210364572 0.76644735648340 0.72331872771229 0.68279404636095 0.64561505583677 0.61256800166995 0.58448274831666 0.56224232244653 0.54680502867801 0.53924172732245 0.54079161115864 0.55294101057354 0.57753158272634 0.61690696339694 0.67411095530971 0.75315614217452 0.85939026577182 8.P100,例4.4.6 function[err,a,b]=nlfit(x,y) ifnargin<2, x=[1: 8]'/10; y=[0.61.11.61.82.01.91.71.3]'; end c=fminsearch(@fitfun,[0;0],optimset,x,y); fprintf('Thenonlinearleastaquarefittingy=a*sin(b*x)fordata\n\n'); fprintf('%6.1f',x); fprintf('\n'); fprintf('%6.1f',y); fprintf('\n\nis\n\ty=%7.4f*sin(%7.4f*x)\n\n',c (1),c (2)); z=linspace(x (1),x(end),100);%求等差数列 plot(x,y,'r+',z,c (1)*sin(c (2)*z),'b-.') functionerr=fitfun(c,x,y) a=c (1); b=c (2); err=y-a*sin(b*x); err=err'*err; 运行结果: Thenonlinearleastaquarefittingy=a*sin(b*x)fordata 0.10.20.30.40.50.60.70.8 0.61.11.61.82.01.91.71.3 is y=1.9751*sin(3.0250*x) 9.p108.数值试验四.第1题. 次序年份人口(亿) a)19535.82 b)19646.95 c)198210.08 d)190011.34 e)200012.66 解: 把f(x)(1+cx)≈(a+bx)变成f(x)≈a+bx-cxf(x)则近似看成基函数是1,x,-x*f(x)而数据是(xi,f(xi))的最小二乘拟合问题,程序如下: x=[19531964198219002000]'; y=[5.826.9510.0811.3412.66]'; A=[ones(5,1)x-x.*y]; Z=A\y; a=Z (1) b=Z (2) c=Z(3) 结果: a= 11.5250 b= -0.0059 c= -5.0979e-004 10.p149.数值试验五.第4题.(讨论Logisitc模型) 年份 1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 人口 76.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.5 251.4 试分别用两点公式和三点公式计算美国人口20世纪的年增长率。 解
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 计算 上机 报告 上海 电力 学院