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

计算方法（B）上机作业一、三次样条拟合某通信公司在一次施工中，需要在水面宽度为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.8110.8010.93(1)请用合适的曲线拟合所测数据点；(2)估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图；解：1、算法实现的思想及依据题目（1）为曲线拟合问题多项式插值、分段插值和最小二乘法。多项式插值，随着插值数据点的数目增多，误差也会随之增大，因此不选用。最小二乘法适于数据点较多的场合，在此也不适用。故选用分段插值。分段插值又分为分段线性插值、分段二次插值、三次样条插值及更高阶的多项式插值。由本题的物理背景知，光缆正常工作时各点应该是平滑过渡，因此至少选用三次样条插值法。对于更高阶的多项式插值，由于“龙格现象”而不选用。题目（2）求光缆长度，即求拟合曲线在0到20的长度，对弧长进行积分即可。光缆长度的第一型线积分表达式为。2、算法实现的结构参考教材给出的SPLINEM算法和TTS算法，在选定边界条件和选定插值点等距分布后，可以先将数据点的二阶差商求出来并赋值给右端向量d，再根据TSS解法求解M。光缆长度的第一型线积分表达式为。3、程序运行结果及分析图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;% 第一个点和最后一个点的二阶差商为0endif I=2 d(1)=6*y1(1); d(21)=-6*y1(21); a(1)=1;c(1)=1;endif I=3 d(1)=-12*y3(1); d(21)=-12*y3(21); a(21)=-2;c(1)=-2;%endfor 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);endM=xg;%所有点的二阶导数值% 求函数表达式并积分t=1;h=1;N=1000x1=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)*f13/6/h+M(t)*f23/6/h+(y(t-1)-M(t-1)*h2/6)*f1+(y(t)-M(t)*h2/6)*f2)/h; Sp(i)=-M(t-1)*f12/2/h+M(t)*f22/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);%第一类线积分endfigure(1);plot(x1,-S,r-)%深度曲线griddisp(第,num2str(I),种边界条件下长度,num2str(length(N+1),米)axis fill;xlabel(测点/米);ylabel(深度/米);title(三次样条曲线拟合);legend(数据点,拟合曲线,3);二、最小二乘近似假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温变化的规律;试计算这一天的平均气温,并试估计误差。时刻0123456789101112平均气温15141414141516182020232528时刻131415161718192021222324平均气温313431292725242220181716解：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; %求一天的平均温度averagefor i=1:N sum=sum+y(i);endaverage=sum/N;fprintf(平均温度T为：%fn,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); endendfor i=1:1:N %将y作为矩阵G中第(N+1)列元素 G(i,n+1)=y(i);end% 形成矩阵Qka=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=h1X(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); endend% 显示拟合函数各阶系数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; endE2=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列，只需要计算到，每一行也只需从到(i从k+1到k+q)；在第n-c+1到n行，第k列，计算到，每一行只需要从计算到(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;cleardata_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,n2,float); %非压缩格式，需要读取n2个浮点数，以一维格式存储到中间变量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); endend% 若为压缩矩阵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); endend% 部分解输出disp(方程的前5个根:); %输出5个根，用于与正确解对比for j = 1:5 fprintf(%5.5f ,x (j) %输出小数点后保留5位数的浮点数endfprintf(n);4.本专业中的实际问题，提炼一个使用方程组进行求解的例子机械系统中常用的组合弹簧问题就是典型的用到带状稀疏矩阵方程组进行求解的问题。如图所示组合弹簧系统，已知：刚度系数k1=100N/mm,k2=200N/mm,k3=100N/mm;各点受力F1=-200N(以向右为正方向)，F2=100N，F3=500N,f4=300N.求各节点处的位移u1,u2,u3,u4。根据叠加原理，直接得组合弹簧系统的刚度矩阵：K=100-100000300-20000-200300-10000-100100，显然，刚度矩阵为带状稀疏，对称矩阵。根据胡克定律，有F=K*U,其中:F=F1,F2,F3,F4T=-200,100,500,300T; U=u1,u2,u3,u4T 解方程组：-200100500300=100-100000300-20000-200300-10000-100100u1u2u3u4得各点位移：u1 =7.0000mm; u2=9.0000mm; u3=13.0000mm; u4=16.0000mm