大连理工矩阵与数值分析上机作业.docx
《大连理工矩阵与数值分析上机作业.docx》由会员分享，可在线阅读，更多相关《大连理工矩阵与数值分析上机作业.docx（16页珍藏版）》请在冰点文库上搜索。
大连理工矩阵与数值分析上机作业
矩阵与数值分析
（上机作业）
Matrixandnumericalanalysis
学院（系）：
~~~~~~~~~~
学生姓名：
~~
指导老师：
~~~
学号：
~~~~~~~~~
完成日期：
2012.11.24
大连理工大学
DalianUniversityofTechnology
1.给定n阶方程组
，其中
，
则方程组有解
。
对
和
，分别用Gauss消去法和列主元消去法解方程组，并比较计算结果。
1）程序代码：
输入A，b
function[X,A,b,ep]=gauss（A,b）
n=length（A）;
%选主元
fork=1:
n-1
[~,p]=max（abs（A（k:
n,k）））;
p=p+k-1;
ifp>k
y=A（k,:
）;A（k,:
）=A（p,:
）;A（p,:
）=y;
y2=b（k）;b（k）=b（p）;b（p）=y2;
end
ifabs（A（k,k））<1e-10
disp（'高斯消去法求解失效'）
break
end
fori=k+1:
n
m=A（i,k）/A（k,k）;
A（i,k:
n）=A（i,k:
n）-m.*A（k,k:
n）;
b（i）=b（i）-m*b（k）;
end
end
%求解X
X（n）=b（n）/A（n,n）;
fork=n-1:
-1:
1
X（k）=（b（k）-sum（A（k,k+1:
n）.*X（k+1:
n）））/A（k,k）;
end
end
输入数组，A，b
调用函数[X,a,c,ep]=gauss（A,b）;
2）结果分析：
n=10，两种方法运行结果均近似为x=[1,1,…,1];
n=84，选主元的gauss消去法运算后的右下角元素变得很小，约为10^-24,矩阵变得开始病态。
最大误差为2.80*10^-6
不选主元，求解的最大误差为5.37*10^8,可见求解结果误差极大。
2.设有方程组
，其中
，
（a）选取不同的初始向量
和不同的右端向量
，给定迭代误差要求，用Jacobi和Gauss-Seidel迭代法计算，观测得出的迭代向量序列是否收敛。
若收敛，记录迭代次数，分析计算结果并得出你的结论。
（b）选定初始向量初始向量
和不同的右端向量
，如取
。
将
的主对角线元素成倍增长若干次，非主对角元素不变，每次用Jacobi法计算，要求迭代误差满足
，比较收敛速度，分析现象并得出你的结论。
1）程序代码：
function[X,k,t]=e2（A,b,ep,X0）
U=-triu（A,1）;
L=-tril（A,-1）;
D=diag（diag（A,0））;
B=D^-1*（L+U）;
f=D^-1*b;
%B=（D-L）^-1*U;
%f=（D-L）^-1*b;
k=0;
t=1;
whilet>ep
X=B*X0+f;
t=max（abs（X-X0））;
X0=X;
k=k+1;
ifk>1000
disp（'不收敛'）
break
end
end
end
2）结果分析：
情况一：
对于b=（1,1,…,1），X0=（0,0,…，0）;Jacobi迭代收敛，迭代次数20；G-S收敛，迭代次数13；X0=（1，2，3，…，10），Jacobi迭代收敛，迭代次数23；G-S收敛，迭代次数16
情况二：
对于b=（1，2，…，10），X0=（0,0,…，0）；Jacobi迭代收敛，迭代次数23；G-S收敛，迭代次数16；X0=（1，2，3，…，10），Jacobi迭代收敛，迭代次数22；G-S收敛，迭代次数15
在这两种情况下，两种方法均收敛，且G-S收敛速度较快
3.用迭代法求方程
的全部根，要求误差限为
。
1）程序代码：
二分法：
clear
t=1;
ep=0.5*10^-8;
%求根区间
a=-10000;b=10000;
m=0;p=3;c（10000）=0;cc（10000）=0;c
（1）=a;c
（2）=b;c（3）=（a+b）/2;d（6）=0;
whilem<3
m=0;
fori=1:
p-1
cc（2*i-1）=c（i）;
cc（2*i）=（c（i）+c（i+1））/2;
end
cc（2*i+1）=c（p）;
c=cc;
p=2*p-1;
fori=1:
p-1
if（f（c（i））*f（c（i+1））<0）
m=m+1;
d（2*m-1:
2*m）=[c（i）,c（i+1）];
end
end
t=t+1;
ift>1e3
disp（'error'）
break
end
end
%求根
fori=1:
3
a=d（2*i-1）;b=d（2*i）;
z=（a+b）/2;
whileabs（a-b）>2*ep
iff（z）*f（a）<0
b=z;
else
a=z;
end
z=（a+b）/2;
end
x（i）=（a+b）/2
end
2）结果分析：
通过二分法，先求根区间，再求根。
求得的结果为：
x=
-2.8794-0.65270.5321
4.给定数据如下表：
0
1
2
3
4
5
6
7
8
9
10
0.0
0.79
1.53
2.19
2.71
3.03
3.27
2.89
3.06
3.19
3.29
编制程序求三次样条插值函数在插值中点的样条函数插值，并作点集
和样条插值函数的图形，满足的边界条件为
（1）
（2）
1）程序代码：
y=[0.0;0.79;1.53;2.19;2.71;3.03;3.27;2.89;3.06;3.19;3.39];
x=[0;1;2;3;4;5;6;7;8;9;10];
%形成矩阵方程组
%边界条件1
%fori=1:
9
%g（i）=1.5*（（y（i+2）-y（i）））;
%end
%d=ones（1,8）.*0.5;
%A=diag（d,-1）+diag（d,1）+diag（2.*ones（1,9））;
%g
（1）=g
（1）-0.5*0.8;g（9）=g（9）-0.5*0.2;
%m（2:
10）=A^-1*g';
%m
（1）=0.8;m（11）=0.2;
%边界条件2
fori=2:
10
g（i）=1.5*（（y（i+1）-y（i-1）））;
end
g
（1）=3*（y
（2）-y
（1））;
g（11）=3*（y（11）-y（10））;
d=ones（1,10）.*0.5;
A=diag（d,-1）+diag（d,1）+diag（2.*ones（1,11））;
A（1,2）=1;A（10,9）=1;
m=A^-1*g';
fori=1:
10
s（1,i）=（x（i）+x（i+1））/2;
s（2,i）=y（i）*0.5+m（i）*（s（1,i）-x（i））*0.25+y（i+1）*0.5+m（i+1）*（s（1,i）-x（i+1））*0.25;
xx（i*100-99:
i*100）=linspace（x（i）,x（i+1）,100）;
dx=xx（i*100-99:
i*100）-ones（1,100）.*x（i）;
dx1=xx（i*100-99:
i*100）-ones（1,100）.*x（i+1）;
yy（i*100-99:
i*100）=y（i）.*（1+2*dx）.*dx1.^2+m（i）.*dx.*dx1.^2+y（i+1）.*（1-2.*dx1）.*dx.^2+m（i+1）.*dx1.*dx.^2;
end
plot（s（1,:
）,s（2,:
）,'o',xx,yy）
2）结果分析
边界条件1：
S=
Columns1through8
0.39861.16841.87152.47822.87323.21393.08352.9233
Columns9through10
3.13703.2823
边界条件2
S=
Columns1through8
0.39841.16851.87152.47822.87333.21373.08462.9193
Columns9through10
3.15183.2596
5.对下列数据作三次多项式拟合，取权系数
，给出拟合多项式的系数、平方误差，并作离散数据
和拟合多项式的图形。
-1.0
-0.5
0.0
0.5
1.0
1.5
2.0
-4.447
-0.452
0.551
0.048
-0.447
0.549
4.552
1）程序代码：
clear;
x=[-1.0,-0.5,0.0,0.5,1.0,1.5,2.0];
y=[-4.447,-0.452,0.551,0.048,-0.447,0.549,4.552];
A=zeros（4,4）;b（4）=0;
fori=1:
4
forj=1:
4
fork=1:
7
A（i,j）=A（i,j）+x（k）^（i+j-2）;
end
end
fork=1:
7
b（i）=b（i）+y（k）*x（k）^（i-1）;
end
end
a=A^-1*b';
xx=linspace（-1,2,100）;
yy=a
（1）+xx.*a
（2）+xx.^2.*a（3）+xx.^3.*a（4）;
yi=a
（1）+x.*a
（2）+x.^2.*a（3）+x.^3.*a（4）;
dy2=sum（（y-yi）.^2）
plot（x,y,'o',xx,yy）;
2）结果分析：
各项系数：
a=
0.5491
-0.0000
-2.9977
1.9991
平方误差
dy2=
2.1762e-005
图像：
6.用复化梯形公式、复化Simpson公式和复化Gauss-Legendre公式计算下列积分的近似值，使绝对误差限为
，并将计算结果与精确解作比较以及比较各种算法的计算量。
（1）
，
（2）
.
1）程序代码
%编写目标函数f
clc
clear
ep=0.5e-7;n=1;g=1;
while（abs（g-pi/4）>ep）
n=n+1;
x=linspace（0,1,n+1）;
dx=1/n;
%g=（f1（x
（1））+2.*sum（f1（x（2:
n）））+f1（x（n+1）））/（2*n）;
%k1=n+1;
%d1=abs（g-pi/4）;
%x50=x（2:
n+1）-ones（1,n）.*dx/2;
%g=（f1（x
（1））+4.*sum（f1（x50））+2*sum（f1（x（2:
n）））+f1（x（n+1）））/（6*n）;
%k2=2*n+1;
%d2=abs（g-pi/4）;x25=x（2:
n+1）-ones（1,n）.*dx*3/4;x50=x（2:
n+1）-ones（1,n）.*dx/2;x75=x（2:
n+1）-ones（1,n）.*dx/4;g=（7*f1（x
（1））+32*sum（f1（x25））+12*sum（f1（x50））+32*sum（f1（x75））+14*sum（f1（x（2:
n）））+7*f1（x（n+1）））/（90*n）;
d3=abs（g-pi/4）;
k3=4*n+1;
end
2）结果分析：
a）用梯形公式计算，积分间距，h=8.94^-4，达到要求精度，误差为4.99*10^-8，需要计算1120个函数值；
用simpson公式计算，积分间距h=0.0667,达到要求精度，误差为3.85*10^-8,需要计算31个函数值；
用Cotes公式计算，积分间距h=0.25,达到要求精度。
误差为1.37*10^-8，需要计算17个函数值；
b）用梯形公式计算，积分间距，h=0.0011，达到要求精度，误差为5.00*10^-8，需要计算914个函数值；
用simpson公式计算，积分间距h=0.25,达到要求精度，误差为3.78*10^-8,需要计算9个函数值；
用Cotes公式计算，积分间距h=0.3333,达到要求精度。
误差为1.10*10^-8，需要计算13个函数值；
综上可以看出，梯形公式计算收敛最慢，计算量非常大；应用cotes公式计算，收敛最快，但计算量并不总是最少，在b）中由于每个区间中需要计算三个函数值，计算量反而比simpson公式大，并且它代码复杂。
在两例中，simpson均表现较好的收敛性和较小的计算量，并且代码不很复杂。
7.用4阶Runge-Kutta法求解微分方程：
分别用
和
计算，并比较两个近似值与精确解
1）程序代码
%4阶Runge-Kutta法：
h=0.1;yy=zeros（1,2/h+1）;
yy
（1）=1;
fori=1:
（2/h）
yk=yy（i）;xk=（i-1）*h;
k1=-yk+cos（2*xk）-2*sin（2*xk）+2*xk*exp（-xk）;
k2=-（yk+k1*h/2）+cos（2*（xk+h/2））-2*sin（2*（xk+h/2））+2*（xk+h/2）*exp（-（xk+h/2））;
k3=-（yk+k2*h/2）+cos（2*（xk+h/2））-2*sin（2*（xk+h/2））+2*（xk+h/2）*exp（-（xk+h/2））;
k4=-（yk+k3*h）+cos（2*（xk+h））-2*sin（2*（xk+h））+2*（xk+h）*exp（-（xk+h））;
yy（i+1）=yk+h*（k1+2*k2+2*k3+k4）/6;
end
x=linspace（0,2,2/h+1）;
y=x.^2.*exp（-x）+cos（2.*x）;
dy=y-yy;
max（abs（dy））
plot（x,y,'-',x,yy,'--'）
2）结果分析
h1=0.1
最大误差为1.01*10^-6。
实线为真实值，圆圈为近似值
h2=0.05
最大误差为6.34*10^-8，可见步长缩短一倍，结果精度提高O（（h1/h2）4）倍。
实线为真实值，圆圈为近似值
升级会员 