北理工数值分析大作业.docx
《北理工数值分析大作业.docx》由会员分享，可在线阅读，更多相关《北理工数值分析大作业.docx（26页珍藏版）》请在冰点文库上搜索。
北理工数值分析大作业
数值分析上机作业
第1章
1.1计算积分
，n=9。
（要求计算结果具有6位有效数字）
程序：
n=1:
19;
I=zeros（1,19）;
I（19）=1/2*（（exp（-1）/20）+（1/20））;
I（18）=1/2*（（exp（-1）/19）+（1/19））;
fori=2:
10
I（19-i）=1/（20-i）*（1-I（20-i））;
end
formatlong
disp（I（1:
19））
结果截图及分析：
在MATLAB中运行以上代码，得到结果如下图所示：
当计算到数列的第10项时，所得的结果即为
n=9时的准确积分值。
取6位有效数字可得
.
1.2分别将区间[-10.10]分为100,200,400等份，利用mesh或surf命令画出二元函数
z=
的三维图形。
程序：
>>x=-10:
0.1:
10;
y=-10:
0.1:
10;
[X,Y]=meshgrid（x,y）;
Z=exp（-abs（X））+cos（X+Y）+1./（X.^2+Y.^2+1）;
subplot（2,2,1）;
mesh（X,Y,Z）;
title（'步长0.1'）
>>x=-10:
0.2:
10;
y=-10:
0.2:
10;
[X,Y]=meshgrid（x,y）;
Z=exp（-abs（X））+cos（X+Y）+1./（X.^2+Y.^2+1）;
subplot（2,2,1）;
mesh（X,Y,Z）;
title（'步长0.2'）
>>x=-10:
0.05:
10;
y=-10:
0.05:
10;
[X,Y]=meshgrid（x,y）;
Z=exp（-abs（X））+cos（X+Y）+1./（X.^2+Y.^2+1）;
subplot（2,2,1）;
mesh（X,Y,Z）;
title（'步长0.05'）
结果截图及分析：
由图可知，步长越小时，绘得的图形越精确。
第2章
试用MATLAB编程实现追赶法求三对角方程组的算法，并考虑梯形电路电阻问题：
电路中的电流
满足下列线性方程组：
设
，求各段电路的电流量。
处理思路：
观察该方程的系数矩阵可知，它是一个三对角矩阵，故可运用追赶法对其进行求解。
程序:
fori=1:
8
a（i）=-2;b（i）=5;c（i）=-2;d（i）=0;
end
a
（1）=0;b
（1）=2;c（8）=0;d
（1）=220/27;
fori=2:
8
a（i）=a（i）/b（i-1）;
b（i）=b（i）-c（i-1）*a（i）;
d（i）=d（i）-a（i）*d（i-1）;
end
d（8）=d（8）/b（8）;
fori=7:
-1:
1
d（i）=（d（i）-c（i）*d（i+1））/b（i）;
end
fori=1:
8
x（i）=d（i）;
end
x
结果截图及分析：
在MATLAB中运行以上代码，得到结果如下图所示：
图中8个值依次为
的数值。
第3章
试分别用
（1）Jacobi迭代法；
（2）Gauss-Seidel解线性方程组
迭代初始向量取
.
3.1Jacobi迭代法
程序：
>>A=[101234;
19-12-3;
2-173-5;
32312-1;
4-3-5-115];
b=[12;-27;14;-17;12];
x0=[0;0;0;0;0];
D=diag（diag（A））;
I=eye（5）;
L=-tril（A,-1）;
B=I-D\A;
g=D\b;
y=B*x0+g;
n=1;
whilenorm（y-x0）>=1.0e-6
x0=y;
y=B*x0+g;
n=n+1;
end
fprintf（'%8.6f\n',y）;
n
结果截图及分析：
得到此结果时迭代次数为67次，达到精度要求。
3.2Gauss-Seidel迭代法：
程序：
>>A=[101234;
19-12-3;
2-173-5;
32312-1;
4-3-5-115];
b=[12;-27;14;-17;12];
x0=[0;0;0;0;0];
D=diag（diag（A））;
U=-triu（A,1）;
L=-tril（A,-1）;
M=（D-L）\U;
g=（D-L）\b;
y=M*x0+g;
n=1;
whilenorm（y-x0）>=1.0e-6
x0=y;
y=M*x0+g;
n=n+1;
end
fprintf（'%8.6f\n',y）;
结果截图及分析：
Gauss-Seidel迭代法只需要迭代38次即可满足精度要求。
第4章
设A=
，取
先用幂法迭代3次，得到A的按模最大特征值的近似值，取
为其整数部分，再用反幂法计算A的按模最大特征值的更精确的近似值，要求误差小于
.
程序：
A=[126-6;
6162;
-6216];
x0=[1;1;1];
y=x0;b=max（abs（x0））;k=1;
while（k<4）
x=A*y;b=max（abs（x））;y=x./b;
k=k+1;
fprintf（'eig1equals%6.4f\n',b）;
end
>>bb0=fix（b）;
I=eye（3,3）;
x0=[1;1;1];
y=x0;l=0;bb=max（abs（x0））;k=1;
while（abs（bb-l）>=1.0e-10）
l=bb;
x=（A-bb0*I）\y;bb=max（abs（x））;y=x./bb;
eig=l+b;
>>fprintf（'eig2（%d）equals%12.10f\n',k,eig）;
k=k+1;
end
实验截图及分析：
由图可知，由幂法3次迭代后得到的特征值为19.4，而由反幂法得到的特征值为20.3999999999.误差小于
第5章
试编写MATLAB函数实现Newton插值，要求能输出插值多项式。
对函数f（x）=
在区间[-5,5]上实现10次多项式插值。
要求：
（1）输出插值多项式。
（2）在区间[-5,5]内均匀插入99个节点，计算这些节点上函数f（x）的近似值，并在同一张图上画出原函数和插值多项式的图形。
（3）观察龙格现象，计算插值函数在各节点处的误差，并画出误差图。
5.1输出插值多项式
程序：
x=-5:
1:
5;
y=1./（1+4*（x.^2））;
newpoly（x,y）
function[c,d]=newpoly（x,y）
n=length（x）;
d=zeros（n,n）;
d（:
1）=y';
forj=2:
n
fork=j:
n
d（k,j）=（d（k,j-1）-d（k-1,j-1））/（x（k）-x（k-j+1））;
end
end
c=d（n,n）;
fork=（n-1）:
-1:
1
c=conv（c,poly（x（k）））;
m=length（c）;
c（m）=c（m）+d（k,k）;
end
end
结果及分析：
ans=
Columns1through2
-0.000049595763049
Columns3through4
0.0027401659084830.000000000000000
Columns5through6
.0514********
Columns7through8
0.3920149852823120.000000000000000
Columns9through10
-1.1432840483510250.000000000000001
Column11
1.000000000000000
10次插值多项式由高到低系数为Columns1至Column11
5.2原函数与插值多项式的图形
程序：
x=-5:
1:
5;
y=1./（1+4*（x.^2））;
n=newpoly（x,y）;
x0=-5:
0.1:
5;
y0=1./（1+4*（x0.^2））;
vn=polyval（n,x0）;
plot（x0,vn,'-r',x0,y0,'--b'）;
xlabel（'x'）;ylabel（'y'）;
实验结果截图：
原函数与插值多项式的图形如上图所示，蓝色为原函数的图形，红色为插值多项式的图形。
5.3各节点的误差及误差图
程序：
formatlong;
x=-5:
1:
5;
y=1./（1+4*（x.^2））;
n=newpoly（x,y）;
x0=-5:
0.1:
5;
y0=1./（1+4*（x0.^2））;
vn=polyval（n,x0）;
plot（x0,y0-vn,'-r'）;
xlabel（'x'）;ylabel（'y'）;
实验结果截图：
误差图如上图所示。
第6章
炼钢厂出钢时所用的盛钢水的钢包，在使用过程中由于钢液及炉渣对包衬耐火材料的腐蚀，使其容积不断加大。
经试验，钢包的容积与相应的使用次数的数据列表如下：
使用次数x
容积y
使用次数x
容积y
2
106.42
11
110.59
3
108.26
12
110.60
5
109.58
14
110.72
6
109.50
16
110.90
7
109.86
17
110.76
9
110.00
19
111.10
10
109.93
20
111.30
选用双曲线
对数据进行拟合，使用最小二乘法求出拟合函数，作出拟合曲线图。
处理思路：
用Y替代1/y，用X替代1/x，原曲线化为Y=a+bx，双曲线转化为一次线性方程，使用最小二乘法求出该一次方程的系数。
程序：
x=[2356791011121416171920];
y=[106.42108.26109.58109.5109.86110109.93110.59110.60110.72110.9110.76111.1111.3];
k1=0;k2=0;k3=0;k4=0;
fori=1:
14
k1=k1+1/x（i）;
end
fori=1:
14
k2=k2+1/y（i）;
end
fori=1:
14
k3=k3+1/（x（i））^2;
end
fori=1:
14
k4=k4+1/（x（i）*y（i））;
end
b=（k1*k2-14*k4）/（k1^2-14*k3）
a=k2/14-k1*b/14
plot（x,y,'r*'）
holdon
x=2:
0.01:
20;
y=1./（a+b./x）;
plot（x,y）
xlabel（'x'）
ylabel（'y'）
gridon
实验结果截图与分析：
即最小二乘法求出拟合函数为：
=0.008973+0.000842
拟合曲线图为：
第7章
考纽螺线的形状象钟表的发条，也称回旋曲线，它在直角坐标系中的参数方程为
曲线关于原点对称，取a=1，参数s的变化范围[-5,5]，容许误差限分别是
和
。
选取适当的节点个数，利用数值积分方法计算曲线上点的坐标，并画出曲线的图形。
误差限为
时：
程序：
x=zeros（101,1）;
y=zeros（101,1）;
f1=inline（'cos（1/2*（t.^2））'）;
f2=inline（'sin（1/2*（t.^2））'）;
i=1;
fors=-5:
0.1:
5
x（i,1）=quad（f1,0,s,1e-6）;
y（i,1）=quad（f2,0,s,1e-6）;
i=i+1;
end
plot（x,y,'r-'）;
title（'误差限-1e-6'）;
xlabel（'x（s）'）;
ylabel（'y（s）'）;
实验截图：
误差限为
时得到曲线如下图所示：
误差限为
时：
程序：
x=zeros（101,1）;
y=zeros（101,1）;
f1=inline（'cos（1/2*（t.^2））'）;
f2=inline（'sin（1/2*（t.^2））'）;
i=1;
fors=-5:
0.1:
5
x（i,1）=quad（f1,0,s,1e-10）;
y（i,1）=quad（f2,0,s,1e-10）;
i=i+1;
end
plot（x,y,'r-'）;
title（'误差限-1e-10'）;
xlabel（'x（s）'）;
ylabel（'y（s）'）;
实验结果截图：
第8章
求方程
的非零根。
程序：
x=sym（'x'）;
f=sym（'log（（513+0.6651*x）/（513-0.6651*x））-x/（1400*0.0981）'）;
df=diff（f,x）;
FX=x-f/df;
Fx=inline（FX）;
formatlong;
i=1;
x0=760;
fori=1:
10
disp（x0）;
x0=feval（Fx,x0）;
end
x0
结果截图及分析：
取初值为760，迭代9次，两个非零根为765.48及-765.48.
第9章
设有常微分方程的初值问题
其精确解为
。
选取步长h使四阶Adams预测-校正算法和经典RK法均稳定，分别用这两种方法求解微分方程，将数值解与精确解进行比较，输出结果。
其中多步法需要的初值由经典RK法提供。
程序：
f=inline（'-y+2*cos（x）','x','y'）;
n=20;
Y=zeros（1,n+1）;
Y
（1）=1;
h=0.05*pi;
X=0:
h:
pi;
fori=1:
20
K1=f（X（i）,Y（i））;
K2=f（X（i）+h/2,Y（i）+（h*K1）/2）;
K3=f（X（i）+h/2,Y（i）+（h*K2）/2）;
K4=f（X（i）+h,Y（i）+h*K3）;
Y（i+1）=Y（i）+h*（K1+2*K2+2*K3+K4）/6;
end;
Y1=cos（X）+sin（X）;
disp（'步长经典RK法精确解'）;
disp（[X'Y'Y1']）;
plot（X,Y1,'k*',X,Y,'-r'）;
gridon;
title（'经典RK法解非线性方程'）;
legend（'精确解','经典RK法'）;
实验结果截图：
四阶Adams预测-校正算法
程序：
f=inline（'-y+2*cos（x）','x','y'）;
n=20;
Y=zeros（1,n+1）;
Y
（1）=1;
h=pi/n;
X=0:
h:
pi;
fori=1:
3
K1=f（X（i）,Y（i））;
K2=f（X（i）+h/2,Y（i）+（h*K1）/2）;
K3=f（X（i）+h/2,Y（i）+（h*K2）/2）;
K4=f（X（i）+h,Y（i）+h*K3）;
Y（i+1）=Y（i）+h*（K1+2*K2+2*K3+K4）/6;
end;
Y1=Y;
fori=4:
n
f1=f（X（i-3）,Y（i-3））;
f2=f（X（i-2）,Y（i-2））;
f3=f（X（i-1）,Y（i-1））;
f4=f（X（i）,Y（i））;
Y（i+1）=Y（i）+h*（55*f4-59*f3+37*f2-9*f1）/24;
Yn=f（X（i+1）,Y（i+1））;
Y1（i+1）=Y（i）+h*（9*Yn+19*f4-5*f3+f2）/24;
end
Y2=cos（X）+sin（X）;
disp（'步长四阶Adams预测-校正算法精确解'）;
disp（[X'Y1'Y2']）;
plot（X,Y2,'k*',X,Y1,'-r'）;
gridon;
title（'四阶Adams预测-校正算法解微分方程'）;
legend（'精确解','四阶Adams预测-校正算法'）;
实验结果截图：
升级会员 