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

1、实用标准计算方法（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请用合适的曲线拟合所测数据点；估算所需光缆

2、长度的近似值，并作出铺设河底光缆的曲线图;解：1、算法实现的思想及依据题目（1）为曲线拟合问题多项式插值、分段插值和最小二乘法。多项式插 值，随着插值数据点的数目增多，误差也会随之增大，因此不选用。最小二乘法 适于数据点较多的场合，在此也不适用。故选用分段插值。分段插值又分为分段线性插值、分段二次插值、三次样条插值及更高阶的多 项式插值。由本题的物理背景知，光缆正常工作时各点应该是平滑过渡， 因此至 少选用三次样条插值法。对于更高阶的多项式插值，由于“龙格现象”而不选用。题目（2）求光缆长度，即求拟合曲线在 0到20的长度，对弧长进行积分小f '(x)2dx。19 k i即可。光缆长度

3、的第一型线积分表达式为l文案大全实用标准2、算法实现的结构参考教材给出的SPLINEM算法和TTS算法，在选定边界条件和选定插值点等距分布后，可以先将数据点的二阶差商求出来并赋值给不端叫量d ,再根据TSS解法求解M。光缆长度的第一型线积分表达式为l :f'(x)2dx oKk 03、程序运行结果及分析0狂点；- o 111J?-jr 11 I? 则篇索一前一稗电霹聿启三茏弹条的底nr苜三二住年由式阳台图1.1三种边界条件下三次样条插值你想使用第几种边界条件Y请物入1、2、口之一：1LODCI后第1神边界条件下长度38米图1.2光缆长度4、MATLAB 代码:1)自己编程实现代码 cl

4、ear;文案大全实用标准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.229.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);fo门=2:1:21y1(i)=(y(i)-y(i-1)/(x(i)-x(i

5、-1);end%计算二阶差商y2=ones(1,21);for i=3:1:21y2(i)=(y1(i)-y1(i-1)/(x(i)-x(i-2);end%计算三阶差商y3=ones(1,21);for i=4:1:21y3(i)=(y2(i)-y2(i-1)/(x(i)-x(i-3);end文案大全实用标准%选择边界条件(I)if I=1d(1)=0;d(21)=0;a(21)=0;c(1)=0;%第一个点和最后一个点的二阶差商为endif I=2d(1)=6*y1(1);d(21)=-6*y1(21);a(1)=1;c(1)=1;endif I=3d(1)=-12*y3(1);d(21)=

6、-12*y3(21);a(21)=-2;c(1)=-2;%endfor i=2:20d(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:21l(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:

7、-1:1xg(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:Nfor 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)*f1A3/6/h+M(t)*f2A3/6/h+(y(t-1)-M(t-1)*hA2/6)*f1+(y(t)-M(t)*hA2/6)* f2)/h;Sp(i

8、)=-M(t-1)*f1A2/2/h+M(t)*f2A2/2/h+(y(t)-y(t-1)/h-(M(t)-M(t-1)*h/6;length(i+1)=sqrt(1+Sp(i)A2)*(20/(N-1)+length(i);%第一类线积分end figure(1); plot(x1,-S,'r-')% 深度曲线 grid disp('第',num2str(I),'种边界条件下长度,num2str(length(N+1),' 米') axis fill; xlabel('测点/米');ylabel('深度/米

9、9;)； title('三次样条曲线拟合,);文案大全实用标准legend（'数据点,拟合曲线,3）;二、最小二乘近似假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气 温变化的规律；试计算这一天的平均气温，并试估计误差。时刻0123456789101112平均气温11111111222225444456800358时刻111111122222345678901234平均气温333222222111141975420876解:1、算法实现的思想及依据本题共有25个数据点，数据量较大，因此选用多项式插值会导致较大的误 差，插值样条函数虽然有较好的数值性质，但是 M

10、i的计算量不小，而且没有统 一的公式来表示。另外，本题要求找出温度变化规律，插值方法要求p（x）必须通 已知数据点，只会导致拟合曲线失去本有的数据所表示的规律；从表中的数据点 可以看出，温度并不准确（温度只测量到整数位），因此没有必要要求拟合曲线 通过数据点。因此，选取最小二乘近似。2、算法实现的结构算法参考课本LSS算法。利用已有数据来生成了 G,将y作为其最后一列 （第n+1歹1）存放。先形成矩阵Qk,再变换Gk-1到Gk；求解三角方程Rx=h1 ； 最后根据定义计算误差。平均温度采用数据点相加求和求平均，避免数值积分的 繁琐。文案大全实用标准3、程序运行结果及分析区副科卜史2尸fl -悬

11、金幔«也曲古跳=J牵1呼工 口 葺 意 声Dr平均温里二为:21.200000道输入多项式近偎更I翻最高项的次梆U 2各项批经系数为：吗酸系敷乳和的1次项系赶之6Q642次项廉iftm;弱艮墓大小为433上片麻,sjikh，二K.1. n -(=, 7A+,正 46：jli101邑IX*H- Q p#lrwr 牛彝分行雷匚平均温度为：21. 2000(10得如工绛取，迂份函数既高底的次覆上：3各项樨合系数为，小对页系致1* 388I次血天数-U 22733；次项系触匚之。；的之次啖泵顿- OOS2SJ1误差走小为1LN在平均瑞度T为：21.200000诱输久多项式近设函孙最高项的次数

12、五：4 备臣职合系却为；(次顶系翻，工弱1次喷系数7.而32次咬系机L醺成：N次顼系数-L 03起受4次项系物。009"&3误差大小为7.48384、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:Nsum=sum+y(i);endaverage=sum/N;fprintf(' 平均

13、温度 T 为：%fn',average);h=input('请输入多项式近似函数最高项的次数h :');n=h+1;%参课本LSS算法G=zeros(N,n+1);%建立一个 N行，n+1列的矩阵 G%生成矢I阵G各个元素for i=1:n%矩阵G中前N列元素for j=1:NG(j,i)=x(1,jF(i-1);endendfor i=1:1:N%将丫作为矩阵 G中第(N+1)列元素G(i,n+1)=y(i);end%形成矢I阵Qka=zeros(1,n);%建立存放o-的矩阵ab=zeros(1,N); % 建立存放3的矩阵bfor k=1:nfor i=k:N文案

14、大全实用标准a(k)=G(i,k)A2+a(k);enda(k)=-sign(G(k,k)*(sqrt(a(k);w(k)=G(k,k)-a(k);%建立存放3的矩阵 wfor j=k+1:1:Nw(j)=G( j,k);endb(k)=a(k)*w(k);%变换Gk-1到GkG(k,k)=a(k);for j=k+1:1:n+1sum_wg=0;for i=k:Nsum_wg=w(i)*G(i,j)+sum_wg;endt=sum_wg/b(k);for i=k:1:NG(i,j)=t*w(i)+G(i,j);endendend%解三角方程 Rx=h1X(n)=G(n,n+1)/G(n,n)

15、;for i=(n-1):-1:1sum_gx=0;文案大全实用标准for j=i+1:nsum_gx=G(i,j)*X(j)+sum_gx;endX(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:Nfor i=1:nP(j)=P(j)+X(i)*x(j)A(i-1);endend%显示拟合函数各阶系数disp('各项拟合系数为：,);for i=1:ndisp(num2str(i-1),'次项系数',num2str(X(i);end%误差E

16、2求解E2=0;%多项式计算误差for i=n+1:NE2=G(i,n+1)A2+E2;endE2=sqrt(E2);disp('误差大小为,);文案大全实用标准disp（E2）;%绘制曲线，显示结果plot（x,y,'r*',x,p,'k-','LineWidth',1.5）;xlabel（'时刻/h'）;ylabel（'平均气温 / '）;title（'最小二乘法拟合的气温变化曲线图'）xlim（0 24）;legend（'数据点',拟合曲线）三、线性方程组求解。（1）编

17、写程序实现大规模方程组的高斯消去法程序，并对所附的方程组进行求 解。所附方程组的类型为对角占优的带状方程组。（2）针对本专业中所碰到的实际问题，提炼一个使用方程组进行求解的例子， 并对求解过程进行分析、求解。解：1、算法实现的思想及依据和算法实现的结构高斯消去法主要分为两大步骤，即消去过程和回带过程。算法借鉴课本 GAUSSPP算法。由于题目中给出的系数矩阵是对角占优的矩阵，因此可以不选 主元直接进行高斯消去；另外在非压缩格式带状矩阵中，存在着大量的非零元素， 非零元素的运算毫无意义并且占用大量机器资源和时间，因此对课本中给出的 GAUSSPP算法进行优化。对于n阶、上带宽q、下带宽p的带状矩

18、阵，选取p与q较大者（赋值给c）, 在第1行到第n-c行，第k列，只需要计算ak9到ac,k,每一行也只需从ai,k 1 到ai,k q（i从k+1到k+q）；在第n-c+1到n行，第k歹1，计算ak 1,k到a1k ,每文案大全实用标准一行只需要从ai,k 1计算到4n (i从k+1到n),这样可以减小运算量。对于压缩带状矩阵，其消去过程和非压缩带状矩阵基本一致，不同之处在于：压缩格式忽略了零元素，因此需要建立压缩格式带状矩阵与非压缩格式带状矩阵 的对应关系。主元素对应关系：B(k,p+1)=A(k,k)(B是压缩格式带状矩阵，A是非压缩格式带状矩阵)，编写程序时需要根据此关系建立其元素间的

19、对应关系。2、程序运行结果对比及分析命令行窗口dat 1. dat为非压缩矩阵方程的前5个根一1.00000 1. 00000 1. 00000 L 00000 L 00000 程序运行耗时/S：0. Ci058图3.1求解dat61文件结果命令行窗口dat62. dat为压缩审时方程的前5个根：1.00000 1. OODOO L OOOOC 1.00000 1.00000Syiste* time: 1B-Dsc-2O16 10:17:35程序运行耗时0.0058573 S图3.2求解dat62文件结果命令¥丁寮口dat63r dat为非王痂矩降方程的前5个相：1.00000 1.

20、 00000 1. oioco 1.00000 t. toooo程序运行耗时245. 1247 I图3.3求解dat63文件结果文案大全实用标准命令行蜀口dat64. dat为压缩后阵方程的前5个根；2.01600 2.01600 2.01600 2.01600 2. 01600System tiae.18-Dtc-2016 10:19 25程序运行耗时.4674 SA»图3.4求解dat64文件结果3、Matlab代码、%改进前的程序clc;cleardata_fname='dat64.dat'% 文件名file_id = fopen(data_fname, 

21、9;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),'

22、102')%比较字符串，确定数据格式str= data_fname ' 为非压缩矩阵'文案大全实用标准disp(str);d = fread(file_id,nA2,'float'); %非压缩格式，需要读取nA2个浮点数，以一维格式存储到中间变量db = fread(file_id,n,'float'); %再读取右端向量的n个元素%将读取到的数据放到阶数为n,上带宽为q,下带宽为p的系数矩阵A中k = 1;for i = 1:nforj = 1:nA(i,j) = d(k);k = k+1; end endfclose(file_id

23、);%消去过程for k=1:n-1for i=k+1:nA(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);endb(i)=b(i)-A(i,k)*b(k);end文案大全实用标准end%回带过程求解方程组的根x(n)=b(n)/A(n,n);for k=n-1:-1:1S=b(k);for i=k+1:nS=S-A(k,i)*x(i);endx(k)=S/A(k,k);end end%若为压缩矩阵if strcmp(dec2hex(ver),'202')str= data_fname ' 为压缩矩

24、阵'%压缩格式一共要读取%再读取右端向量的n个元素disp(str);d = fread(file_id,n*(p+q+1),'float');n*(p+q+1)个数b = fread(file_id,n,'float');%将读取到的数据放到n行、p+q+1列的系数矩阵中k=1;for i = 1:nfor j = 1:(q+p+1)文案大全实用标准A(i,j) = d(k);k = k+1;endendfclose(file_id);%消去过程for k=1:n-qfor i=1:pA(k+i,p+1-i)=A(k+i,p+1-i)/A(k,p+1);for j=1:qA(k+i,p+1-i+j)=A(k+i,p+1-i+j)-A(k+i,p+1-i)*A(k,p+1+j);endb(k+i)=b(k+i)-b(k)*A(k+i,p+1-i);endendfor k = n-q+1:n-1for i = 1:n-kA(k+i,p+1-i)=A(k+i,p+1-i)/A(k,p+1);for j = 1:n-kA(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);文案大全实用标准ende