西安交大 计算方法B上机作业

一、三次样条拟合

某通信公司在一次施工中，需要在水面宽度为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)估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图； 解：

1、算法实现的思想及依据

题目（1）为曲线拟合问题多项式插值、分段插值和最小二乘法。多项式插值，随着插值数据点的数目增多，误差也会随之增大，因此不选用。最小二乘法适于数据点较多的场合，在此也不适用。故选用分段插值。

分段插值又分为分段线性插值、分段二次插值、三次样条插值及更高阶的多项式插值。由本题的物理背景知，光缆正常工作时各点应该是平滑过渡，因此至少选用三次样条插值法。对于更高阶的多项式插值，由于“龙格现象”而不选用。

题目（2）求光缆长度，即求拟合曲线在0到20的长度，对弧长进行积分即可。光缆长度的第一型线积分表达式为lk019k1k1f'(x)2dx。

2、算法实现的结构

参考教材给出的SPLINEM算法和TTS算法，在选定边界条件和选定插值点等距分布后，可以先将数据点的二阶差商求出来并赋值给右端向量d，再根据TSS19k11f'(x)2dx。 解法求解M。光缆长度的第一型线积分表达式为lk0k3、程序运行结果及分析

图1.1三种边界条件下三次样条插值

图1.2光缆长度

4、MATLAB代码：

1）自己编程实现代码

clear;

clc;

I=input('你想使用第几种边界条件？请输入1、2、3之一: '); x=0:20;

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.8 10.93];

plot(x,-y,'k.','markersize',15)%y为深度，取负号 hold on

%% 计算一阶差商

y1=ones(1,21); for i=2:1:21

y1(i)=(y(i)-y(i-1))/(x(i)-x(i-1)); end

%% 计算二阶差商

y2=ones(1,21); for i=3:1:21

y2(i)=(y1(i)-y1(i-1))/(x(i)-x(i-2)); end

%% 计算三阶差商

y3=ones(1,21); for i=4:1:21

y3(i)=(y2(i)-y2(i-1))/(x(i)-x(i-3)); end

%% 选择边界条件（I）

if I==1

d(1)=0;d(21)=0;a(21)=0;c(1)=0;% 第一个点和最后一个点的二阶差商为0 end if I==2

d(1)=6*y1(1); d(21)=-6*y1(21); a(1)=1;c(1)=1; end if I==3

d(1)=-12*y3(1); d(21)=-12*y3(21); a(21)=-2;c(1)=-2;% end

for i=2:20

d(i)=6*y2(i+1); end

%% 构造带状矩阵求解（追赶法） b=2*ones(1,21);

a=0.5*ones(1,21);%a(21)=-2; c=0.5*ones(1,21);%c(1)=-2; u(1)=b(1);r(1)=c(1); %% 追 yz(1)=d(1); for i=2:21

l(i)=a(i)/u(i-1); u(i)=b(i)-l(i)*r(i-1); r(i)=c(i);

yz(i)=d(i)-l(i)*yz(i-1); end %% 赶

xg(21)=yz(21)/u(21); for i=20:-1:1

xg(i)=(yz(i)-r(i)*xg(i+1))/u(i); end

M=xg;%%所有点的二阶导数值 %% 求函数表达式并积分

t=1; h=1; N=1000

x1=0:20/(N-1):20; length=0; for i=1:N

for j=2:20

if x1(i)<=x(j) t=j; break; else

t=j+1; end end

f1=x(t)-x1(i); f2=x1(i)-x(t-1);

S(i)=(M(t-1)*f1^3/6/h+M(t)*f2^3/6/h+(y(t-1)-M(t-1)*h^2/6)*f1+(y(t)-M(t)*h^2/6)*f2)/h; Sp(i)=-M(t-1)*f1^2/2/h+M(t)*f2^2/2/h+(y(t)-y(t-1))/h-(M(t)-M(t-1))*h/6; length(i+1)=sqrt(1+Sp(i)^2)*(20/(N-1))+length(i);%第一类线积分 end

figure(1);

plot(x1,-S,'r-')%深度曲线

grid

disp(['第',num2str(I),'种边界条件下长度',num2str(length(N+1)),'米']) axis fill;

xlabel('测点/米');ylabel('深度/米'); title('三次样条曲线拟合'); legend('数据点','拟合曲线',3);

二、最小二乘近似

假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温变化的规律;试计算这一天的平均气温,并试估计误差。

时刻 0 1 2 3 4 5 6 7 8 9 10 11 12 平均气温 15 14 14 14 14 15 16 18 20 20 23 25 28 时刻 13 14 15 16 17 18 19 20 21 22 23 24 平均气温 31 34 31 29 27 25 24 22 20 18 17 16 解：

1、算法实现的思想及依据

本题共有25个数据点，数据量较大，因此选用多项式插值会导致较大的误差，插值样条函数虽然有较好的数值性质，但是Mi的计算量不小，而且没有统一的公式来表示。另外，本题要求找出温度变化规律，插值方法要求p(x)必须通已知数据点，只会导致拟合曲线失去本有的数据所表示的规律；从表中的数据点可以看出，温度并不准确（温度只测量到整数位），因此没有必要要求拟合曲线通过数据点。因此，选取最小二乘近似。

2、算法实现的结构

算法参考课本LSS算法。利用已有数据来生成了G，将y作为其最后一列（第n+1列）存放。先形成矩阵Qk，再变换Gk-1到Gk；求解三角方程

Rx=h1；最后根据定义计算误差。平均温度采用数据点相加求和求平均，避免数值积分的繁琐。

3、程序运行结果及分析

4、MATLAB代码

%%最小二乘法拟合气温变化规律 clear; clc; x=0: 24;

y=[15,14,14,14,14,15,16,18,20,20,23,25,28,31,34,31,29,27,25,24,22,20,18,17,16]; N=length(y);

sum=0; %求一天的平均温度average for i=1:N sum=sum+y(i); end

average=sum/N;

fprintf('平均温度T为：%f\\n',average);

h=input('请输入多项式近似函数最高项的次数h：'); n=h+1;

%%参课本LSS算法

G=zeros(N,n+1); %建立一个N行，n+1列的矩阵G %% 生成矩阵G各个元素

for i=1:n %矩阵G中前N列元素 for j=1:N

G(j,i)=x(1,j)^(i-1); end end

for i=1:1:N %将y作为矩阵G中第(N+1)列元素 G(i,n+1)=y(i); end

%% 形成矩阵Qk

a=zeros(1,n); %建立存放σ的矩阵a b=zeros(1,N); %建立存放β的矩阵b for k=1:n for i=k:N

a(k)=G(i,k)^2+a(k); end

a(k)=-sign(G(k,k))*(sqrt(a(k)));

w(k)=G(k,k)-a(k); %建立存放ω的矩阵w for j=k+1:1:N w(j)=G(j,k); end

b(k)=a(k)*w(k); %% 变换Gk-1到Gk G(k,k)=a(k);

for j=k+1:1:n+1 sum_wg=0; for i=k:N

sum_wg=w(i)*G(i,j)+sum_wg; end

t=sum_wg/b(k); for i=k:1:N

G(i,j)=t*w(i)+G(i,j); end end end

%% 解三角方程 Rx=h1 X(n)=G(n,n+1)/G(n,n); for i=(n-1):-1:1 sum_gx=0; for j=i+1:n

sum_gx=G(i,j)*X(j)+sum_gx; end

X(i)=(G(i,n+1)-sum_gx)/G(i,i); sum_gx=0; end

%% 参数X(i)的回代

p=zeros(1,N); %建立一维最小二乘近似数组p(x) for j=1:N for i=1:n

p(j)=p(j)+X(i)*x(j)^(i-1); end end

%% 显示拟合函数各阶系数 disp('各项拟合系数为：'); for i=1:n

disp([num2str(i-1),'次项系数',num2str(X(i))]); end

%% 误差E2求解 E2=0;%多项式计算误差

for i=n+1:N

E2=G(i,n+1)^2+E2; end E2=sqrt(E2); disp('误差大小为'); disp(E2);

%% 绘制曲线，显示结果 plot(x,y,'r*',x,p,'k-','LineWidth',1.5); xlabel('时刻/h'); ylabel('平均气温/°C');

title('最小二乘法拟合的气温变化曲线图') xlim([0 24]);

legend('数据点','拟合曲线')

三、线性方程组求解。

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

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

1、算法实现的思想及依据和算法实现的结构

高斯消去法主要分为两大步骤，即消去过程和回带过程。算法借鉴课本GAUSSPP算法。由于题目中给出的系数矩阵是对角占优的矩阵，因此可以不选主元直接进行高斯消去；另外在非压缩格式带状矩阵中，存在着大量的非零元素，非零元素的运算毫无意义并且占用大量机器资源和时间，因此对课本中给出的GAUSSPP算法进行优化。

对于n阶、上带宽q、下带宽p的带状矩阵，选取p与q较大者(赋值给c)，在第1行到第n-c行，第k列，只需要计算到

ak1,k到

akc,k，每一行也只需从

ai,k1ai,kq(i从k+1到k+q)；在第n-c+1到n行，第k列，计算

ak1,k到

an,k，每一

行只需要从

ai,k1计算到

ai,n(i从k+1到n)，这样可以减小运算量。

对于压缩带状矩阵，其消去过程和非压缩带状矩阵基本一致，不同之处在于：压缩格式忽略了零元素，因此需要建立压缩格式带状矩阵与非压缩格式带状矩阵的对应关系。主元素对应关系：B(k,p+1)=A(k,k)（B是压缩格式带状矩阵，A是非压缩格式带状矩阵），编写程序时需要根据此关系建立其元素间的对应关系。

2、程序运行结果对比及分析

图3.1求解dat61文件结果

图3.2 求解dat62文件结果

图3.3求解dat63文件结果

图3.4求解dat64文件结果

3、Matlab代码`

%% 改进前的程序 clc; clear

data_fname='dat64.dat';%文件名

file_id = fopen(data_fname, 'rb');%以二进制格式打开源文件

Id = fread(file_id, 1, 'uint32');%数据文件标示 id=dec2hex(Id);

ver = fread(file_id, 1, 'int32');%版本号 n = fread(file_id,1,'int32');%方程组的阶数 p = fread(file_id,1,'int32');%带状矩阵上带宽 q = fread(file_id,1,'int32');%带状矩阵下带宽

%% 确定数据格式

if strcmp(dec2hex(ver),'102') %比较字符串，确定数据格式 str=[ data_fname ' 为非压缩矩阵']; disp(str);

d = fread(file_id,n^2,'float'); %非压缩格式，需要读取n^2个浮点数，以一维格式存储到中间变量d

b = fread(file_id,n,'float'); %再读取右端向量的n个元素

%% 将读取到的数据放到阶数为n，上带宽为q，下带宽为p的系数矩阵A中 k = 1; for i = 1:n for j = 1:n

A(i,j) = d(k); k = k+1; end end

fclose(file_id); %% 消去过程 for k=1:n-1 for i=k+1:n

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

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

%% 回带过程求解方程组的根 x(n)=b(n)/A(n,n); for k=n-1:-1:1 S=b(k); for i=k+1:n

S=S-A(k,i)*x(i); end

x(k)=S/A(k,k); end end

%% 若为压缩矩阵

if strcmp(dec2hex(ver),'202')

str=[ data_fname ' 为压缩矩阵']; disp(str);

d = fread(file_id,n*(p+q+1),'float'); %压缩格式一共要读取n*(p+q+1)个数 b = fread(file_id,n,'float'); %再读取右端向量的n个元素 %% 将读取到的数据放到n行、p+q+1列的系数矩阵中 k=1;

for i = 1:n

for j = 1:(q+p+1) A(i,j) = d(k); k = k+1; end end

fclose(file_id); %% 消去过程 for k=1:n-q for i=1:p

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

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

for k = n-q+1:n-1 for i = 1:n-k

A(k+i,p+1-i)=A(k+i,p+1-i)/A(k,p+1); for j = 1: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)-A(k+i,p+1-i)*b(k); end end

%% 回带过程

x(n)=b(n)/A(n,p+1); for k=n-1:-1:n-q+1 S=b(k);

for i = k+1:n

S=S-A(k,p+1+i-k)*x(i); end

x(k)=S/A(k,p+1); end

for k=n-q:-1:1 S=b(k);

for j = k+1:k+q

S=S-A(k,p+1+j-k)*x(j); end

x(k)=S/A(k,p+1); end end

%% 部分解输出

disp('方程的前5个根:'); %输出5个根，用于与正确解对比 for j = 1:5

fprintf('%5.5f ',x (j)) %输出小数点后保留5位数的浮点数 end

fprintf('\\n');

4.本专业中的实际问题，提炼一个使用方程组进行求解的例子

机械系统中常用的组合弹簧问题就是典型的用到带状稀疏矩阵方程组进行求解的问题。

如图所示组合弹簧系统，已知：

刚度系数k1=100N/mm,k2=200N/mm,k3=100N/mm;

各点受力F1=-200N(以向右为正方向)，F2=100N，F3=500N,f4=300N. 求各节点处的位移u1,u2,u3,u4。

根据叠加原理，直接得组合弹簧系统的刚度矩阵： 1000K=⌈

00

阵。

−100300−2000

0−200300−100

00

⌉，显然，刚度矩阵为带状稀疏，对称矩−100100

根据胡克定律，有F=K*U,其中:

F=[F1,F2,F3,F4]T=[-200,100,500,300]T; U=[u1,u2,u3,u4]T

−2001001000

解方程组：[]=[

50003000−1000

300−200−2003000−1000𝑢1

0𝑢2

][]得各点位移： −100𝑢3100𝑢4

u1 =7.0000mm; u2=9.0000mm; u3=13.0000mm; u4=16.0000mm