西安交大计算方法B大作业.docx
- 文档编号:4807747
- 上传时间:2022-12-09
- 格式:DOCX
- 页数:27
- 大小:288.92KB
西安交大计算方法B大作业.docx
《西安交大计算方法B大作业.docx》由会员分享,可在线阅读,更多相关《西安交大计算方法B大作业.docx(27页珍藏版)》请在冰豆网上搜索。
西安交大计算方法B大作业
计算方法B上机报告
姓名:
学号:
班级:
学院:
任课教师:
2017年12月29日
题目一:
1.1题目内容
某通信公司在一次施工中,需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。
在铺设光缆之前需要对沟底的地形进行初步探测,从而估计所需光缆的长度,为工程预算提供依据。
已探测到一组等分点位置的深度数据(单位:
米)如下表所示:
分点
0
1
2
3
4
5
6
深度
9.01
8.96
7.96
7.97
8.02
9.05
10.13
分点
7
8
9
10
11
12
13
深度
11.18
12.26
13.28
13.32
12.61
11.29
10.22
分点
14
15
16
17
18
19
20
深度
9.15
7.90
7.95
8.86
9.81
10.80
10.93
(1)请用合适的曲线拟合所测数据点;
(2)估算所需光缆长度的近似值,并作出铺设河底光缆的曲线图;
1.2实现题目的思想及算法依据
首先在题目
(1)中要实现的是数据的拟合,显然用到的是我们在第三章中数据近似的知识内容。
多项式插值时,这里有21个数据点,则是一个20次的多项式,但是多项式插值随着数据点的增多,会导致误差也会随之增大,插值结果会出现龙格现象,所以不适用于该题目中点数较多的情况。
为了避免结果出现大的误差,同时又希望尽可能多地使用所提供的数据点,提高数据点的有效使用率,这里选择分段插值方法进行数据拟合。
分段插值又可分为分段线性插值、分段二次插值和三次样条插值。
由于题目中所求光缆的现实意义,而前两者在节点处的光滑性较差,因此在这里选择使用三次样条插值。
根据课本SPLINEM算法和TSS算法,采用第三种真正的自然边界条件,在选定边界条件和选定插值点等距分布后,可以先将数据点的二阶差商求出并赋值给右端向量d,再根据TSS解法求解三对角线线性方程组从而解得M值。
求出M后,对区间进行加密,计算200个点以便于绘图以及光缆长度计算。
对于问题
(2),使用以下的公式
20
0
(x)ds
1.3算法结构
1.Fori0,1,2,,n
1.1
yi
Mi
For
k
1,2
2.1
For
i
n,n
1,,k
2.1.
1(MiM
i1)/(XiXi
x
h1
For
i
1,2,
n-
1
4.1
X1
X
hi1
4.2
hi1/(hi
hi1)
Ci;1C
4.3
6Mi
1
di
do
M0
;dn
Mn;
0C0
2
b0;
n
an;2
bn
3
1,d1
1
For
k
2,3,
m
!
获取1
7.1
aJ
k1
lk
7.2
bk-
lkC
k1
k
7.3
dk-
lk
k1
k
m/
m
Mr
n
For
k
m1,m2,
1
9.1
(k
Ck
Mk1)/
kMk
1
k
!
Fori
12,
s1
11.1
if
X
Xithenik
break
ai;2b
k)Mi
M的矩阵元素个数,存入m
获取x的元素个数存入s
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.xk
elsei1k
x
xk1
〜一〜
h;x 3x xh2 h"〜 [Mk 1_ Mk(yk1Mk1)x(yk Mk—)xVhy 6 66 6 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)); end if(k==2)%计算2阶差商 d(2: n)=6*M(3: n+1);%给d赋值 end if(k==3) d (1)=(-12)*h (2)*M(4);%计算边界条件d(0),d(n),采用的是第三类边界条件 d(n+1)=12*h(n+1)*M(n+1); end end 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); end p(n+1)=q(n+1)/u(n+1); fork=n: -1: 1 p(k)=(q(k)-r(k)*p(k+1))/u(k); end 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); end end end plot(x,-y,'x')%画出插值点 holdon plot(x1,-s)%画出三次样条插值拟合曲线 holdon title('三次样条插值法拟合电缆曲线'); xlabel('河流宽度/m'); ylabel('河流深度/m'); Length=0; fori=1: n1 L=sqrt((x1(i+1)-x1(i))A2+(s(i+1)-s(i))A2);%计算电缆长度 Length=Length+L; end fprintf('电缆长度(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题目内容: 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作为根;stop 101 5.x+xx 2 6.Ifx1x1x1then输出x作为根;stop 7.fxf 8.iff2then输出x作为根; 9.if£f0then 9.1x Else 9.2xx1 10.goto5 2.4Matlab源程序 ************************ functiony=theta(i)y=i*h; end ************************* functiony=g(x,theta) %为了计算关于theta的数值积分,先令 g=cos(x*sin(theta)) y=cos(x*sin(theta)); end ************************** functionf=hsz(x)%计算数值积分 n=10000; h=pi/n;%将区间分成n份 f=0; fori=1: n f=f+h/(2*pi)*(g(x,i+1)+g(x,i)); end end ************************** error=1e-4;a=2; b=3; f0=hsz(a); f1=hsz(b);iff0*f1>0 %误差允许值 %初始区间 %判断方程是否有解 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<0 b0=half;%定义新区间,为原区间的一半 else a0=half;%定义新区间 end end x=(b0+a0)/2;%方程组的解 end fprintf('方程组的解为: ') 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 2. Forik1,k2,…,n 2.2Forjk1,k2,…,n 2.3 nn Xn 3.Forkn1,n2,…,1 n (kkjXj)/kkxk jk1 3.4Matlab源程序 %首先计算非压缩带状方程组 fname='dat20171.dat: %读取文件,可改为dat20173.dat fprintf('目标文件为: '); disp(fname); fid=fopen(fname,'r'); %二进制模式读取数据 fhead=fread(fid,5,'long'); %以longint形式读取文件前5个数据,获取系数矩阵信息 bbh=dec2hex(fhead (2)); %将版本号转换为16进制 %判断所读取的文件是否为目标文件 ifstrcmp(bbh,'1O2') %比较版本号,相同则说明为非压缩带状矩阵 disp('目标文件系数矩阵为非压缩带状矩阵,读取正确') else disp('读取错误') end n=fhead(3); %读取矩阵阶数n q=fhead(4); %读取上带宽q p=fhead(5); %读取下带宽p fprintf('系数矩阵阶数为: '); disp(n) fprintf('上带宽为: '); disp(q) fprintf('下带宽为: '); disp(q) d=fread(fid,nA2,'float');b=fread(fid,n,'float'); %生成系数矩阵 A=zeros(n,n); k=1; fori=1: n forj=1: n A(i,j)=d(k); (i,j) k=k+1; end end fclose(fid); %非压缩格式,需要读取nA2个浮点数, %读取右端系数 %因为系数矩阵按行方式顺序存储在d %读取结束 n=length(b); 按行方式顺序存储 中,所以依次形成A %高斯消去法,书本33页到34页算法GAUSSPP %消去过程 fork=1: 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);%更新各列 end b(i)=b(i)-A(i,k)*b(k);%更新右端系数b列 end end %回代过程 x(n)=b(n)/A(n,n); fork=n-1: -1: 1 sum=0; forj=k+1: n sum=sum+A(k,j)*x(j); end x(k)=(b(k)-sum)/A(k,k); endm=10; fprintf('方程的前%d个根: \n',m);%输出前m个根 forj=1: m fprintf(1,'%0.5f',x(j));%小数点后保留5位小数 end %将方程出的解输出为txt文件 ifstrcmp(fname,'dat20173.dat') file=strcat('f: \\',fname,'矩阵方程组的解','.txt');%生成要创建文件的文件路径和文件 名 fid=fopen(file,'a+');%以读写方式打开指定文件,将文件指针指向文件末尾 fprintf(fid,'%s的所有解\n',fname); fprintf(fid,'%0.5f',x);%向指定文件写入数据,保留5位小数 fclose(fid); endfprintf('\n'); %第二步处理压缩带状方程组, 即数据文件 (2)和(4), (2)供调试使用 fname='dat20172.dat'; %读取文件,可改为dat20174.dat fprintf('目标文件为: '); disp(fname); fid=fopen(fname,'r'); %二进制模式读取数据 fhead=fread(fid,5,'long'); %以longint形式读取文件前5个数据,获取系数矩阵信息 bbh=dec2hex(fhead (2)); %将版本号转换为16进制 %判断所读取的文件是否为目标文件 ifstrcmp(bbh,'2O2') %比较版本号,相同则说明文件数据是压缩带状矩阵 disp('目标文件系数矩阵为压缩带状矩阵,读取正确') else disp('读取错误') end n=fhead(3); %读取矩阵阶数n q=fhead(4); %读取上带宽q p=fhead(5); %读取下带宽p fprintf('系数矩阵阶数为disp(n) fprintf('上带宽为: '); disp(q) fprintf('下带宽为: '); disp(q) d=fread(fid,n*(p+q+1),'float'); %压缩格式,以n*(p+q+1)矩阵存储,同样按行存储 b=fread(fid,n,'float'); %读取右端系数 %生成系数矩阵 A=zeros(n,p+q+1);k=1; fori=1: n forj=1: p+q+1 A(i,j)=d(k);%因为系数矩阵按行方式顺序存储在d中,所以依次形成A (i,j) k=k+1; end end fclose(fid);%读取结束 n=length(b); %高斯消去法,书本33页到34页算法GAUSSPP %消去过程 fork=1: n-1 fori=1: min(p,n-k)%用第k行消去 A(k+i,p+1-i)=A(k+i,p+1-i)/A(k,p+1);%压缩格式p+1列即为原来第k列 forj=1: 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);%更新各列 end b(k+i)=b(k+i)-b(k)*A(k+i,p+1-i);%更新右端系数b列 end end %回代过程 x(n)=b(n)/A(n,p+1); fori=n-1: -1: 1 sum=O; forj=1: min(q,n-i) sum=sum+x(i+j)*A(i,p+1+j); end x(i)=(b(i)-sum)/A(i,p+1); end m=20; fprintf('方程的前%d个根: \n',m);%输出前m个根 forj=1: m fprintf(1,'%0.5f',x(j));%小数点后保留5位小数 end %将方程出的解输出为txt文件 ifstrcmp(fname,'dat20174.dat') file=strcat('f: \\',fname,'矩阵方程组的解','.txt');%生成要创建文件的文件路径和文件 名 fid=fopen(file,'a+');%以读写方式打开指定文件,将文件指针指向文件末尾 fprintf(fid,'%s的所有解\n',fname); fprintf(fid,'%0.5f',x);%向指定文件写入数据,保留5位小数 fclose(fid); end 3.5结果与说明 目标dat 目标女件系埶年阵淘非压缩帯状拒阵」读取IF确 系数世解阶數为10 上带宽为: 3 下帯克为.3 方程的前汕吓更 1.000001,COOOO1.OOCOOL000D0LOOOCO 目标文件为: dat^Ul,v.dat 目标交件系数柜阵为圧缩帚我矩陳,读則正确 垂皴柜解撕數为20 1.00000 1.00000 1.D0000 LdOOOOL00000 上带宽为: 3 下带竞为: 3 方程的前创个艰 1,00000l.COOOOL.000001.C00501.OOOCO 1.00000 1.00000 1.00000 1.000001.000001.04000 首先使用dat20171和dat20172验证程序的正确性, 如图 3.1 图3.1验证dat20171和dat20172的正确性 由上图可知 dat20171.dat中为上带宽为3,下带宽为3,阶数为10的非压缩矩阵,解都为 1 dat20172.dat中为上带宽为3,下带宽为3,阶数为20的压缩矩阵,解都为1 与题目中所给的条件相符 图3.2dat20173和dat20174的方程组解 目棕文件为: dat201T3.dat 目标文件系数矩阵为非压縮带状矩0车・读现止确 慕埶拒晦翫數为: 29E0 上帯宽为: 4 下带竟为: 4 方程的前1Q吓很. 1.019001.013001.01900L019001.019001.01S&01.019001.019001.019&01.01900 目标
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 西安 交大 计算方法 作业