数值计算方法上机实习题答案.docx
《数值计算方法上机实习题答案.docx》由会员分享,可在线阅读,更多相关《数值计算方法上机实习题答案.docx(14页珍藏版)》请在冰豆网上搜索。
![数值计算方法上机实习题答案.docx](https://file1.bdocx.com/fileroot1/2022-12/13/85969140-1860-4f70-ad64-01f6c701bc8f/85969140-1860-4f70-ad64-01f6c701bc8f1.gif)
数值计算方法上机实习题答案
1.设
,
(1)由递推公式
,从
的几个近似值出发,计算
;
解:
易得:
ln6-ln5=0.1823,
程序为:
I=0.182;
forn=1:
20
I=(-5)*I+1/n;
end
I
输出结果为:
=-3.0666e+010
(2)粗糙估计
,用
,计算
;
因为
所以取
程序为:
I=0.0087;
forn=1:
20
I=(-1/5)*I+1/(5*n);
end
I
(3)分析结果的可靠性及产生此现象的原因(重点分析原因)。
首先分析两种递推式的误差;设第一递推式中开场时的误差为
,递推过程的舍入误差不计。
并记
,那么有
。
因为
,所此递推式不可靠。
而在第二种递推式中
,误差在缩小,所以此递推式是可靠的。
出现以上运行结果的主要原因是在构造递推式过程中,考虑误差是否得到控制,即算法是否数值稳定。
2.求方程
的近似根,要求
,并比拟计算量。
(1)在[0,1]上用二分法;
程序:
a=0;b=1.0;
whileabs(b-a)>5*1e-4
c=(b+a)/2;
ifexp(c)+10*c-2>0
b=c;
elsea=c;
end
end
c
结果:
c=
(2)取初值
,并用迭代
;
程序:
x=0;
a=1;
whileabs(x-a)>5*1e-4
a=x;
x=(2-exp(x))/10;
end
x
结果:
x=
(3)加速迭代的结果;
程序:
x=0;
a=0;b=1;
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;
end
x
结果:
x=
(4)取初值
,并用牛顿迭代法;
程序:
x=0;
a=0;b=1;
whileabs(b-a)>5*1e-4
a=x;
x=x-(exp(x)+10*x-2)/(exp(x)+10);
b=x;
end
x
结果:
x=
(5)分析绝对误差。
solve('exp(x)+10*x-2=0')
3.钢水包使用次数多以后,钢包的容积增大,数据如下:
x
2
3
4
5
6
7
8
9
y
10
10
11
12
13
14
15
16
试从中找出使用次数和容积之间的关系,计算均方差。
〔注:
增速减少,用何种模型〕
设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
z=lny
10
11
12
13
14
15
16
程序:
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=
由此可得A=2.4578,B=-1.1107,
方程即为
计算均方差编程:
x=[2:
16];y=[6.428.29.589.59.7109.939.9910.4910.5910.6010.810.610.910.76];
f(x)=11.6791*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=
4.设
,
,
分析以下迭代法的收敛性,并求
的近似解及相应的迭代次数。
(1)JACOBI迭代;
程序:
functiony=jacobi(a,b,x0)
D=diag(diag(a));
U=-triu(a,1);
L=-tril(a,-1);
B=D\(L+U);
f=D\b;
y=B*x0+f;n=1;
whilenorm(y-x0)>1e-4
x0=y;
y=B*x0+f;n=n+1;
end
y
n
以文件名jacobi.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]';
jacobi(a,b,x0);
运行结果为:
y=
n=
28
(2)GAUSS-SEIDEL迭代;
程序:
functiony=seidel(a,b,x0)
D=diag(diag(a));
U=-triu(a,1);
L=-tril(a,-1);
G=(D-L)\U;
f=(D-L)\b;
y=G*x0+f;n=1;
whilenorm(y-x0)>10^(-4)
x0=y;
y=G*x0+f;n=n+1;
end
y
n
以文件名deisel.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]';
jacobi(a,b,x0);
运行结果为:
y=
n=
15
(3)SOR迭代〔
〕。
程序:
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=
n=
13
y=
n=
241
y=
n=
17
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=[1-2-3];maxtr(A,11,0.001)
结果为:
特征值:
mt=
特征向量:
my=
6.用经典R-K方法求解初值问题
〔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〕
,
,
。
和准确解
比拟,分析结论。
程序:
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
结论:
R-K方法求解的结果精度较高。
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=
8.拟合形如f〔x〕≈〔a+bx〕/(1+cx)的函数的一种快速方法是将最小二乘法用于以下问题:
f〔x〕(1+cx)≈〔a+bx〕,试用这一方法拟合表4-4给出的中国人口数据。
表4-4
次序年份人口〔亿〕
a)
b)
c)
d)
e)
解:
把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=
b=
c=