数值计算方法上机实习题NEW.docx
- 文档编号:5403263
- 上传时间:2023-05-08
- 格式:DOCX
- 页数:27
- 大小:155.67KB
数值计算方法上机实习题NEW.docx
《数值计算方法上机实习题NEW.docx》由会员分享,可在线阅读,更多相关《数值计算方法上机实习题NEW.docx(27页珍藏版)》请在冰点文库上搜索。
数值计算方法上机实习题NEW
数值计算方法上机实习题
1.设
,
(1)由递推公式
,从
出发,计算
;
程序如下:
functionI=myhs(I0,n)
ifn>=1
I=myhs(I0,n-1)*(-5)+1/n;
elseifn==0
I=I0;
end
命令行窗口输入:
I0=0.1822;
I1=myhs(I0,20);
I0=0.1823;
I2=myhs(I0,20);
运行结果:
当I0=0.1822时,I20=-1.1593e+10。
当I0=0.1823时,I20=-2.0558e+9。
(2)
,
用
,计算
;
程序如下:
functionI=myhs2(I20,n)
if(n<20)
I=myhs2(I20,n+1)*(-1)/5+1/(5*(n+1));
elseifn==20
I=I20;
end
命令行窗口输入:
I20=0;
I1=myhs2(I20,20);
I20=10000;
I2=myhs2(I20,20);
运行结果:
当I20=0时,I0=0.182********3955。
当I20=10000时,I0=0.182********8812。
(3)分析结果的可靠性及产生此现象的原因(重点分析原因)。
根据上式,
假设In的真值为I*,误差为en,即e=I*-In。
综合递推式,有en=-5*en-1,这意味着哪怕开始只有一点点误差,只要n足够大,按照这种计算一步误差增长5倍的方式,所得的结果总是不可信的,因此整个算法是不稳定的。
根据上式,假设In的真值为I*,误差为en,即e=I*-In。
综合递推式,有en-1.=(-1/5)*en,按照这种计算误差会以每步缩小到1/5的方式进行,根据
(2)得到的结果而言,该算法是相对稳定的。
2.求方程
的近似根,要求
,并比较计算量。
(1)在[0,1]上用二分法;
输入程序如下:
a=0;
b=1;
n=0;
while(abs(b-a)>5*1e-4)
c=(a+b)/2;
y=exp(c)+10*c-2;
if(y>0)
b=c;
else
a=c;
end
n=n+1;
end
c
得到结果:
c=0.090332,n=11
(2)取初值
,并用迭代
;
输入程序如下:
xk=0;
xki=(2-exp(xk))/10;
n=1;
while(abs(xki-xk)>5*1e-4)
xk=xki;
xki=(2-exp(xk))/10;
n=n+1;
end
xki
程序运行结果为:
xki=0.090513,n=4
(3)加速迭代的结果;
输入程序如下
xk=0;
yk=(2-exp(xk))/10;
zk=(2-exp(yk))/10;
xki=xk-(yk-xk)^2/(zk-2*yk+xk);
n=1;
while(abs(xki-xk)>5*1e-4)
xk=xki;
yk=(2-exp(xk))/10;
zk=(2-exp(yk))/10;
xki=xk-(yk-xk)^2/(zk-2*yk+xk);
n=n+1;
end
xki
程序运行的结果为
xki=0.090525,n=2
(4)取初值
,并用牛顿迭代法;
输入程序如下:
xk=0;
xki=xk-(exp(xk)+10*xk-2)/(exp(xk)+10);
n=1;
while(abs(xki-xk)>5*1e-4)
xk=xki;
xki=xk-(exp(xk)+10*xk-2)/(exp(xk)+10);
n=n+1;
end
xki
程序运行的结果为
xki=0.090525,n=2
(5)分析绝对误差。
通过程序x=solve('exp(x)+10*x-2=0');
解得x=0.090525
用二分法求得x的值为0.090332,迭代次数为11,绝对误差为0.000193。
用不动点迭代法求得x的值为0.090513,迭代次数为4,绝对误差为0.000012。
用加速迭代法和牛顿迭代法求得x的值均为0.090525,迭代次数均为2,绝对误差为0。
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=(ax+b)/(c+x),将原模型变为(c+x)y=(ax+b),采用非线性最小二乘法。
应选取参数a、b、c使得
达到极小值。
具体方法:
对Z分别求a、b、c的偏导,并使偏导为0,相应的方程组如下
对上述方程组进行整理,写成AX=B的形式,利用MATLAB求解X,对应的就是a、b、c三个参数的值。
具体程序如下:
%方程组求解和画拟合曲线
C=xlsread('test3.xlsx');
x=C(1,:
);
y=C(2,:
);
A=[sum(x.^2),sum(x),-sum(x.*y);sum(x),15,-sum(y);sum(x.*y),sum(y),-sum(y.^2)];
B=[sum(x.^2.*y);sum(x.*y);sum(x.*y.^2)];
Z=A\B;
a=Z
(1);
b=Z
(2);
c=Z(3);
X=linspace(1,50,50);
Y=zeros(1,50);
fori=1:
50
Y(i)=(a*X(i)+b)/(c+X(i));
end
plot(x,y,'*',X,Y,'-');
%计算方差
z1=0;
fori=2:
16
z1=z1+Y(i);
end
z1=z1/15;
sum=0;
fori=2:
16
sum=sum+(Y(i)-z1)^2;
end
fc=sqrt(sum/15);
运行结果为:
a=11.3400,b=-12.4953,c=-0.3403
拟合方程为y=(11.34x-12.4953)/(-0.3403+x),均方差为1.2285
4.设
,
,
分析下列迭代法的收敛性,并求
的近似解及相应的迭代次数。
(1)JACOBI迭代;
解:
输入程序如下
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];
D=diag(diag(A));
L=D-tril(A);
U=D-triu(A);
Bj=D\(L+U);
[V,landa]=eig(Bj);
Maxlanda=abs(landa(1,1));
fori=2:
size(landa)
if(abs(landa(i,i))>Maxlanda)
Maxlanda=abs(landa(i,i));
end
end
运行结果为:
Maxlanda=0.6830
即雅克比迭代矩阵半径=0.683<1,雅克比迭代收敛。
%自定义的雅克比迭代函数
function[x,iter]=jacobi(A,b,x0,tol)
D=diag(diag(A));
L=D-tril(A);
U=D-triu(A);
x=D\(b+L*x0+U*x0);
iter=1;
error=norm(x-x0);
while(error>0.0001)
fprintf('%d%f%f%f%f%f%f%f\n',iter,x
(1),x
(2),x(3),x(4),x(5),x(6),error);
if(error break; end x0=x; x=D\(b+L*x0+U*x0); error=norm(x-x0); iter=iter+1; end fprintf('%d%f%f%f%f%f%f%f\n',iter,x (1),x (2),x(3),x(4),x(5),x(6),error); 命令行窗口输入以下程序即可得到近似解: formatshort; 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]; tol=0.0001; x0=zeros(size(B)); jacobi(A,B,x0,tol); 运行结果为: 10.0000001.250000-0.5000001.250000-0.5000001.5000002.423840 20.6250001.0000000.5000001.0000000.5000001.2500001.605654 30.5000001.6562500.3125001.6562500.3125001.7500001.094196 40.8281251.5312500.7656251.5312500.7656251.6562500.747228 50.7656251.8398440.6796881.8398440.6796881.8828130.510360 60.9199221.7812500.8906251.7812500.8906251.8398440.348582 70.8906251.9252930.8505861.9252930.8505861.9453130.238086 80.9626461.8979490.9489751.8979490.9489751.9252930.162616 90.9489751.9651490.9302981.9651490.9302981.9744870.111069 100.9825741.9523930.9761961.9523930.9761961.9651490.075861 110.9761961.9837420.9674841.9837420.9674841.9880980.051814 120.9918711.9777910.9888951.9777910.9888951.9837420.035390 130.9888951.9924150.9848311.9924150.9848311.9944480.024172 140.9962081.9896390.9948201.9896390.9948201.9924150.016510 150.9948201.9964620.9929231.9964620.9929231.9974100.011276 160.9982311.9951670.9975831.9951670.9975831.9964620.007702 170.9975831.9983490.9966991.9983490.9966991.9987920.005260 180.9991751.9977450.9988731.9977450.9988731.9983490.003593 190.9988731.9992300.9984601.9992300.9984601.9994360.002454 200.9996151.9989480.9994741.9989480.9994741.9992300.001676 210.9994741.9996410.9992821.9996410.9992821.9997370.001145 220.9998201.9995090.9997551.9995090.9997551.9996410.000782 230.9997551.9998320.9996651.9998320.9996651.9998770.000534 240.9999161.9997710.9998861.9997710.9998861.9998320.000365 250.9998861.9999220.9998441.9999220.9998441.9999430.000249 260.9999611.9998930.9999471.9998930.9999471.9999220.000170 270.9999471.9999640.9999271.9999640.9999271.9999730.000116 280.9999821.9999500.9999751.9999500.9999751.9999640.000079 线性方程组的雅克比迭代近似解为: x=[0.9999821.9999500.9999751.9999500.9999751.999964] 误差为0.000079 (2)GAUSS-SEIDEL迭代; 解: 首先,求GAUSS-SEIDEL迭代的收敛性,输入如下程序: 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]; D=diag(diag(A)); L=D-tril(A); U=D-triu(A); Bg=(D+L)\U; [V,landa]=eig(Bg); Maxlanda=abs(landa(1,1)); fori=2: size(landa) if(abs(landa(i,i))>Maxlanda) Maxlanda=abs(landa(i,i)); end end 运行结果为 Maxlanda=0.4806 即GAUSS-SEIDEL迭代矩阵半径=0.4806<1,GAUSS-SEIDEL迭代收敛。 接下来,求GAUSS-SEIDEL迭代的近似解: %自定义的GAUSS-SEIDEL迭代函数 function[x,iter]=gauseidel(A,b,x0,tol) D=diag(diag(A)); L=D-tril(A); U=D-triu(A); x=(D+L)\(b+U*x0); iter=1; error=norm(x-x0); while(error>0.0001) fprintf('%d%f%f%f%f%f%f%f\n',iter,x (1),x (2),x(3),x(4),x(5),x(6),error); if(error break; end x0=x; x=(D+L)\(b+U*x0); error=norm(x-x0); iter=iter+1; end fprintf('%d%f%f%f%f%f%f%f\n',iter,x (1),x (2),x(3),x(4),x(5),x(6),error); 命令行窗口输入如下程序: formatshort; 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]; tol=0.0001; x0=zeros(size(B)); gauseidel(A,B,x0,tol); 运行结果为: 10.0000001.250000-0.8125001.453125-1.1757811.9970703.115282 20.6757810.5839840.2165530.732971-0.3299711.5283551.847408 30.3292391.139336-0.2195021.140073-0.6877641.7268160.975595 40.5698520.880720-0.0034580.936461-0.5225911.6315120.499274 50.4542951.004914-0.1092351.033087-0.6016221.6777140.240174 60.5095000.944911-0.0585270.986851-0.5635121.6555100.115337 70.4829400.973755-0.0828491.009099-0.5818361.6661710.055438 80.4957140.959900-0.0711580.998402-0.5730331.6610480.026645 90.4895760.966559-0.0767771.003542-0.5772631.6635100.012805 100.4925250.963359-0.0740771.001072-0.5752301.6623270.006154 110.4911080.964896-0.0753741.002259-0.5762071.6628950.002957 120.4917890.964157-0.0747511.001689-0.5757381.6626220.001421 130.4914620.964513-0.0750501.001963-0.5759631.6627530.000683 140.4916190.964342-0.0749061.001831-0.5758551.6626900.000328 150.4915430.964424-0.0749761.001894-0.5759071.6627210.000158 160.4915800.964384-0.0749421.001864-0.5758821.6627060.000076 线性方程组的雅克比迭代近似解为: x=[0.4915800.964384-0.0749421.001864-0.5758821.662706] 误差为0.000076 (3)SOR迭代(取 ,找到迭代步数最少的 )。 解: 首先看SOR迭代的收敛性,输入以下程序: formatshort; 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]; D=diag(diag(A)); L=D-tril(A); U=D-triu(A); w=0.1: 0.1: 1.9; n=length(w); l=1; fori=1: n Bs=(D+w(i)*L)\(D-w(i)*D+w(i)*U); [V,landa]=eig(Bs); Maxlanda=abs(landa(1,1)); forj=2: size(landa) if(abs(landa(j,j))>Maxlanda) Maxlanda=abs(landa(j,j)); end end if(Maxlanda<1) W(l)=w(i); l=l+1; end end fori=length(W) fprintf('%f\n',W(i)); end 运行结果为: W=[0.10000.20000.30000.40000.50000.60000.70000.80000.90001.00001.1000] 即SOR迭代的w*为W内的值时,该迭代方法是收敛的。 然后求SOR迭代次数最少时的w*,输入以下程序: %自定义的sor迭代函数 function[x,iter,error]=sor(A,b,x0,w,tol) D=diag(diag(A)); L=D-tril(A); U=D-triu(A); x=(D-w*L)\(D*x0+w*(b-D*x0+U*x0)); iter=1; error=norm(x-x0); while(error>tol) x0=x; x=(D-w*L)\(D*x0+w*(b-D*x0+U*x0)); error=norm(x-x0); iter=iter+1; end 在命令行输入程序: formatshort; 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]; tol=0.0001; x0=zeros(size(B)); w=[0.10.20.30.40.50.60.70.80.91.01.1]; [X,Iter,Error]=sor(A,B,x0,0.1,tol); fori=2: length(w) x0=zeros(size(B)); [x,iter,error]=sor(A,B,x0,w(i),tol); if(Iter>iter) X=x; Iter=iter; W=w(i); Error=error; end end 运行得到: Iter=12X=[1.00002.00001.00002.00001.00002.0000] Error=6.8688e-05W=1.1000 即SOR迭代收敛的条件下,当w*=1.1时,迭代次数最少为12,对应的近似解为x=[1.00002.00001.00002.00001.00002.0000],误差为6.8688e-05 6.用经典R-K方法求解初值问题 (1) , , ; (2) , , 。 和精确解 比较,进行误差分析得到结论,图形显示精确解和数值解。 解: 1)利用欧拉公式求解式 (1)的数值解,然后将得到的数值解与精确解通过图形显示,步长设置为0.1,输入如下程序: clear; a=0;b=10;N=100;y1(1,: )=[2,3]; h=(b-a)/N;x=a: h: b; forn=1: N y1(n+1,1)=y1(n,1)+h*(-2*y1(n,1)+y1(n,2
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 计算方法 上机 实习 NEW