数值分析上机作业
第 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));
for i=2:10
I(19-i)=1/(20-i)*(1-I(20-i));
end
format long
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 编程实现追赶法求三对角方程组的算法，并考虑梯形电
路电阻问题：电路中的电流128{,,,}i i i L 满足下列线性方程组：
12123
234
345
456
567
6787822/25202520252025202520
2520
250
i i V R i i i i i i i i i i i i i i i i i i i i -=-+-=-+-=-+-=-+-=-+-=-+-=-+=
设220,27V V R ==Ω，求各段电路的电流量。
处理思路：观察该方程的系数矩阵可知，它是一个三对角矩阵，故可运用追赶法对其进行求解。 程序:
for i=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; for i=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); for i=7:-1:1
d(i)=( d(i)-c(i)*d(i+1) )/b(i); end for i=1:8 x(i)=d(i); end
x
结果截图及分析：在MATLAB 中运行以上代码，得到结果如下图所示：图中8
个值依次为128{,,,}i i i L 的数值。
第 3 章
试分别用（1）Jacobi 迭代法；（2）Gauss-Seidel 解线性方程组
123451012
34121912327217351432312117435
11512x x x x x ??????
??????---???
???
??????--=???
???--??????
?????
?---??????
迭代初始向量取(0)(0,0,0,0,0)T x =. 3.1 Jacobi 迭代法 程序：
>> A=[10 1 2 3 4;
1 9 -1
2 -3; 2 -1 7
3 -5; 3 2 3 12 -1;
4 -3 -
5 -1 15]; 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;
while norm(y-x0)>=1.0e-6 x0=y;
y=B*x0+g; n=n+1; end
fprintf('%8.6f\n',y); n
结果截图及分析：
得到此结果时迭代次数为67次，达到精度要求。
3.2 Gauss-Seidel迭代法：
程序：
>> A=[10 1 2 3 4;
1 9 -1
2 -3;
2 -1 7
3 -5;
3 2 3 12 -1;
4 -3 -
5 -1 15];
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;
while norm(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=???
?
?
?????--162621666612，取先用幂法迭代3次，得到
A 的按模最大特征值的近似值，取为其整数部分，再用反幂
法计算A 的按模最大特征值的更精确的近似值，要求误差小于
.
程序：
A=[12 6 -6;
6 16 2; -6 2 16]; 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('eig1 equals %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';
for j=2:n
for k=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);
for k=(n-1):-1:1
c=conv(c,poly(x(k)));
m=length(c);
c(m)=c(m)+d(k,k);
end
end
结果及分析：
ans =
Columns 1 through 2
-0.3049
Columns 3 through 4
0.8483 0.0000
Columns 5 through 6
-0.6720 0.0000
Columns 7 through 8
0.2312 0.0000
Columns 9 through 10
-1.1025 0.0001
Column 11
1.0000
10次插值多项式由高到低系数为Columns 1至Column 11
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');
实验结果截图：
y
x
原函数与插值多项式的图形如上图所示，蓝色为原函数的图形，红色为插值多项式的图形。
5.3各节点的误差及误差图
程序：
format long;
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');
实验结果截图：
y
x 误差图如上图所示。
第 6 章
炼钢厂出钢时所用的盛钢水的钢包，在使用过程中由于钢液及炉渣对包衬耐火材料的腐蚀，使其容积不断加大。经试验，钢包的容积与相应的使用次数的数据列表如下：
选用双曲线x
b a y
11+=对数据进行拟合，使用最小二乘法求出拟合函数，作出拟合曲线图。
处理思路：用Y 替代1/y ，用X 替代1/x ，原曲线化为Y=a+bx ，双曲线转化为一次线性方程，使用最小二乘法求出该一次方程的系数。
程序：
x=[2 3 5 6 7 9 10 11 12 14 16 17 19 20];
y=[106.42 108.26 109.58 109.5 109.86 110 109.93 110.59 110.60 110.72 110.9 110.76 111.1 111.3];
k1=0; k2=0; k3=0; k4=0;
for i=1:14
k1=k1+1/x(i);
end
for i=1:14
k2=k2+1/y(i);
end
for i=1:14
k3=k3+1/(x(i))^2;
end
for i=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*')
hold on
x=2:0.01:20;
y=1./(a+b./x);
plot(x,y)
xlabel('x')
ylabel('y')
grid on
实验结果截图与分析：
即最小二乘法求出拟合函数为：
=0.008973+0.000842
拟合曲线图为：
第 7 章
考纽螺线的形状象钟表的发条，也称回旋曲线，它在直角坐标系中的
参数方程为???
???
?
==
??s
s
dt at s y dt at s x 0
2
2
21sin )(21cos
)( 曲线关于原点对称，取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;
for s= -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;
for s= -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)');
实验结果截图：