西安交大计算方法B大作业Word格式.docx
- 文档编号:17755881
- 上传时间:2022-12-09
- 格式:DOCX
- 页数:27
- 大小:288.92KB
西安交大计算方法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
目标
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 西安 交大 计算方法 作业