西安交大计算方法B大作业Word格式.docx
《西安交大计算方法B大作业Word格式.docx》由会员分享，可在线阅读，更多相关《西安交大计算方法B大作业Word格式.docx（27页珍藏版）》请在冰豆网上搜索。
（2）估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图;
1.2实现题目的思想及算法依据
首先在题目
（1）中要实现的是数据的拟合，显然用到的是我们在第三章中数据近似的知识内容。
多项式插值时，这里有21个数据点，则是一个20次的多项式，但是多项式插值随着数据点的增多，会导致误差也会随之增大，插值结果会出现龙格现象，所以不适用于该题目中点数较多的情况。
为了避免结果出现大的误差，同时又希望尽可能多地使用所提供的数据点，提高数据点的有效使用率，这里选择分段插值方法进行数据拟合。
分段插值又可分为分段线性插值、分段二次插值和三次样条插值。
由于题目中所求光缆的现实意义，而前两者在节点处的光滑性较差，因此在这里选择使用三次样条插值。
根据课本SPLINEM算法和TSS算法，采用第三种真正的自然边界条件，在选定边界条件和选定插值点等距分布后，可以先将数据点的二阶差商求出并赋值给右端向量d，再根据TSS解法求解三对角线线性方程组从而解得M值。
求出M后，对区间进行加密，计算200个点以便于绘图以及光缆长度计算。
对于问题
（2），使用以下的公式
（x）ds
1.3算法结构
1.Fori0,1,2,,n
1.1
yi
Mi
For
k
1,2
2.1
i
n,n
1,,k
2.1.
1（MiM
i1）/（XiXi
x
h1
1,2,
n-
4.1
X1
X
hi1
4.2
hi1/（hi
hi1）
Ci;
1C
4.3
6Mi
di
do
M0
;
dn
Mn;
0C0
b0;
n
an;
bn
1,d1
2,3,
m
!
获取1
7.1
aJ
k1
lk
7.2
bk-
lkC
7.3
dk-
m/
m
Mr
m1,m2,
1
9.1
（k
Ck
Mk1）/
kMk
Fori
12,
s1
11.1
if
Xithenik
break
ai；
2b
k）Mi
M的矩阵元素个数，存入m
获取x的元素个数存入s
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.xk
elsei1k
xk1
〜一〜
h;
x<
xx;
xx<
3x
xh2
h"
〜
[Mk
1_
Mk（yk1Mk1）x（yk
Mk—）xVhy
66
1.4matlab源程序
n=20;
x=O:
n;
y=[9.018.967.967.978.029.0510.1311.1812.2613.2813.3212.6111.2910.22
9.157.907.958.869.8110.8010.93];
M=y;
%用于存放差商，此时为零阶差商
h=zeros（1,n+1）;
c=zeros（1,n+1）;
d=zeros（1,n+1）;
a=zeros（1,n+1）;
b=2*ones（1,n+1）;
h
（2）=x
（2）-x
（1）;
fori=2:
n%书本110页算法SPLINEM
h（i+1）=x（i+1）-x（i）;
c（i）=h（i+1）/（h（i）+h（i+1））;
a（i）=1-c（i）;
end
a（n+1）=-2;
%计算边界条件c（0）,a（n+1）,采用的是第三类边界条件
c
（1）=-2;
fork=1:
3%计算k阶差商
fori=n+1:
-1:
k+1
M（i）=（M（i）-M（i-1））/（x（i）-x（i-k））;
if（k==2）%计算2阶差商
d（2:
n）=6*M（3:
n+1）;
%给d赋值
if（k==3）
d
（1）=（-12）*h
（2）*M（4）;
%计算边界条件d（0）,d（n）,采用的是第三类边界条件
d（n+1）=12*h（n+1）*M（n+1）;
l=zeros（1,n+1）;
r=zeros（1,n+1）;
u=zeros（1,n+1）;
q=zeros（1,n+1）;
u
（1）=b
（1）;
r
（1）=c
（1）;
q
（1）=d
（1）;
fork=2:
n+1%利用书本49页算法TSS求解三对角线性方程组
r（k）=c（k）;
l（k）=a（k）/u（k-1）;
u（k）=b（k）-l（k）*r（k-1）;
q（k）=d（k）-l（k）*q（k-1）;
p（n+1）=q（n+1）/u（n+1）;
fork=n:
p（k）=（q（k）-r（k）*p（k+1））/u（k）;
fprintf（'
三对角线性方程组的解为：
'
）;
disp（p）;
%求拟合曲线
x1=0:
0.1:
20;
%首先对区间进行加密，增加插值点
n1=10*n;
x2=zeros（1,n1+1）;
x3=zeros（1,n1+1）;
s=zeros（1,n1+1）;
fori=1:
n1+1
forj=1:
n
ifx1（i）>
=x（j）&
&
x1（i）<
=x（j+1）%利用书本111页算法EVASPLINE求解拟合
曲线s（x）
h（j+1）=x（j+1）-x（j）;
x2（i）=x（j+1）-x1（i）;
x3（i）=x1（i）-x（j）;
s（i）=（p（j）.*（x2（i））A3/6+p（j+1）.*（x3（i））A3/6+（y（j）-p（j）.*（（h（j+1））A2/6））.*x2（i）+…
（y（j+1）-p（j+1）.*（h（j+1））A2/6）.*x3（i））/h（j+1）;
plot（x,-y,'
x'
）%画出插值点
holdon
plot（x1,-s）%画出三次样条插值拟合曲线
title（'
三次样条插值法拟合电缆曲线'
xlabel（'
河流宽度/m'
ylabel（'
河流深度/m'
Length=0;
n1
L=sqrt（（x1（i+1）-x1（i））A2+（s（i+1）-s（i））A2）;
%计算电缆长度
Length=Length+L;
电缆长度（m）='
disp（Length）;
1.5结果与说明
图1.1三次样条插值法拟合海底光缆曲线
山舅10-0.70091,922&
0.8703-山24斫0.3520-0.9224-1,8224
电缆长l（n）=26.6656
图1.2海底光缆长度结果
铺设海底光缆的曲线如图1.1所示
由上图可以看出，所得到的曲线光滑，能够较好得反映实际的河沟底部地势
形貌。
电缆长度计算结果为26.6656m（图1.2）。
题目二
2.1题目内容：
已知非线性方程一0cos（xsin）d0在［2,3］中有根。
请设计算法，求出该根,
并使求出的根的误差不超过104。
2.2实现题目的思想及算法依据
对于该题的非线性方程，可以将其分解成两个部分：
（1）求解数值积分；
（2）求解非线性方程。
首先求解数值积分，令g（）cos（xsin），则利用最简单的梯形公式可以得到
nh
f（x）—［（g（i1）g（i）］，其中h—,iih。
于是就有了f（x）0形式的非i12n
线性方程，这里选择二分法进行求解。
算法参考课本BISECTION算法。
2.3算法结构
2.iff0f10thenstop
3.iff02then输出x0作为根；
stop
4.iff12then输出x作为根；
101
5.x+xx
6.Ifx1x1x1then输出x作为根；
7.fxf
8.iff2then输出x作为根；
9.if£
f0then
9.1x
Else
9.2xx1
10.goto5
2.4Matlab源程序
************************
functiony=theta（i）y=i*h;
*************************
functiony=g（x,theta）
%为了计算关于theta的数值积分，先令
g=cos（x*sin（theta））
y=cos（x*sin（theta））;
**************************
functionf=hsz（x）%计算数值积分
n=10000;
h=pi/n;
%将区间分成n份
f=0;
f=f+h/（2*pi）*（g（x,i+1）+g（x,i））;
error=1e-4;
a=2;
b=3;
f0=hsz（a）;
f1=hsz（b）;
iff0*f1>
%误差允许值
%初始区间
%判断方程是否有解
disp（'
该方程在［a,b］上无解'
elseiffO==O
x=a;
elseiff1==0
x=b;
%判断方程解是否在区间两边界上
else
%二分法求解方程得解
a0=a;
b0=b;
whileabs（（b0-a0）/2）>
=error
half=（a0+b0）/2;
fa=hsz（aO）;
fb=hsz（bO）;
fhalf=hsz（half）;
%计算中点处的函数值，用以判断解的位置
iffhalf==0
x=half;
break;
elseiffa*fhalf<
b0=half;
%定义新区间，为原区间的一半
a0=half;
%定义新区间
x=（b0+a0）/2;
%方程组的解
方程组的解为：
）
disp（x）
2.5结果与说明
n的取值是很关
由于是利用梯形公式来求解，需要将区间划分为n个区间，而键的，需要取得适当的n值才能满足误差精度。
»
th2
方程的解为；
2.4048
这里将区间分成1OOOO份，可以得到结果图2.1,x=2.4048
图2.1n=10000方程的解
由于变量的嵌套很多及自身水平的不足，所以我将很多变量通过调用子程序的方式实现，分别为theta.m:
g.m;
hsz.m。
3.1题目内容：
编写程序实现大规模方程组的高斯消去法程序，并对所附的方程组进行求解。
所附方程组的类型为对角占优的带状方程组。
数据说明：
⑴dat20171.dat为非压缩带状方程组，阶数为10阶，该方程组供调
试程序使用，该方程组的根都为1;
（2）dat20172.dat为压缩带状方程组，阶数为20阶，该方程组供调试
程序使用，该方程组的根都为1;
⑵dat20173.dat为非压缩带状方程组，阶数在3000阶左右；
⑶dat20174.dat为压缩带状方程组，阶数在50000阶左右；
3.2实现题目的思想及算法依据
高斯消去法主要分为消去和回代过程，消去过程形成上三角矩阵，回代过程就是自下而上实现方程组求解。
题目中系数矩阵是严格对角占优矩阵，所以无需进行列主元，极大的提高了程序的效率。
对于n阶、上带宽q、下带宽p的非压缩格式带状矩阵，通过比较k+p、k+q及n的大小，对于行消去，选择k+p和n中的最小值作为循环终点。
在用第k行进行消去时，只需要计算第k列的ak+1，k到3k+p，k，每一行只需计算ak，k+1到ak，k+q；
从第n-c+1到n行，第k列计算ak+1,k到an，k，每一行只需计算ak，k+1到3k，n，于是减小运算量。
压缩格式忽略了零元素，存储方式与非压缩格式有所不同，所以需要建立两种格式中元素之间的对应关系，即压缩Ak，p+1=非压缩人，k，根据此关系进行程序编写。
3.3算法结构
1.A的阶数n
Forik1,k2,…,n
2.2Forjk1,k2,…,n
2.3
nn
Xn
3.Forkn1,n2,…,1
（kkjXj）/kkxk
jk1
3.4Matlab源程序
%首先计算非压缩带状方程组
fname='
dat20171.dat：
%读取文件，可改为dat20173.dat
目标文件为:
disp（fname）;
fid=fopen（fname,'
r'
%二进制模式读取数据
fhead=fread（fid,5,'
long'
%以longint形式读取文件前5个数据，获取系数矩阵信息
bbh=dec2hex（fhead
（2））;
%将版本号转换为16进制
%判断所读取的文件是否为目标文件
ifstrcmp（bbh,'
1O2'
%比较版本号，相同则说明为非压缩带状矩阵
目标文件系数矩阵为非压缩带状矩阵，读取正确'
读取错误'
n=fhead（3）;
%读取矩阵阶数n
q=fhead（4）;
%读取上带宽q
p=fhead（5）;
%读取下带宽p
系数矩阵阶数为:
disp（n）
上带宽为:
disp（q）
下带宽为:
d=fread（fid,nA2,'
float'
b=fread（fid,n,'
%生成系数矩阵
A=zeros（n,n）;
k=1;
A（i,j）=d（k）;
（i,j）
k=k+1;
fclose（fid）;
%非压缩格式，需要读取nA2个浮点数,
%读取右端系数
%因为系数矩阵按行方式顺序存储在d
%读取结束
n=length（b）;
按行方式顺序存储
中，所以依次形成A
%高斯消去法，书本33页到34页算法GAUSSPP
%消去过程
n-1
fori=k+1:
min（k+p,n）
%仅对不为零的区域作高斯消去，边界时k+p有可能大
A（i,k）=A（i,k）/A（k,k）;
%用第k行消去
forj=k+1:
min（k+q,n）
%只需处理右边k+q个数据，边界时k+q可能大于n
A（i,j）=A（i,j）-A（i,k）*A（k,j）;
%更新各列
b（i）=b（i）-A（i,k）*b（k）;
%更新右端系数b列
%回代过程
x（n）=b（n）/A（n,n）;
fork=n-1:
sum=0;
sum=sum+A（k,j）*x（j）;
x（k）=（b（k）-sum）/A（k,k）;
endm=10;
方程的前％d个根:
\n'
m）;
%输出前m个根
fprintf（1,'
%0.5f'
x（j））;
%小数点后保留5位小数
%将方程出的解输出为txt文件
ifstrcmp（fname,'
dat20173.dat'
file=strcat（'
f:
\\'
fname,'
矩阵方程组的解'
'
.txt'
%生成要创建文件的文件路径和文件
名
fid=fopen（file,'
a+'
%以读写方式打开指定文件，将文件指针指向文件末尾
fprintf（fid,'
%s的所有解\n'
fname）;
x）;
%向指定文件写入数据，保留5位小数
endfprintf（'
%第二步处理压缩带状方程组，
即数据文件
（2）和（4），
（2）供调试使用
dat20172.dat'
%读取文件，可改为dat20174.dat
2O2'
%比较版本号，相同则说明文件数据是压缩带状矩阵
目标文件系数矩阵为压缩带状矩阵，读取正确'
系数矩阵阶数为disp（n）
d=fread（fid,n*（p+q+1）,'
%压缩格式，以n*（p+q+1）矩阵存储，同样按行存储
b=fread（fid,n,'
A=zeros（n,p+q+1）;
k=1;
p+q+1
%因为系数矩阵按行方式顺序存储在d中，所以依次形成A
%读取结束
min（p,n-k）%用第k行消去
A（k+i,p+1-i）=A（k+i,p+1-i）/A（k,p+1）;
%压缩格式p+1列即为原来第k列
min（q,n-k）%只需处理右边q个数据，但边界时q可能大于
（n-k）
A（k+i,p+1-i+j）=A（k+i,p+1-i+j）-A（k+i,p+1-i）*A（k,p+1+j）;
b（k+i）=b（k+i）-b（k）*A（k+i,p+1-i）;
x（n）=b（n）/A（n,p+1）;
fori=n-1:
sum=O;
min（q,n-i）
sum=sum+x（i+j）*A（i,p+1+j）;
x（i）=（b（i）-sum）/A（i,p+1）;
m=20;
方程的前%d个根:
dat20174.dat'
3.5结果与说明
目标dat
目标女件系埶年阵淘非压缩帯状拒阵」读取IF确
系数世解阶數为10
上带宽为：
3
下帯克为.3
方程的前汕吓更
1.000001,COOOO1.OOCOOL000D0LOOOCO
目标文件为：
dat^Ul,v.dat
目标交件系数柜阵为圧缩帚我矩陳，读則正确
垂皴柜解撕數为20
1.00000
1.D0000
LdOOOOL00000
下带竞为：
方程的前创个艰
1,00000l.COOOOL.000001.C00501.OOOCO
1.00000
1.000001.000001.04000
首先使用dat20171和dat20172验证程序的正确性,
如图
3.1
图3.1验证dat20171和dat20172的正确性
由上图可知
dat20171.dat中为上带宽为3，下带宽为3，阶数为10的非压缩矩阵，解都为
dat20172.dat中为上带宽为3，下带宽为3，阶数为20的压缩矩阵，解都为1
与题目中所给的条件相符
图3.2dat20173和dat20174的方程组解
目棕文件为:
dat201T3.dat
目标文件系数矩阵为非压縮带状矩0车・读现止确
慕埶拒晦翫數为:
29E0
上帯宽为：
4
下带竟为：
方程的前1Q吓很.
1.019001.013001.01900L019001.019001.01S&
01.019001.019001.019&
01.01900
目标
升级会员 