计算方法上机作业

1、计算方法上机报告 上机实习题目1.某通信公司在一次施工中，需要在水面宽度为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)估算所需光缆长度的

2、近似值，并作出铺设河底光缆的曲线图；（1）算法思想分段多项式是由一些在相互连接的区间上的不同多项式连接而成的一条连续曲线，其中三次样条插值方法是一种具有较好“光滑性”的分段插值方法。在本题中，假设所铺设的光缆足够柔软，在铺设过程中光缆触地走势光滑，紧贴地面，并且忽略水流对光缆的冲击。计算光缆长度时，用如下公式： 本题采取三次样条插值的方法，因为三次样条插值方法是一种具有较好“光滑性”的分段插值方法。根据提供的数据，只用x,y值，不包含导数值，因此采用第三类三次插值多项式进行插值编程。设计算法如下：1.For 1.1 2.For 2.1 For 2.1.1 3.4.For 4.1 4.2 4.3

3、 5.6.7.获取M的矩阵元素个数，存入m8.For 8.1 8.2 8.3 9.10.For 10.1 11.获取x的元素个数存入s12.13.For 13.1 if then ；break else 14.（3） 源程序clear;clc; x=0:1:20; %产生从0到20含21个等分点的数组X=0:0.2: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.80,10.93; %等分点位置的深度数据n=length(

4、x); %等分点的数目N=length(X);% 求三次样条插值函数s(x) M=y; for k=2:3; %计算二阶差商并存放在M中 for i=n:-1:k; M(i)=(M(i)-M(i-1)/(x(i)-x(i-k+1); endendh(1)=x(2)-x(1); %计算三对角阵系数a,b,c及右端向量dfor i=2:n-1; h(i)=x(i+1)-x(i); c(i)=h(i)/(h(i)+h(i-1); a(i)=1-c(i); b(i)=2; d(i)=6*M(i+1);end M(1)=0; %选择自然边界条件M(n)=0;b(1)=2;b(n)=2;c(1)=0;a(

5、n)=0; d(1)=0; d(n)=0; u(1)=b(1); %对三对角阵进行LU分解y1(1)=d(1);for k=2:n; l(k)=a(k)/u(k-1); u(k)=b(k)-l(k)*c(k-1); y1(k)=d(k)-l(k)*y1(k-1);endM(n)=y1(n)/u(n); %追赶法求解样条参数M(i)for k=n-1:-1:1; M(k)=(y1(k)-c(k)*M(k+1)/u(k);ends=zeros(1,N);for m=1:N; k=1; for i=2:n-1 if X(m)<=x(i); k=i-1; break; else k=i; end

6、end H=x(k+1)-x(k); %在各区间用三次样条插值函数计算X点处的值 x1=x(k+1)-X(m); x2=X(m)-x(k); s(m)=(M(k)*(x13)/6+M(k+1)*(x23)/6+(y(k)-(M(k)*(H2)/6)*x1+(y(k+1)-(M(k+1)*(H2)/6)*x2)/H;end% 计算所需光缆长度L=0; %计算所需光缆长度for i=2:N L=L+sqrt(X(i)-X(i-1)2+(s(i)-s(i-1)2);enddisp('所需光缆长度为 L=');disp(L);figureplot(x,y,'*',X,s

7、,'-') %绘制铺设河底光缆的曲线图xlabel('位置X/测量点','fontsize',16); %标注坐标轴含义ylabel('深度Y/m','fontsize',16);title('铺设河底光缆的曲线图','fontsize',16);grid;（4）结果与分析 铺设海底光缆的曲线图如下图所示：拟合结果表明，运用分段三次样条插值所得的拟合曲线能较准确地反映铺设光缆的走势图。计算出所需光缆的长度为 L=26.4844m。可以用Newton法进行拟合，求得拟合曲线和光缆长度，

8、也可以用三次样条法来拟合，则精度会更高，根据实际光缆的现实特性，三次样条法来拟合。2.假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温变化的规律;试计算这一天的平均气温,并试估计误差。时刻0123456789101112平均气温15141414141516182020232528时刻131415161718192021222324平均气温313431292725242220181716（1）解题思想和数学原理： 对于具体实验时，通常不是先给出函数的解析式，再进行实验，而是通过实验的观察和测量给出离散的一些点，再来求出具体的函数解析式。又因为测量误差的存在，实际真实

9、的解析式曲线并不一定通过测量给出的所有点。最小二乘法，形成法方程是求解这一问题的很好的方法，本实验运用这一方法实现对给定数据的拟合。（2）matlab源程序：x=0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24; %给出题目数据（时间）y=15 14 14 14 14 15 16 18 20 20 23 25 28 31 32 31 29 27 25 24 22 20 18 17 16; %给出题目数据（温度）plot(x,y, 'm*') %画出各个离散数据点hold onfor n=2:4; %

10、2、3、4代表拟合函数最高项次数alltemp=25; % alltemp代表数据点总共有25个A=zeros(n+1,n+1); %定义初始正规方程组的系数矩阵AC=ones(n+1,1); %定义初始正规方程组的系数向量CD=zeros(n+1,1); %定义初始正规方程组的向量Dfor i=1:n+1 for j=1:n+1 for k1=1: alltemp A(i,j)=A(i,j)+(x(k1).(i-1+j-1); end end for k2=1: alltemp D(i,1)=D(i,1)+(x(k2).(i-1).*(y(k2); endend %以上为计算出正规方程组矩阵

11、A、D的所有元素的程序tol=1.0e-12;maxit=1000;C=bicg(A,D,tol,maxit); %使用bicg迭代算出正规方程组的系数向量Cp=0; %误差分量E=0; %误差总量if n=2 b=0:24; f=C(1)+C(2).*b+C(3).*(b.2); p=y(b+1)-f; for v=1:25 E=E+(p(v).2; end plot(b,f, 'r-')end %以上是对2阶拟合函数的图形处理与误差计算if n=3 b=0:24; f=C(1)+C(2).*b+C(3).*(b.2)+C(4).*(b.3); p=y(b+1)-f; for

12、 v=1:25 E=E+(p(v).2; end plot(b,f, 'g-')end %以上是对3阶拟合函数的图形处理与误差计算if n=4 b=0:24; f=C(1)+C(2).*b+C(3).*(b.2)+C(4).*(b.3)+C(5).*(b.4); p=y(b+1)-f; for v=1:25 E=E+(p(v).2; end plot(b,f, 'b-')end %以上是对4阶拟合函数的图形处理与误差计算C,Eendn=2; %重新对n赋值，进行指数函数拟合A=zeros(n+1,n+1); %重新对A矩阵赋初值C=zeros(n+1,1); %

13、重新对C向量赋初值D=zeros(n+1,1); %重新对D向量赋初值for i=1:n+1 for j=1:n+1 for k=1: alltemp A(i,j)=A(i,j)+(x(k).(i-1+j-1); end end for l=1: alltemp D(i,1)=D(i,1)+(x(l).(i-1).*(log(y(l); endend %计算出A矩阵、D向量各元素数值C=bicg(A,D,tol,maxit); %利用bicg迭代求解系数b=0:24;p=0;E=0;f=exp(C(1)+C(2).*b+C(3).*(b.2); p=y(b+1)-f; for v=1:25 E

14、=E+(p(v).2; endplot(b,f, 'c-') %对指数拟合函数进行图形处理和误差计算b=-C(3);c=C(2)/(2*b);a=exp(b*(c2)+C(1); %算出题设要求的指数拟合函数的各个系数a,b,c,Egrid onlegend('测量数据','二次函数','三次函数','四次函数','指数拟合','Location','NorthWest')hold off %hold on与hold off联合使用，表示将各个曲线画在同一个图中二次曲

15、线拟合系数与2范数误差 C = 8.4274 2.5606 -0.0920 E = 253.6476三次曲线拟合系数与2范数误差 C = 13.4072 -0.2164 0.2032 -0.0082 E = 110.2954四次曲线拟合系数与2范数误差 C = 16.6766 -3.5547 0.8592 -0.0513 0.0009 E = 43.9298 a = 26.0247 b = 0.0044 c = 14.0969 E = 191.67713.分析总结：通过上述几种拟合可以发现，多项式的次数越高，计算拟合的效果越好，误差越小，说明结果越准确；同时，指数多项式拟合的次数虽然不高，但误

16、差最小，说明结果最准确。 3.线性方程组求解。（1）编写程序实现大规模方程组的高斯消去法程序，并对所附的方程组进行求解。所附方程组的类型为对角占优的带状方程组。（1）算法思想 高斯消去法是利用现行方程组初等变换中的一种变换，即用一个不为零的数乘一个方程后加只另一个方程，使方程组变成同解的上三角方程组，然后再自下而上对上三角方程组求解。列主元消去法是当高斯消元到第步时，从列的以下（包括）的各元素中选出绝对值最大的，然后通过行交换将其交换到的位置上。交换系数矩阵中的两行（包括常数项），只相当于两个方程的位置交换了，因此，列选主元不影响求解的结果。 程序的核心就是高斯列主元消去法。根据教材提供的算法

17、，编写列主元消去法的子函数与适应于超大规模超出系统内存的方程组的改编程序。同时，在Gauss消去过程中，适当交换方程的顺序对保证消去过程能顺利进行及计算解的精确度都是有必要的，交换方程的原则是使中，绝对值最大的一个换到（k,k)位置而成为第k步消去的主元，这就是列主元Gauss消去法。（2）源程序clear; %清除工作空间变量clc; %清除命令窗口命令 % 读取系数矩阵f,p=uigetfile('*.dat','选择数据文件'); %读取数据文件num=5; %输入系数矩阵文件头的个数name=strcat(p,f);file=fopen(name,

18、9;r');head=fread(file,num,'uint'); %读取二进制头文件id=dec2hex(head(1); %读取标识符fprintf('文件标识符为');idver=dec2hex(head(2); %读取版本号fprintf('文件版本号为');vern=head(3); %读取阶数fprintf('矩阵A的阶数');nq=head(4); %上带宽fprintf('矩阵A的上带宽');q p=head(5); %下带宽fprintf('矩阵A的下带宽');p dis

19、t=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

20、-1, %列主元高斯消去法 m=k; for i=k+1:n, %寻找主元 if abs(a(m,k)<abs(a(i,k) m=i; end end 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)

21、-a(i,k)*a(k,j); end 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); endend% 对压缩带状矩阵进行求解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); %求a(i,j) end end b=zeros(n,1); for i=

22、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)<n st1=k+p; end 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);

23、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); st2=n; if (k+q)<n st2=k+q; end 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; endenddisp('方程组的的解为：') %输出方程组的解disp(x)1.文件dat61.dat 读取 2.文件dat62.dat 读取3.文件dat63.dat 读取4. 文件dat64.dat

24、 读取分析文件dat61.dat所得的结果是15维列向量，所得的根均是1.0000。dat62.dat所得的结果是20阶，上带宽是4，下带宽也是4，所得的根是 1.0000。dat63.dat所得的结果是2480阶，上带宽是4，下带宽也是4，所得的根均是1.0000，此计算结果耗时间很长，大约为400秒。dat64.dat所得的结果是46848阶，上带宽是5，下带宽也是5，所得的根均是2.0160，此计算结果耗时间为10秒。（3）分析心得：采用Gauss消去法时，如果在消元时对角线上的元素始终较大，那么本方法不需要进行列主元计算，计算结果一般就可以达到要求，否则必须进行列主元这一步，以减少机器

25、误差带来的影响，使方法得出的结果正确。（2） 针对本专业中所碰到的实际问题，提炼一个使用方程组进行求解的例子，并对求解过程进行分析、求解。我涉及船舶方向，需要对海域的进行测绘。在某海域测得一些点(x,y)处的水深z由下表给出，船的吃水深度为5英尺，在矩形区域（75，200）×（-50，150）里的哪些地方船要避免进入。xyz129 140 103.5 88 185.5 195 1057.5 141.5 23 147 22.5 137.5 85.5 4 8 6 8 6 8 8xyz 157.5 107.5 77 81 162 162 117.5 -6.5 -81 3 56.5 -66.

26、5 84 -33.5 9 9 8 8 9 4 9（1） 输入插值基点数据；（2） 在矩形区域（75，200）×（-50，150）作二维插值；（3） 作海底曲面图；作出水深小于5的海域范围，即z=5的等高线程序：%输入插值基点数据x=129 140 103.5 88 185.5 195 105 157.5 107.5 77 81 162 162 117.5;y=7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5;z=4 8 6 8 6 8 8 9 9 8 8 9 4 9;z=-z;%在矩形区域（75，200）×（-50，150）作二维插值cx=75:0.5:200;cy=-50:0.5:150;cz=griddata(x,y,z,cx,cy','cubic');%作海底曲面图subplot(1,2,1),meshz(cx,cy,cz)xlabel('x'),ylabel('y'),zla