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

计算方法B上机报告

姓名：

学号：

班级：

学院：

任课教师：

2017年12月29日

题目一：

1.1

题目内容

某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。已探测到一组等分点位置的深度数据(单位:

(1)(2)估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图；

1.2 实现题目的思想及算法依据

首先在题目（1）中要实现的是数据的拟合，显然用到的是我们在第三章中数据近似的知识内容。多项式插值时，这里有21个数据点，则是一个20次的多项式，但是多项式插值随着数据点的增多，会导致误差也会随之增大，插值结果会出现龙格现象，所以不适用于该题目中点数较多的情况。为了避免结果出现大的误差，同时又希望尽可能多地使用所提供的数据点，提高数据点的有效使用率，这里选择分段插值方法进行数据拟合。分段插值又可分为分段线性插值、分段二次插值和三次样条插值。由于题目中所求光缆的现实意义，而前两者在节点处的光滑性较差，因此在这里选择使用三次样条插值。 根据课本SPLINEM 算法和TSS 算法，采用第三种真正的自然边界条件，在选定边界条件和选定插值点等距分布后，可以先将数据点的二阶差商求出并赋值给右端向量d ，再根据TSS 解法求解三对角线线性方程组从而解得M 值。求出M 后，对区间进行加密，计算200个点以便于绘图以及光缆长度计算。 对于问题（2），使用以下的公式:

20

=()L f x ds ?

20

(f x =? 19

1

0(k k

k f x +==∑?

1.3 算法结构

1.For n i ,,2,1,0???=

1.1 i i M y ?

2. For 2,1=k

2.1 For k n n i ,,1, -=

2.1.1 i k i i i i M x x M M ?----)/()(1

3. 101h x x ?-

4. For 1-,,2,1n i =

4.1 11++?-i i i h x x

4.2 b a c c h h h i i i i i i ??-?+++2;1;)/(11 4.3 i i d M ?+16

5. 0000;;c M d M d n n ???λ

n n n b a b ???2;;20μ

6. 1111,γμ??d b

7. For m k ,,3,2 =! 获取M 的矩阵元素个数，存入m

7.1

k k k l a ?-1/μ

7.2 k k k k c l b μ??-1- 7.3 k k k k l d γγ??-1- 8. m m m M ?μγ/

9. For 1,,2,1 --=m m k

9.1 k k k k k M M c ??-+μγ/)(1 10. k ?1! 获取x 的元素个数存入s 11. For 1,,2,1-=s i

11.1if i x x ≤~ then k i ?；break

else k i ?+1

12. x

x x x x x h x x k k k k ?~;~;

11?-?-?--- y h x h M y x h M y x M x M k k k k k k ~/]?)6

()6(6?6[2

211331

?-+-++---

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,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); end

a(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); %计算边界条件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);

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);

end

p(n+1)=q(n+1)/u(n+1);

for k=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);

for i=1:n1+1

for j=1:n

if x1(i)>=x(j)&&x1(i)<=x(j+1) %利用书本111页算法EV ASPLINE求解拟合曲线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)).^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

end

end

plot(x,-y,'x') %画出插值点

hold on

plot(x1,-s) %画出三次样条插值拟合曲线

hold on

title('三次样条插值法拟合电缆曲线');

xlabel('河流宽度/m');

ylabel('河流深度/m');

Length=0;

for i=1:n1

L=sqrt((x1(i+1)-x1(i))^2+(s(i+1)-s(i))^2); %计算电缆长度

Length=Length+L;

end

fprintf('电缆长度(m)=');

disp(Length);

1.5 结果与说明

铺设海底光缆的曲线如图1.1所示

图1. 1三次样条插值法拟合海底光缆曲线

由上图可以看出，所得到的曲线光滑，能够较好得反映实际的河沟底部地势形貌。电缆长度计算结果为26.6656m （图1.2）。

图1. 2海底光缆长度结果

题目二

2.1 题目内容： 已知非线性方程

0)sin cos(1

=?

π

θθπd x 在[2,3]中有根。请设计算法，求出该根，

并使求出的根的误差不超过410-。

2.2 实现题目的思想及算法依据

对于该题的非线性方程，可以将其分解成两个部分： （1）求解数值积分；（2）求解非线性方程。

首先求解数值积分，令()cos(sin )g x θθ=，则利用最简单的梯形公式可以得到

11()[(()()]2n

i i i h

f x

g g θθπ

-=≈+∑

，其中,i h ih n πθ==。于是就有了()0f x =形式的非

线性方程，这里选择二分法进行求解。算法参考课本BISECTION 算法。

2.3 算法结构

1. ()()00f x f ?()()

1

1f x f ?

2. if 010f f > then stop

3. if 02f ε< then 输出()0

x 作为根；stop

4. if 12f ε< then 输出()1

x 作为根；stop

5.

()()

()

0112

x +x x ? 6. If ()()11

1x x x ε-< then 输出x 作为根；stop

7. ()f x f ?

8. if 2f ε< then 输出x 作为根； 9. if 10f f < then 9.1 ()0

x x ? ;0f f ?

Else

9.2 ()1

x x ? ; 1f f ?

10. go to 5

2.4Matlab源程序

************************

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));

end

end

**************************

error=1e-4; %误差允许值

a=2;

b=3; %初始区间

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*fhalf<0

b0=half; %定义新区间，为原区间的一半else

a0=half; %定义新区间

end

end

x=(b0+a0)/2; %方程组的解

end

fprintf('方程组的解为：')

disp(x)

2.5结果与说明

由于是利用梯形公式来求解，需要将区间划分为n个区间，而n的取值是很关键的，需要取得适当的n值才能满足误差精度。

这里将区间分成10000份，可以得到结果图2.1，x=2.4048。

图2. 1n=10000方程的解

由于变量的嵌套很多及自身水平的不足，所以我将很多变量通过调用子程序的方式实现，分别为theta.m : g.m; hsz.m。

题目三

3.1 题目内容：

编写程序实现大规模方程组的高斯消去法程序，并对所附的方程组进行求解。所附方程组的类型为对角占优的带状方程组。 数据说明：

(1)dat20171.dat 为非压缩带状方程组，阶数为10阶，该方程组供调试程序使用,该方程组的根都为1;

(2)dat20172.dat 为压缩带状方程组，阶数为20阶，该方程组供调试程序使用,该方程组的根都为1;

(2)dat20173.dat 为非压缩带状方程组，阶数在3000阶左右; (3)dat20174.dat 为压缩带状方程组，阶数在50000阶左右;

3.2 实现题目的思想及算法依据

高斯消去法主要分为消去和回代过程，消去过程形成上三角矩阵，回代过程就是自下而上实现方程组求解。题目中系数矩阵是严格对角占优矩阵，所以无需进行列主元，极大的提高了程序的效率。

对于n 阶、上带宽q 、下带宽p 的非压缩格式带状矩阵，通过比较k+p 、k+q 及n 的大小，对于行消去，选择k+p 和n 中的最小值作为循环终点。在用第k 行进行消去时，只需要计算第k 列的k+1k a ，到

k+p k

a ，，每一行只需计算k k+1a ，到

k k+q

a ，；

从第n-c+1到n 行，第k 列计算k+1k a ，到n k a ，，每一行只需计算k k+1a ，到k n a ，，于是减小运算量。

压缩格式忽略了零元素，存储方式与非压缩格式有所不同，所以需要建立两种格式中元素之间的对应关系，即压缩k p+1

A ，=非压缩k k A ，，根据此关系进行程序编

写。

3.3算法结构 1.A 的阶数? n

2． For 12i k ,k ,,n =++… 2.1

ik kk

ik

ααα? 2.2For 12j k ,k ,,n =++… 2.2.1 ij kj ik ij αααα?-

2.3 i k ik i ββαβ?-

n nn n x ?αβ/

3. For 121k n ,n ,,=--…

∑+=?-

n

k j k kk j kj

k x x 1

/)(αα

β

3.4 Matlab 源程序

%首先计算非压缩带状方程组

fname='dat20171.dat'; %读取文件，可改为dat20173.dat fprintf('目标文件为:'); disp(fname);

fid=fopen(fname,'r'); %二进制模式读取数据

fhead=fread(fid,5,'long'); %以long int 形式读取文件前5个数据，获取系数矩阵信息 bbh=dec2hex(fhead(2)); %将版本号转换为16进制 %判断所读取的文件是否为目标文件

if strcmp(bbh,'102') %比较版本号，相同则说明为非压缩带状矩阵 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^2,'float'); %非压缩格式，需要读取n^2个浮点数，按行方式顺序存储 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; end

end

fclose(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 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列

end

end

%回代过程

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);

end

m=10;

fprintf('方程的前%d个根:\n',m); %输出前m个根

for j=1:m

fprintf(1,'%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); %向指定文件写入数据，保留5位小数

fclose(fid);

end

fprintf('\n');

%第二步处理压缩带状方程组，即数据文件（2）和（4），（2）供调试使用

fname='dat20172.dat'; %读取文件，可改为dat20174.dat

fprintf('目标文件为:');

disp(fname);

fid=fopen(fname,'r'); %二进制模式读取数据

fhead=fread(fid,5,'long'); %以long int形式读取文件前5个数据，获取系数矩阵信息

bbh=dec2hex(fhead(2)); %将版本号转换为16进制

%判断所读取的文件是否为目标文件

if strcmp(bbh,'202') %比较版本号，相同则说明文件数据是压缩带状矩阵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;

for i=1:n

for j=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

%消去过程

for k=1:n-1

for i=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列

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列

end

end

%回代过程

x(n)=b(n)/A(n,p+1);

for i=n-1:-1:1

sum=0;

for j=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个根

for j=1:m

fprintf(1,'%0.5f ',x(j)); %小数点后保留5位小数

end

%将方程出的解输出为txt文件

if strcmp(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 结果与说明

首先使用dat20171和dat20172验证程序的正确性，如图3.1

图3. 1 验证dat20171和dat20172的正确性

由上图可知

dat20171.dat中为上带宽为3，下带宽为3，阶数为10的非压缩矩阵，解都为1 dat20172.dat中为上带宽为3，下带宽为3，阶数为20的压缩矩阵，解都为1 与题目中所给的条件相符。

图3. 2 dat20173和dat20174的方程组解

接下来求解dat20172和dat20174数据，如图3.2

除此之外，还将程序运行的全部结果输出在txt文件中分别给出部分截图如图3.3

图3. 3 dat20173和dat20174输出txt文件中的结果

题目四

4.1题目内容

针对本专业中所碰到的实际问题，提炼一个使用方程组进行求解的例子，并对求解过程进行分析、求解。

对一维稳态无源项的对流扩散方程，取其边界条件为00L x==x=L =φφφφ，，，， 。试在x/L=0~1的范围内取15个节点，采用中心差分、一阶迎风、混合格式以及QUICK 格式，对1510P =,,△ 三种情形，画出00L -x φφφφ（）/(-)与/L 的变化图线 4.2解题思路

对一维稳态无源项的对流扩散方程采用中心差分、一阶迎风、混合格式和QUICK 格式的得到离散方程（数值传热学教材P173），之后根据其边界条件，进行求解

4.3matlba 源程序

Pa=input('请输入Pa=')%输入P △的值

x=0:1/15:1 Pe=15*Pa;

y=(exp(Pe*x)-1)/(exp(Pe)-1) plot(x,y,'-*k') %精确解 hold on

y(1)=0,y(16)=1; for i=2:15

y(i)=((1+0.5*Pa)*y(i-1)+(1-0.5*Pa)*y(i+1))/2; end

plot(x,y(1:16),'-or') %中心差分 hold on for i=2:15

y(i)=((1+Pa)*y(i-1)+y(i+1))/(2+Pa); end

plot(x,y(1:16),'-.>g') %一阶迎风 hold on for i=2:15 if Pa==1

y(i)=((1+0.5*Pa)*y(i-1)+(1-0.5*Pa)*y(i+1))/2; else y(i)=y(i-1) end end

plot(x,y(1:16),'-+y')

%混合格式

hold on

for i=2:15

if i==2

y(i)=y(i+1)/(2+Pa)+(1+Pa)*y(i-1)/(2+Pa)+(Pa/(2+Pa))*(6*y(i)-3*y(i-1)-3*y(i+1))/8

else

y(i)=y(i+1)/(2+Pa)+(1+Pa)*y(i-1)/(2+Pa)+(Pa/(2+Pa))*(5*y(i)-y(i-1)-y(i-2)-3*y(i+1))/8 end

end

plot(x, y(1:16),'-

%QUICK格式

hold on

title(['Pa=',num2str(Pa)])

legend('精确解','中心差分','一阶迎风','混合格式','QUICK格式')

4.4结果与说明

运行程序分别得到1510

的曲线图

P=,,

△

由图可以看出，QUICK和中心差分格式的计算精度较高，但两种格式都只是条件稳定；一阶迎风和混合格式绝对稳定，但一阶迎风的精度较低。