西安交大计算方法B2017大作业

1、计算方法B上机报告姓 名： 学 号：班 级： 学 院： 任课教师： 2017年12月29日题目一：1.1 题目内容某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。已探测到一组等分点位置的深度数据(单位:米)如下表所示：分点0123456深度9.018.967.967.978.029.0510.13分点78910111213深度11.1812.2613.2813.3212.6111.2910.22分点14151617181920深度9.157.907.958.869.811

2、0.8010.93(1)请用合适的曲线拟合所测数据点；(2)估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图；1.2 实现题目的思想及算法依据 首先在题目（1）中要实现的是数据的拟合，显然用到的是我们在第三章中数据近似的知识内容。多项式插值时，这里有21个数据点，则是一个20次的多项式，但是多项式插值随着数据点的增多，会导致误差也会随之增大，插值结果会出现龙格现象，所以不适用于该题目中点数较多的情况。为了避免结果出现大的误差，同时又希望尽可能多地使用所提供的数据点，提高数据点的有效使用率，这里选择分段插值方法进行数据拟合。分段插值又可分为分段线性插值、分段二次插值和三次样条插值。由于题目中

3、所求光缆的现实意义，而前两者在节点处的光滑性较差，因此在这里选择使用三次样条插值。 根据课本SPLINEM算法和TSS算法，采用第三种真正的自然边界条件，在选定边界条件和选定插值点等距分布后，可以先将数据点的二阶差商求出并赋值给右端向量d，再根据TSS解法求解三对角线线性方程组从而解得M值。求出M后，对区间进行加密，计算200个点以便于绘图以及光缆长度计算。 对于问题（2），使用以下的公式:1.3 算法结构1. For 1.1 2. For 2.1 For 2.1.1 3. 4. For 4.1 4.2 4.3 5. 6. 7. For ! 获取M的矩阵元素个数，存入m7.1 7.2 7.3

4、8. 9. For 9.1 10. ! 获取x的元素个数存入s11. For 11.1 if then ；break else 12. 1.4 matlab源程序n=20;x=0:n;y=9.01 8.96 7.96 7.97 8.02 9.05 10.13 11.18 12.26 13.28 13.32 12.61 11.29 10.22 9.15 7.90 7.95 8.86 9.81 10.80 10.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,

5、n+1);h(2)=x(2)-x(1);for i=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);enda(n+1)=-2; %计算边界条件c(0),a(n+1),采用的是第三类边界条件c(1)=-2;for k=1:3 %计算k阶差商 for i=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

6、); %计算边界条件d(0),d(n),采用的是第三类边界条件 d(n+1)=12*h(n+1)*M(n+1); endendl=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);for k=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);endp(n+1)=q(n+1)/u(n+1);for k=n:-1:1

7、 p(k)=(q(k)-r(k)*p(k+1)/u(k);endfprintf('三对角线性方程组的解为：');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);for i=1:n1+1 for j=1:n if x1(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

8、); x3(i)=x1(i)-x(j); s(i)=(p(j).*(x2(i).3/6+p(j+1).*(x3(i).3/6+(y(j)-p(j).*(h(j+1).2/6).*x2(i)+. (y(j+1)-p(j+1).*(h(j+1).2/6).*x3(i)/h(j+1); end endendplot(x,-y,'x') %画出插值点hold onplot(x1,-s) %画出三次样条插值拟合曲线hold ontitle('三次样条插值法拟合电缆曲线');xlabel('河流宽度/m');ylabel('河流深度/m');

9、Length=0;for i=1:n1 L=sqrt(x1(i+1)-x1(i)2+(s(i+1)-s(i)2); %计算电缆长度 Length=Length+L;endfprintf('电缆长度(m)=');disp(Length);1.5 结果与说明图1. Error! Main Document Only.三次样条插值法拟合海底光缆曲线铺设海底光缆的曲线如图1.1所示图1. Error! Main Document Only. 海底光缆长度结果由上图可以看出，所得到的曲线光滑，能够较好得反映实际的河沟底部地势形貌。电缆长度计算结果为26.6656m（图1.2）。题目二2.

10、1 题目内容：已知非线性方程在2,3中有根。请设计算法，求出该根，并使求出的根的误差不超过。2.2 实现题目的思想及算法依据对于该题的非线性方程，可以将其分解成两个部分：（1）求解数值积分；（2）求解非线性方程。首先求解数值积分，令，则利用最简单的梯形公式可以得到，其中。于是就有了形式的非线性方程，这里选择二分法进行求解。算法参考课本BISECTION算法。2.3 算法结构1. 2. if then stop 3. if then 输出 作为根；stop4. if then 输出 作为根；stop5. 6. If then 输出x作为根；stop7.8. if then 输出x作为根；9. i

11、f then 9.1 ; Else9.2 ; 10. go to 52.4 Matlab源程序*function y = theta( i ) y=i*h;end*function y = g( x,theta )%为了计算关于theta的数值积分，先令g=cos(x*sin(theta) y=cos(x*sin(theta);end*function f = hsz(x) %计算数值积分n=10000; h=pi/n; %将区间分成n份f=0;for i=1:n f=f+h/(2*pi)*(g(x,i+1)+g(x,i);endend*error=1e-4; %误差允许值a=2; b=3;

12、%初始区间f0=hsz(a);f1=hsz(b); if f0*f1>0 %判断方程是否有解 disp('该方程在a,b上无解'); elseif f0=0 x=a;elseif f1=0 x=b; %判断方程解是否在区间两边界上else %二分法求解方程得解 a0=a; b0=b; while abs(b0-a0)/2)>=error half=(a0+b0)/2; fa=hsz(a0); fb=hsz(b0); fhalf=hsz(half); %计算中点处的函数值，用以判断解的位置 if fhalf=0 x=half; break; elseif fa*fha

13、lf<0 b0=half; %定义新区间，为原区间的一半 else a0=half; %定义新区间 end end x=(b0+a0)/2; %方程组的解endfprintf('方程组的解为：')disp(x) 2.5 结果与说明由于是利用梯形公式来求解，需要将区间划分为n个区间，而n的取值是很关键的，需要取得适当的n值才能满足误差精度。这里将区间分成10000份，可以得到结果图2.1，x=2.4048。图2. 1 n=10000方程的解由于变量的嵌套很多及自身水平的不足，所以我将很多变量通过调用子程序的方式实现，分别为theta.m : g.m; hsz.m。题目三3.

14、1 题目内容：编写程序实现大规模方程组的高斯消去法程序，并对所附的方程组进行求解。所附方程组的类型为对角占优的带状方程组。数据说明：(1)dat20171.dat 为非压缩带状方程组，阶数为10阶，该方程组供调试程序使用,该方程组的根都为1;(2)dat20172.dat 为压缩带状方程组，阶数为20阶，该方程组供调试程序使用,该方程组的根都为1;(2)dat20173.dat 为非压缩带状方程组，阶数在3000阶左右;(3)dat20174.dat 为压缩带状方程组，阶数在50000阶左右;3.2 实现题目的思想及算法依据高斯消去法主要分为消去和回代过程，消去过程形成上三角矩阵，回代过程就是

15、自下而上实现方程组求解。题目中系数矩阵是严格对角占优矩阵，所以无需进行列主元，极大的提高了程序的效率。对于n阶、上带宽q、下带宽p的非压缩格式带状矩阵，通过比较k+p、k+q及n的大小，对于行消去，选择k+p和n中的最小值作为循环终点。在用第k行进行消去时，只需要计算第k列的 到，每一行只需计算到；从第n-c+1到n行，第k列计算到，每一行只需计算到，于是减小运算量。压缩格式忽略了零元素，存储方式与非压缩格式有所不同，所以需要建立两种格式中元素之间的对应关系，即压缩=非压缩，根据此关系进行程序编写。3.3算法结构1. A的阶数 n2 For 2.1 2.2 For 2.2.1 2.3 3. F

16、or 3.4 Matlab源程序%首先计算非压缩带状方程组fname='dat20171.dat' %读取文件，可改为dat20173.datfprintf('目标文件为:');disp(fname);fid=fopen(fname,'r'); %二进制模式读取数据fhead=fread(fid,5,'long'); %以long int形式读取文件前5个数据，获取系数矩阵信息bbh=dec2hex(fhead(2); %将版本号转换为16进制%判断所读取的文件是否为目标文件if strcmp(bbh,'102')

17、 %比较版本号，相同则说明为非压缩带状矩阵 disp('目标文件系数矩阵为非压缩带状矩阵，读取正确')else disp('读取错误') endn=fhead(3); %读取矩阵阶数nq=fhead(4); %读取上带宽qp=fhead(5); %读取下带宽pfprintf('系数矩阵阶数为:');disp(n)fprintf('上带宽为:');disp(q)fprintf('下带宽为:');disp(q)d=fread(fid,n2,'float'); %非压缩格式，需要读取n2个浮点数，按行方式

18、顺序存储b=fread(fid,n,'float'); %读取右端系数%生成系数矩阵A=zeros(n,n);k=1;for i=1:n for j=1:n A(i,j)=d(k); %因为系数矩阵按行方式顺序存储在d中，所以依次形成A（i,j） k=k+1; endendfclose(fid); %读取结束n=length(b); %高斯消去法，书本33页到34页算法GAUSSPP%消去过程for k=1:n-1 for i=k+1:min(k+p,n) %仅对不为零的区域作高斯消去，边界时k+p有可能大于n A(i,k)=A(i,k)/A(k,k); %用第k行消去 for

19、 j=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列 endend%回代过程x(n)=b(n)/A(n,n);for k=n-1:-1:1 sum=0; for j=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个根for j=1:m fprintf(1,'%

20、0.5f ',x(j); %小数点后保留5位小数end%将方程出的解输出为txt文件if strcmp(fname,'dat20173.dat') file=strcat('f:',fname,'矩阵方程组的解','.txt');%生成要创建文件的文件路径和文件名 fid=fopen(file,'a+'); %以读写方式打开指定文件，将文件指针指向文件末尾 fprintf(fid,'%s的所有解n',fname); fprintf(fid,'%0.5f ',x); %向指定文

21、件写入数据，保留5位小数 fclose(fid);endfprintf('n');%第二步处理压缩带状方程组，即数据文件（2）和（4），（2）供调试使用fname='dat20172.dat' %读取文件，可改为dat20174.datfprintf('目标文件为:');disp(fname);fid=fopen(fname,'r'); %二进制模式读取数据fhead=fread(fid,5,'long'); %以long int形式读取文件前5个数据，获取系数矩阵信息bbh=dec2hex(fhead(2); %

22、将版本号转换为16进制%判断所读取的文件是否为目标文件if strcmp(bbh,'202') %比较版本号，相同则说明文件数据是压缩带状矩阵 disp('目标文件系数矩阵为压缩带状矩阵,读取正确')else disp('读取错误') endn=fhead(3); %读取矩阵阶数nq=fhead(4); %读取上带宽qp=fhead(5); %读取下带宽pfprintf('系数矩阵阶数为:');disp(n)fprintf('上带宽为:');disp(q)fprintf('下带宽为:');disp(

23、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;for i=1:n for j=1:p+q+1 A(i,j)=d(k); %因为系数矩阵按行方式顺序存储在d中，所以依次形成A（i,j） k=k+1; endendfclose(fid); %读取结束n=length(b); %高斯消去法，书本33页到34页算法GAUSSPP%消去过程for k=1:n-1 for i=1:m

24、in(p,n-k) %用第k行消去 A(k+i,p+1-i)=A(k+i,p+1-i)/A(k,p+1); %压缩格式p+1列即为原来第k列 for j=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列 endend%回代过程x(n)=b(n)/A(n,p+1);for i=n-1:-1:1 sum=0; for j=1:min(q,n-i) sum=s

25、um+x(i+j)*A(i,p+1+j); end x(i)=(b(i)-sum)/A(i,p+1);endm=20;fprintf('方程的前%d个根:n',m); %输出前m个根for j=1:m fprintf(1,'%0.5f ',x(j); %小数点后保留5位小数end%将方程出的解输出为txt文件if strcmp(fname,'dat20174.dat') file=strcat('f:',fname,'矩阵方程组的解','.txt');%生成要创建文件的文件路径和文件名 fid=fo

26、pen(file,'a+'); %以读写方式打开指定文件，将文件指针指向文件末尾 fprintf(fid,'%s的所有解n',fname); fprintf(fid,'%0.5f ',x); %向指定文件写入数据，保留5位小数 fclose(fid); end3.5 结果与说明首先使用dat20171和dat20172验证程序的正确性，如图3.1图3. 1 验证dat20171和dat20172的正确性由上图可知dat20171.dat中为上带宽为3，下带宽为3，阶数为10的非压缩矩阵，解都为1dat20172.dat中为上带宽为3，下带宽为3，

27、阶数为20的压缩矩阵，解都为1图3. Error! Main Document Only. dat20173和dat20174的方程组解与题目中所给的条件相符。接下来求解dat20172和dat20174数据，如图3.2图3. Error! Main Document Only. dat20173和dat20174输出txt文件中的结果除此之外，还将程序运行的全部结果输出在txt文件中分别给出部分截图如图3.3题目四4.1 题目内容针对本专业中所碰到的实际问题，提炼一个使用方程组进行求解的例子，并对求解过程进行分析、求解。对一维稳态无源项的对流扩散方程，取其边界条件为 。试在x/L=01的范围内取15个节点，采用中心差分、一阶迎风、混合格式以及QUICK格式，对 三种情形，画出 4.2 解题思路对一维稳态无源项的对流扩散方程采用中心差分、一阶迎风、混合格式和QUICK格式的得到离散方程（数值传热学教材P173），之后根据其边界条件，进行求解4.3 matl