数值计算方法上机实习题NEW.docx
《数值计算方法上机实习题NEW.docx》由会员分享,可在线阅读,更多相关《数值计算方法上机实习题NEW.docx(27页珍藏版)》请在冰豆网上搜索。
![数值计算方法上机实习题NEW.docx](https://file1.bdocx.com/fileroot1/2023-2/22/e39af2a9-bc94-4da3-a675-c5bfd4a556c2/e39af2a9-bc94-4da3-a675-c5bfd4a556c21.gif)
数值计算方法上机实习题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(errorbreak;
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(errorbreak;
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