矩阵与数值分析课程数值实验大作业.docx
《矩阵与数值分析课程数值实验大作业.docx》由会员分享，可在线阅读，更多相关《矩阵与数值分析课程数值实验大作业.docx（20页珍藏版）》请在冰点文库上搜索。
矩阵与数值分析课程数值实验大作业
2011级工科硕士研究生
《矩阵与数值分析》课程数值实验
班级：
学号：
姓名：
任课教师：
大连理工大学
2011年12月20日
一、对于数列
，有如下两种生成方式
1、首项为
，递推公式为
；
2、前两项为
，递推公式为
；
给出利用上述两种递推公式生成的序列的第50项。
【按第一种递推公式】
clear
clc
a=1;
fori=1:
50-1
a=[aa（i）/3];
end
disp（'数列第50项小数表达为：
'）
formatlong
disp（a（50））
disp（'分数表达为：
'）
formatrat
disp（a（50））
formatshort
[运行结果]
数列第50项小数表达为：
4.178866707295615e-024
分数表达为：
1/2392993292300
【按第二种递推公式】
clear
clc
a=[11/3];
fori=2:
50-1
a=[a10/3*a（i）-a（i-1）];
end
formatrat
disp（'数列第50项为：
'）
disp（a（50））
formatshort
[运行结果]
数列第50项为：
2060436
【分析】
第一种算法数值稳定，计算过程舍入误差被严格控制，且按1/3的公差不断缩小。
但第二种算法数值不稳定。
另外，在第二种算法中，若将递推公式“a=[a10/3*a（i）-a（i-1）]”中的分母移动位置，改写成“a=[a10*a（i）/3-a（i-1）]”，则程序运行结果为-4966040，可以舍入误差被放大的十分严重。
二、利用迭代格式
及Aitken加速后的新迭代格式求方程
在
内的根
【未经加速的代码】
clc
eps=1e-15;
i=1;
x0=1;
formatlong
whilei<100
x1=sqrt（10/（x0+4））;
ifabs（x1-x0）<=eps
break
end
x0=x1;
i=i+1;
end
disp（'方程的解[精度10^（-15）]'）
disp（x1）
disp（'未经加速的迭代次数'）
disp（i）
[运行结果]
方程的解[精度10^（-15）]
1.36523001341410
未经加速的迭代次数
18
【经Aitken加速的代码】
clc
eps=1e-15;
i=1;
x0=1;
formatlong
whilei<100
x1=sqrt（10/（x0+4））;
y=sqrt（10/（x1+4））;
z=sqrt（10/（y+4））;
x1=z-（z-y）^2/（z-2*y+x1）;
ifabs（x1-x0）<=eps
break
end
x0=x1;
i=i+1;
end
disp（'方程的解[精度10^（-15）]'）
disp（x1）
disp（'未经加速的迭代次数'）
disp（i）
[运行结果]
方程的解[精度10^（-15）]
1.36523001341410
未经加速的迭代次数
3
【分析】
Aitken加速能对数列｛xk｝起明显的加速作用，在要求相同方程解精度的条件下，它能将迭代次数显著降低。
实际上，Aitken有时甚至能将发散的数列加速后变为收敛。
三、解线性方程组
1．分别Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组
迭代法计算停止的条件为：
．
2.用Gauss列主元消去法、QR方法求解如下方程组：
【1.Jacobi方法】
clc
i=1;
eps=1e-6;
A=[621-2;
250-2;
-2085;
1327];
b=[47-10]';
x0=zeros（4,1）;
D=diag（diag（A））;
L=-tril（A,-1）;
U=-triu（A,1）;
B=inv（D）*（L+U）;
f=inv（D）*b;
whilei<100
x1=B*x0+f;
ifnorm（x1-x0）<=eps
break
end
x0=x1;
i=i+1;
end
disp（'方程的解[精度10^（-6）]'）
disp（x1）
disp（'迭代次数'）
disp（i）
[运行结果]
方程的解[精度10^（-6）]
0.229
1.528
0.24463239101199
-0.572
迭代次数
16
【1.Gauss-Seidel方法】
clc
i=1;
eps=1e-6;
A=[621-2;
250-2;
-2085;
1327];
b=[47-10]';
x0=zeros（4,1）;
D=diag（diag（A））;
L=-tril（A,-1）;
U=-triu（A,1）;
B=inv（D-L）*（U）;
f=inv（D-L）*b;
whilei<100
x1=B*x0+f;
ifnorm（x1-x0）<=eps
break
end
x0=x1;
i=i+1;
end
disp（'方程的解[精度10^（-6）]'）
disp（x1）
disp（'迭代次数'）
disp（i）
[运行结果]
方程的解[精度10^（-6）]
0.111
1.986
0.24463241176981
-0.579
迭代次数
8
【2.Gauss列主元消去法】
clc
A=[2212;
413-1;
-4-201;
2323];
b=[1210]';
Ab=[Ab];
[N,N]=size（A）;
x=zeros（N,1）;
forp=1:
N-1
[max1,j]=max（abs（A（p:
N,p）））;
temp=Ab（p,:
）;
Ab（p,:
）=Ab（j+p-1,:
）;
Ab（j+p-1,:
）=temp;
ifAb（p,p）==0
disp（'方程无解'）
break
end
fork=p+1:
N
mult=Ab（k,p）/Ab（p,p）;
Ab（k,p:
N+1）=Ab（k,p:
N+1）-mult*Ab（p,p:
N+1）;
end
end
b=Ab（:
N+1）;
xx=zeros（N,1）;
fork=N:
-1:
1
xx（k）=b（k）/Ab（k,k）;
i=（1:
k-1）';
b（i）=b（i）-xx（k）*Ab（i,k）;
end
x=xx
[运行结果]
x=
1.54166666666667
-2.750
0.333
1.66666666666667
【2.QR分解法】
clc
fori=1:
3
A{i}=zeros（5-i）;
Q{i}=eye（4）;
end
x=zeros（4,1）;
y=zeros（4,1）;
R=zeros（4）;
A{1}=[2212;
413-1;
-4-201;
2323;];
b=[1,2,1,0]';
fori=1:
3
I=eye（size（A{i}））;
a=A{i}（:
1）;
w=a-norm（a）*I（:
1）;
hw=I-2/（w'*w）*（w*w'）;
Q{i}（i:
4,i:
4）=hw;
ifi==3
break
end
c=hw*A{i};
A{i+1}=c（2:
max（size（A{i}））,2:
max（size（A{i}）））;
end
Qz=（Q{3}*Q{2}*Q{1}）';
R=Q{3}*Q{2}*Q{1}*A{1};
c=Qz'*b;
x（4）=c（4）/R（4,4）;
fori=3:
-1:
1
x（i）=（c（i）-R（i,i+1:
4）*x（i+1:
4））/R（i,i）;
end
x
[运行结果]
x=
1.54166666666666
-2.750
0.333
1.66666666666666
【分析】
Jacobi迭代法和Gauss-Seidel迭代法都可用来求解线性方程组。
在同等精度下，求解本道题Jacobi迭代法迭代了16次，而Gauss-Seidel仅为8次，是前者的一半。
所以可以认为，在某些情况下，Gauss-Seidel是比较好的解法，但不总如此。
Gauss消去法可能发生小主元做除数，从而导致计算结果严重偏离真实值，所以在消元过程中，每一步都按列选主元，也就是Gauss列主元消去法，它可以有效避免过于严重的舍入误差。
正交矩阵是性态最好的矩阵，它不改变矩阵的条件数，通过QR分解来计算线性方程组，也可以避免误差的放大，保证计算过程具有数值稳定性。
四、已知一组数据点
，编写一程序求解三次样条插值函数
满足
并针对下面一组具体实验数据
0.25
0.3
0.39
0.45
0.53
0.5000
0.5477
0.6245
0.6708
0.7280
求解，其中边界条件为
.
【程序代码，文件名Spline.m】
function[s]=Spline（x,y,f0,fn）
%输入实验数据x,y；边界二阶导数f0,fn
clc
figure
（1）
plot（x,y,'*r'）
holdon
N=max（size（x））;
symsts;
fork=1:
N-1;
h（k）=x（k+1）-x（k）;
end
g
（1）=3*（y
（2）-y
（1））/h
（1）-h
（1）/2*f0;
g（N）=3*（y（N）-y（N-1））/h（N-1）+h（N-1）/2*fn;
fork=2:
N-1
lamda（k）=h（k）/（h（k）+h（k-1））;
miu（k）=h（k-1）/（h（k）+h（k-1））;
g（k）=3*（miu（k）*（y（k+1）-y（k））/h（k）+...
lamda（k）*（y（k）-y（k-1））/h（k-1））;
end
c=[1,miu（2:
N-1）];
b=2*ones（1,N）;
a=[lamda（2:
N-1）,1];
A=diag（c,1）+diag（b,0）+diag（a,-1）;
d=c;
a=[0,a];
u
（1）=b
（1）;
fori=2:
N
l（i）=a（i）/u（i-1）;
u（i）=b（i）-l（i）*c（i-1）;
end
L=eye（N）+diag（l（2:
N）,-1）;
U=diag（u）+diag（d,1）;
z
（1）=g
（1）;
fori=2:
N
z（i）=g（i）-l（i）*z（i-1）;
end
m（N）=z（N）/u（N）;
fori=N-1:
-1:
1
m（i）=（z（i）-c（i）*m（i+1））/u（i）;
end
fori=1:
N-1
s（i）=（h（i）+2*（t-x（i）））*（t-x（i+1））^2*y（i）/（h（i））^3+...
（h（i）-2*（t-x（i+1）））*（t-x（i））^2*y（i+1）/（h（i））^3+...
（t-x（i））*（t-x（i+1））^2*m（i）/（h（i））^2+...
（t-x（i+1））*（t-x（i））^2*m（i+1）/（h（i））^2;
End
digits（8）
s=vpa（expand（s））%输出分段表达式
fori=1:
N-1%绘图
v=x（i）:
0.005:
x（i+1）;
y=subs（s（i）,v）;
plot（v,y）;
holdon;
end
gridon
[命令窗口输入]
x=[0.250.30.390.450.53];
y=[0.5000,0.5477,0.6245,0.6708,0.7280];
f0=0;fn=0;
Spline（x,y,f0,fn）
[运行结果]
ans=
[4.6988737*t^2-.20505552*t+.35547747-6.2651650*t^3,
-2.6329843*t^2+1.9945019*t+.13552173+1.8813439*t^3,
.10638708*t^2+.92614706*t+.27440786-.45999912*t^3,
-3.4093028*t^2+2.5082075*t+.37098798e-1+2.1442156*t^3]
【分析】
运行结果是一个四个元素的矩阵，各个元素依次代表四个分段区间上的三次样条曲线。
例如在第一个区间[0.250.3]上，拟合得到的三次样条曲线方程4.6988737*t^2-0.20505552*t+0.35547747-6.2651650*t^3。
由图像可知，三次样条曲线是很光滑的曲线，拟合效果很好。
五、编写程序构造区间
上的以等分结点为插值结点的Newton插值公式，假设结点数为
（包括两个端点），给定相应的函数值，插值区间和等分的份数，该程序能快速计算出相应的插值公式。
以
，
为例计算其对应的插值公式，分别取不同的
值并画出原函数的图像以及插值函数的图像，观察当
增大时的逼近效果.
【程序代码，文件名Newton.m】
functionf1=Newton（n）
clc
symsxtf
a=-1;b=1;
f=1./（1+25*x.^2）;
y=zeros（1,n）;
x=linspace（a,b,n）;
h=-1:
0.01:
1;
m=n;
c（1:
m）=0.0;
yh=subs（f,h）;
y=subs（f,x）;
yk=y;
f1=y
（1）;
y1=0;k=1;
fori=1:
m-1
forj=i+1:
m
y1（j）=（y（j）-y（i））/（x（j）-x（i））;
end
c（i）=y1（i+1）;
k=k*（t-x（i））;
f1=f1+c（i）*k;
y=y1;
end
f1=simplify（f1）
yfh=subs（f1,h）;
figure
（1）
plot（h,yfh,'--k',x,yk,'*r'）
gridon
holdon
x=-1:
0.01:
1;
y=1./（1+25*x.^2）;
plot（x,yh,'b'）
legend（'插值曲线','插值节点','被插曲线'）;
[运行结果]
Newton（4）（5个插值节点）
拟合方程为：
1-3225/754*t^2+1250/377*t^4
Newton（5）（6个插值节点）
拟合方程为：
（用digits（4）,vpa（ans）处理后）
.5673+.1599e-16*t-1.731*t^2-.4441e-15*t^3+1.202*t^4+.1110e-14*t^5
当n继续增大时有：
n=7n=10
n=12n=15
【分析】
本题说明在不分段拟合的条件下，节点并不是越不越好。
插值节点增多会导致插值多项式的次数增高，而高次多项式的振荡次数增多有可能使插值多项式在非节点处的误差变得很大（Runge现象）。
所以在高次插值时，为了克服这种现象，建议改为分段低次插值。
升级会员 