计算方法上机答案.docx
- 文档编号:18379474
- 上传时间:2023-08-16
- 格式:DOCX
- 页数:28
- 大小:180.02KB
计算方法上机答案.docx
《计算方法上机答案.docx》由会员分享,可在线阅读,更多相关《计算方法上机答案.docx(28页珍藏版)》请在冰点文库上搜索。
计算方法上机答案
上海电力学院
数值分析上机实验报告
题 目:
数值分析上机实验报告
学生姓名:
11111111111
学号:
111111111111111
专业:
1111
2013年12月30日
数值计算方法上机实习题
1.设
,
(1)由递推公式
,从
的几个近似值出发,计算
;
(2)粗糙估计
,用
,计算
;
(3)分析结果的可靠性及产生此现象的原因(重点分析原因)。
(1)解答:
n=0,
这里可以用for循环,while循环,根据个人喜好与习惯:
for循环程序:
While循环程序:
I=0.1823;I=0.1823;
forn=1:
20i=1;
I=(-5)*I+1/n;whilei<21
EndI=(-5)*I+1/i;
Ii=i+1;
fprintf('I20=%f',I)end
I=-2.0558e+009>>I
I20=-2055816073.851284>>I=-2.0558e+009
(2)
粗略估计I20:
Mathcad计算结果:
for循环程序:
While循环程序:
>>I=0.007998;I=0.007998;
>>forn=1:
20n=1;
I=(-0.2)*I+1/(5*n);whilen<21
EndI=(-0.2)*I+1/(5*n);
>>In=n+1;
I=0.0083end
>>I
I=0.0083
(3)算法误差分析:
计算在递推过程中传递截断误差和舍入误差
第一种算法:
(从1——>20)
误差放大了5n倍,算法稳定性很不好;
第二种算法:
(从20——>1)
误差在逐步缩小,算法趋近稳定,收敛。
2.求方程
的近似根,要求
,并比较计算量。
(1)在[0,1]上用二分法;
function[ti]=erfenfa(a,b)
f=@(x)(exp(x)+10*x-2)
t=(a+b)./2;
i=0;
whileabs(f(t))>0.001
iff(a)*f(t)<0
b=t;t=(a+b)/2;
elseiff(b)*f(t)<0
a=t;t=(b+a)/2;
end
i=i+1;
end
结果:
t=
0.0906
i=
11
(2)取初值
,并用迭代
;
functionx=diedai(x0)%x0初值
x=x0;
fori=1:
10000
y=(2-exp(x))./10;x=y;y=(2-exp(x))./10;
ifabs(x-y)<5*10^(-4)
disp('迭代次数');
2*i
break;
end
end
结果:
ans=
6
x=
.0906********
(3)加速迭代的结果(艾特肯Aitken加速方法);
function[ym]=aitken(a)
func=@(x)((2-exp(x))./10)
x
(1)=a;
wucha=1;m=1;
whilewucha>5*10^(-4)
p(m+1)=func(x(m));
q(m+1)=func(p(m+1));
x(m+1)=q(m+1)-((q(m+1)-p(m+1))^2)./(q(m+1)-2*p(m+1)+x(m));
wucha=abs(x(m+1)-x(m));
m=m+1;
ifm>1000
break;
end
end
formatlong
y=x(m-1);
m=m-1;
运行结果y=0.090483741803596
m=2
(4)取初值
,并用牛顿迭代法;
functionx=newtondiedai(x0)
x=x0;
fori=1:
100
y=x-(exp(x)+10*x-2)./(exp(x)+10);x=y;
y=x-(exp(x)+10*x-2)./(exp(x)+10);
ifabs(x-y)<0.0001
disp('迭代次数');
i
break;
end
end
运行结果
迭代次数
i=
2
x=
0.0905
(5)分析绝对误差。
通过指令求得方程精确解的近似解:
>>solve('exp(x)+10*x-2=0')
ans=
-lambertw(1/10*exp(1/5))+1/5
>>formatlong
>>-lambertw(1/10*exp(1/5))+1/5
ans=
0.09052510130725
小结:
所以我们可以看到,在要求同一运算精度的情况下,采用二分法运算了11次,采用题设的迭代方法运算了6次,采用加速迭代法只运算了2次,采用牛顿迭代法需运算2次。
也就是说加速迭代速度都是超线性收敛的,而简单迭代法是线性收敛的。
而二分法收敛速度较慢,所以在工程上不经常使用。
3.钢水包使用次数多以后,钢包的容积增大,数据如下:
x
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
y
6.42
8.2
9.58
9.5
9.7
10
9.93
9.99
10.49
10.59
10.60
10.8
10.6
10.9
10.76
试从中找出使用次数和容积之间的关系,计算均方差。
(注:
增速减少,用何种模型)
解答:
因为不知道其函数形式,可先plot而后确定使用哪种模型比较合适。
函数图形程序:
x=[2345678910111213141516];
y=[6.428.29.589.59.7109.939.9910.4910.5910.610.810.610.910.76];
plot(x,y,‘b*’)
与指数函数趋势类似(但是趋势不同,一个趋于递减,一个趋于递增,使用倒数),故拟合模型为:
两边同时取对数:
用此模型进行数据拟合:
x=[2345678910111213141516];
y=[6.428.29.589.59.7109.939.9910.4910.5910.610.810.610.910.76];
t=1./x;
s=polyfit(t,log(y),1)
s=
-1.11072.4578
即是:
最终拟合曲线为:
将此拟合曲线和源数据进行对比:
主程序:
x=[2345678910111213141516];
y=[6.428.29.589.59.7109.939.9910.4910.5910.610.810.610.910.76];
plot(x,y,'ro')
holdon
lims1=[2,16];
fplot('11.679*exp(-1.1107/x)',lims1)
holdoff
均方差的求解
x=[2:
16];
y=[6.428.29.589.59.7109.939.9910.4910.5910.6010.810.610.910.76];
f(x)=11.679*exp(-1.1107./x);
c=0;
fori=1:
15
a=y(i);
b=x(i);
c=c+(a-f(b))^2;
end
averge=c/15
averge=sqrt(averge)
0.2438
4.设
,
,
分析下列迭代法的收敛性,并求
的近似解及相应的迭代次数。
1)JACOBI迭代;
2)GAUSS-SEIDEL迭代;
3)SOR迭代(
)。
解答:
(1)JACOBI迭代;
定义函数。
functiontx=jacobi(A,b,imax,x0,tol)%jacobi迭代%初始值x0,次数imax,精度tol
del=10^-10;
tx=[x0];n=length(x0);
fori=1:
n
dg=A(i,i);
ifabs(dg) disp('对角元太小'); return end end fork=1: imax%jacobi循环 fori=1: n sm=b(i); forj=1: n ifj~=i sm=sm-A(i,j)*x0(j); end end x(i)=sm/A(i,i); end tx=[tx;x]; ifnorm(x-x0) return else x0=x; end end 保存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]; imax=100; tol=10^-4; tx=jacobi(A,b,imax,x0,tol); forj=1: size(tx,1) fprintf('%d%f%f%f%f%f%f\n',j,tx(j,1),tx(j,2),tx(j,3),tx(j,4),tx(j,5),tx(j,6)); end 运行结果: 10.0000000.0000000.0000000.0000000.0000000.000000 20.0000001.250000-0.5000001.250000-0.5000001.500000 30.6250001.0000000.5000001.0000000.5000001.250000 40.5000001.6562500.3125001.6562500.3125001.750000 50.8281251.5312500.7656251.5312500.7656251.656250 60.7656251.8398440.6796881.8398440.6796881.882813 70.9199221.7812500.8906251.7812500.8906251.839844 80.8906251.9252930.8505861.9252930.8505861.945313 90.9626461.8979490.9489751.8979490.9489751.925293 100.9489751.9651490.9302981.9651490.9302981.974487 110.9825741.9523930.9761961.9523930.9761961.965149 120.9761961.9837420.9674841.9837420.9674841.988098 130.9918711.9777910.9888951.9777910.9888951.983742 140.9888951.9924150.9848311.9924150.9848311.994448 150.9962081.9896390.9948201.9896390.9948201.992415 160.9948201.9964620.9929231.9964620.9929231.997410 170.9982311.9951670.9975831.9951670.9975831.996462 180.9975831.9983490.9966991.9983490.9966991.998792 190.9991751.9977450.9988731.9977450.9988731.998349 200.9988731.9992300.9984601.9992300.9984601.999436 210.9996151.9989480.9994741.9989480.9994741.999230 220.9994741.9996410.9992821.9996410.9992821.999737 230.9998201.9995090.9997551.9995090.9997551.999641 240.9997551.9998320.9996651.9998320.9996651.999877 250.9999161.9997710.9998861.9997710.9998861.999832 260.9998861.9999220.9998441.9999220.9998441.999943 270.9999611.9998930.9999471.9998930.9999471.999922 280.9999471.9999640.9999271.9999640.9999271.999973 290.9999821.9999500.9999751.9999500.9999751.999964 (2)GAUSS-SEIDEL迭代; functiontx=gseidel(A,b,imax,x0,tol)%gseidel迭代%初始值x0,次数imax,精度tol del=10^-10; tx=[x0];n=length(x0); fori=1: n dg=A(i,i); ifabs(dg) disp('对角元太小'); return end end fork=1: imax%G-S循环 x=x0; fori=1: n sm=b(i); forj=1: n ifj~=i sm=sm-A(i,j)*x(j); end end x(i)=sm/A(i,i); end tx=[tx;x]; ifnorm(x-x0) return else x0=x; end end 主程序: 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]; imax=100; tol=10^-4; tx=gseidel(A,b,imax,x0,tol); forj=1: size(tx,1) fprintf('%d%f%f%f%f%f%f\n',j,tx(j,1),tx(j,2),tx(j,3),tx(j,4),tx(j,5),tx(j,6)) end (3)SOR迭代( )。 functiontx=sor(A,b,imax,x0,tol,w)%sor迭代%初始值x0,次数imax,精度tol del=10^-10; tx=[x0];n=length(x0); fori=1: n dg=A(i,i); ifabs(dg) disp('对角元太小'); return end end fork=1: imax%G-S循环 x=x0; fori=1: n sm=b(i); forj=1: n ifj~=i sm=sm-A(i,j)*x(j); end end x(i)=sm/A(i,i); x(i)=w*x(i)+(1-w)*x0(i); end tx=[tx;x]; ifnorm(x-x0) return else x0=x; end end 主程序: 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]; imax=100; tol=10^-4; w=1.334; tx=sor(A,b,imax,x0,tol,w); forj=1: size(tx,1) fprintf('%d%f%f%f%f%f%f\n',j,tx(j,1),tx(j,2),tx(j,3),tx(j,4),tx(j,5),tx(j,6)) end 分析: 从题中可以看出,系数矩阵严格对角占优,则Jacobi迭代法与Gauss-Seidel迭代法均收敛;又因系数矩阵对称正定,且0<w<2,则SOR迭代法也收敛。 三种迭代法的结果分析及比较: (1)Jacobi法,其设计思想是将线性方程组的求解归结为对角方程组求解过程的重复,计算规则简单而易于编写汁算程序。 但,我们通常要设置最大迭代次数N,以防止迭代过程不收敛或者收敛速度过于缓慢。 (2)Gauss-Seidel充分利用新信息进行计算,其迭代的效果比Jacobi迭代好。 (3)SOR迭代法, 计算公式简单、编制程序容易,是求解大型稀疏方程组的一种有效方法,且超松弛法收敛速度比Gauss-Seidel迭代和Jacobi迭代都要快。 从对松弛因子三种不同选择的分析发现,若松弛因子选择得当,能显著地提高收敛速度。 5.用逆幂迭代法求 最接近于11的特征值和特征向量,准确到 。 解答: 调用函数: function[mt,my]=maxtr(A,p,ep) n=length(A); B=A-p*eye(n); v0=ones(n,1); k=1; v=B*v0; whileabs(norm(v,inf)-norm(v0,inf))>ep %norm(v-v0)>ep k=k+1; q=v; u=v/norm(v,inf) v=B*u; v0=q; end mt=1/norm(v,inf)+p my=u 主程序: A=[631;321;111]; maxtr(A,11,0.001) 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)]; 主窗口输入: x=0: 10; y1=2*exp(-x)+sin(x); y2=2*exp(-x)+cos(x); [x,y]=ode45('lorenzeq',[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.57751.57721.27581.2761 21.18021.1800-0.1457-0.1455 30.24060.2407-0.8903-0.8904 4-0.7202-0.7202-0.6170-0.6170 5-0.9454-0.94540.29710.2971 6-0.2745-0.27450.96520.9651 70.65890.65880.75570.7557 80.99010.9900-0.1449-0.1448 90.41240.4124-0.9109-0.9109 10-0.5440-0.5439-0.8389-0.8390 (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)]; 主窗口输入: 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 小结: 四阶RungeKutta方法得到的结果已很接近精确解,证明这种迭代方法精确度很好,是一种有效的算法。 7.用有限差分法求解边值问题(h=0.1): . 解答: 主程序: h=0.1; n=(1-(-1))/h+1; x (1)=-1;x(n-1)=1; y (1)=1;y(n-1)=1; fori=1: n-1 x(i)=x (1)+(i-1)*h; q(i)=(1+x(i)^2); B(i)=2/(h^2)+q(i); end fori=1: n-2 C(i)=-1/(h^2); end H=diag(B)+diag(C,1)+diag(C,-1); g (1)=0+1/(h^2); g(n-1)=0+1/(h^2); fori=2: n-2 g(i)=0; end y=H\g' 运行结果: y= 0.9027 0.8235 0.7592 0.7074 0.6661
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 计算方法 上机 答案