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

《计算方法B》上机

实验报告

学院：班级: 姓名：学号：

机械工程学院

2015年12月22日1

1.计算以下和式：S??1nn?016?211??4????8n?18n?48n?58n?6?，要求： ??(1)若保留11个有效数字，给出计算结果，并评价计算的算法；

(2)若要保留30个有效数字，则又将如何进行计算。

实现思想：

以上问题出现了近似数相减的问题，为了减小误差，可分别求得减数之和以及被减数之和，最后将两者相减。另外，减数与被减数求和均为同号计算，按照绝对值递增顺序相加可减小舍入误差。此题中对有效数字有要求，因而计算时首先需要根据有效数字位数计算得出迭代次数，以保证计算值的精度。

源程序：

m=input('输入有效数字个数m=');

s0=1;s1=0;s2=0;n=0; %判断迭代次数

while s0>=0.5*10^-(m-1)

s0=4/(16^n*(8*n+1))-2/(16^n*(8*n+4))-1/(16^n*(8*n+5))-1/(16^n*(8*n+6)); n=n+1; end

%分别求解各项并求和 for k=n-1:-1:0

a1=4/(16^k*(8*k+1)); a2=2/(16^k*(8*k+4)); a3=1/(16^k*(8*k+5)); a4=1/(16^k*(8*k+6)); s1=a1+s1; s2=a4+a3+a2+s2; end

S=vpa(s1-s2,m)

2

实验结果：11位有效数字计算结果如图1所示；30为有效数字计算结果如图2所示。

图1.11位有效数字计算结果

图2.30为有效数字计算结果

3

1. 某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿

直线走向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。已探测到一组等分点位置的深度数据(单位:米)如下表所示：

分点 深度 分点 深度 分点 深度 0 9.01 7 11.18 14 9.15 1 8.96 8 12.26 15 7.90 2 7.96 9 13.28 16 7.95 3 7.97 10 13.32 17 8.86 4 8.02 11 12.61 18 9.81 5 9.05 12 11.29 19 10.80 6 10.13 13 10.22 20 10.93 (1)请用合适的曲线拟合所测数据点；

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

算法思想：由于题中所给点数为20，若采用高次多项式插值将产生很大的误差，所以拉格朗日或牛顿并不适用。题中光缆为柔性，可光滑铺设于水底，鉴于此特性，采用三次样条插值插值法较为合适。

算法结构：三次样条算法结构见《计算方法教程》P110；

光缆长度计算公式：

20l=源程序： clear; clc; x=0:20;

?01+(f(x))dx?'2??k?k?0019k?120'1+(f(x))2dx

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]; d=y;

plot(x,y,'k.','markersize',15) hold on %%%计算差商

for k=1:2

for i=21:-1:(k+1)

d(i)=(d(i)-d(i-1))/(x(i)-x(i-k)); end end

%%%设定d的边界条件 for i=2:20

4

d(i)=6*d(i+1); end d(1)=0; d(21)=0;

%%%带状矩阵求解（追赶法） a=0.5*ones(1,21); b=2*ones(1,21); c=0.5*ones(1,21); a(1)=0; c(21)=0; u=ones(1,21); u(1)=b(1); r=c; yy(1)=d(1); %%%追

for k=2:21

l(k)=a(k)/u(k-1); u(k)=b(k)-l(k)*r(k-1); yy(k)=d(k)-l(k)*yy(k-1); end

%%%赶

m(21)=yy(21)/u(21); for k=20:-1:1

m(k)=(yy(k)-r(k)*m(k+1))/u(k); end %%%绘制曲线 k=1; nn=100;

xx=linspace(0,20,nn); l=0; for j=1:nn for i=2:20 if xx(j)<=x(i) k=i; break; else

k=i+1; end end h=1;

xbar=x(k)-xx(j); xmao=xx(j)-x(k-1);

s(j)=(m(k-1)*xbar^3/6+m(k)*xmao^3/6+(y(k-1)-m(k-1)*h^2/6)*xbar+(y(k)-m(k)*h^2/6)*xmao)/h;

5

结果对比：将二次函数、三次函数和四次函数拟合结果绘制在同一个坐标内，如图3-4所示。其计算误差结果见表3-1所示。

图3-4 拟合结果对比分析

11

4.设计算法，求出非线性方程6x5?45x2?20?0的所有实根，并使误差不超过10?4。

算法思想：本题可采用牛顿法迭代求解,令f(x)=6x5-45x2+20 ,得带格式为

xk+1=xk-f(xk)f(xk)'

根据函数图像可以找出根的大致分布区间,带入不同的初值即可解出不同的根.

源代码:

function y=f2(x)

y=6*x.^5-45*x.^2+20;%定义原函数 function y=f3(x)

y=30*x^4-90*x;%定义原函数倒数

i=-5:0.1:5; y=f2(i); plot(i,y) hold on

plot(i,0,'-')%画出原函数图像

%%Newton法求根 x1=input('输入初值'); e=10^(-4);%误差设定 Nmax=1000;%迭代最大次数限定 for n=1:Nmax f0=f2(x1); if abs(f2(x1))

fprintf('输出的f(x)已经足够小'); x=x1; break else

F0=f3(x1); x=x1-f0/F0; if abs(x-x1)

12

fprintf('输出方程的根x=/',x)

计算结果:函数图像如图4-1所示。计算结果分别见图4-2所示。

图4-1 函数图像

图4-2 计算结果

根据带入不同的初值，可以求出不同的根，有图4-2可以看出，原函数的根大约有三个，分别是-0.654542、0.681174、1.870799。

13

5.线性方程组求解。

（1）编写程序实现大规模方程组的高斯消去法程序，并对所附的方程组进行求解。所附方程组的类型为对角占优的带状方程组。 （2）针对本专业中所碰到的实际问题，提炼一个使用方程组进行求解的例子，并对求解过程进行分析、求解。 算法思想:

高斯消去法是利用现行方程组初等变换中的一种变换,将一个不为零的数乘到一个方程后加到另一个方程,使方程组变成同解的上三角方程组,然后再自下而上对上三角方程组求解。

算法结构：

1. 读取二进制文件，存入计算矩阵

2. 对矩阵进行初等变换，然后求解(详见计算方法教程第2版高斯消去法以及

列主元高斯消去法算法)

源代码：

clear; clc;

%% 读取系数矩阵

[f,p]=uigetfile('*.dat','选择数据文件'); %读取数据文件 num=5; %输入系数矩阵文件头的个数 name=strcat(p,f); file=fopen(name,'r');

head=fread(file,num,'uint'); %读取二进制头文件 id=dec2hex(head(1)); %读取标识符 fprintf('文件标识符为'); id

ver=dec2hex(head(2)); %读取版本号 fprintf('文件版本号为'); ver

n=head(3); %读取阶数 fprintf('矩阵A的阶数'); n

q=head(4); %上带宽

fprintf('矩阵A的上带宽');

14

q

p=head(5); %下带宽

fprintf('矩阵A的下带宽'); p

dist=4*num;

fseek(file,dist,'bof'); %把句柄值转向第六个元素开头处

[A,count]=fread(file,inf,'float'); %读取二进制文件，获取系数矩阵 fclose(file); %关闭二进制头文件 %% 对非压缩带状矩阵进行求解 if ver=='102', a=zeros(n,n); for i=1:n, for j=1:n,

a(i,j)=A((i-1)*n+j); %求系数矩阵a(i,j) end end

b=zeros(n,1); for i=1:n,

b(i)=A(n*n+i); end

for k=1:n-1, %列主元高斯消去法 m=k;

for i=k+1:n, %寻找主元

if abs(a(m,k))

if a(m,k)==0 %遇到条件终止 disp('错误！') return end

for j=1:n, %交换元素位置得主元 t=a(k,j);

a(k,j)=a(m,j); a(m,j)=t; t=b(k); b(k)=b(m); b(m)=t; end

for i=k+1:n, %计算l(i,k)并将其放到a(i,k)中 a(i,k)=a(i,k)/a(k,k); for j=k+1:n

a(i,j)=a(i,j)-a(i,k)*a(k,j); end

15

b(i)=b(i)-a(i,k)*b(k); end end

x=zeros(n,1); %回代过程 x(n)=b(n)/a(n,n); for k=n-1:-1:1,

x(k)=(b(k)-sum(a(k,k+1:n)*x(k+1:n)))/a(k,k); end end

%% 对压缩带状矩阵进行求解 if ver=='202', %高斯消去法 m=p+q+1;

a=zeros(n,m); for i=1:1:n for j=1:1:m

a(i,j)=A((i-1)*m+j); end end

b=zeros(n,1); for i=1:1:n

b(i)=A(n*m+i); %求b(i) end

for k=1:1:(n-1) %开始消去 if a(k,(p+1))==0 disp('错误！'); break; end st1=n; if (k+p)

for i=(k+1):1:st1

a(i,(k+p-i+1))=a(i,(k+p-i+1))/a(k,(p+1)); for j=(k+1):1:(k+q)

a(i,j+p-i+1)=a(i,j+p-i+1)-a(i,k+p-i+1)*a(k,j+p-k+1); end

b(i)=b(i)-a(i,k+p-i+1)*b(k); end end

x=zeros(n,1); %回代 x(n)=b(n)/a(n,p+1); sum=0;

for k=(n-1):-1:1 sum=b(k);

16

st2=n; if (k+q)

for j=(k+1):1:st2

sum=sum-a(k,j+p-k+1)*x(j); end

x(k)=sum/a(k,p+1); sum=0; end end

disp('方程组的的解为：') %输出解 disp(x’)

求解结果

对数据文件dat51求解，结果如下：

文件标识符为 id = F1E1D1A0 文件版本号为 ver = 102 矩阵A的阶数 n = 15

矩阵A的上带宽 q = 3

矩阵A的下带宽 p = 3

方程组的的解为：

1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000

对数据文件dat52求解，结果如下：

文件标识符为 id = F1E1D1A0 文件版本号为 ver = 202 矩阵A的阶数 n = 20

矩阵A的上带宽 q = 5

矩阵A的下带宽 p = 5

17

方程组的的解为：

1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000

对数据文件dat53求解，结果如下：

文件标识符为 id = F1E1D1A0 文件版本号为 ver = 102 矩阵A的阶数 n = 2160 矩阵A的上带宽 q = 5 矩阵A的下带宽 p = 5

方程组的解截图如图5-1所示（由于矩阵阶数较大，计算结果未能截取完整）。

图5-1数据文件dat53求解结果

对数据文件dat53求解，结果如下：

文件标识符为 id = F1E1D1A0 文件版本号为 ver = 202 矩阵A的阶数 n = 43240 矩阵A的上带宽

18

q = 4 矩阵A的下带宽 p = 4

方程组的解截图如图5-2所示（由于矩阵阶数较大，计算结果未能截取完整）。

图5-2 数据文件dat54求解结果

19

本文来源：https://www.bwwdw.com/article/r085.html