2022年西安交大计算方法b大作业课件

《计算方法B》上机实验报告学院：机械工程学院班级:姓名：学号：2015年12月22日

1.计算以下和式：，要求：

(1)若保留11个有效数字，给出计算结果，并评价计算的算法；

(2)若要保留30个有效数字，则又将如何进行计算。实现思想：以上问题出现了近似数相减的问题，为了减小误差，可分别求得减数之和以及被减数之和，最后将两者相减。另外，减数与被减数求和均为同号计算，按照绝对值递增顺序相加可减小舍入误差。此题中对有效数字有要求，因而计算时首先需要根据有效数字位数计算得出迭代次数，以保证计算值的精度。源程序：m=input('输入有效数字个数m=');s0=1;s1=0;s2=0;n=0;%判断迭代次数whiles0>=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%分别求解各项并求和fork=n-1:-1:0a1=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;endS=vpa(s1-s2,m)实验结果：11位有效数字计算结果如图1所示；30为有效数字计算结果如图2所示。 图1.11位有效数字计算结果 图2.30为有效数字计算结果

某通信公司在一次施工中，需要在水面宽度为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)估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图；算法思想：由于题中所给点数为20，若采用高次多项式插值将产生很大的误差，所以拉格朗日或牛顿并不适用。题中光缆为柔性，可光滑铺设于水底，鉴于此特性，采用三次样条插值插值法较为合适。算法结构：三次样条算法结构见《计算方法教程》P110；光缆长度计算公式：源程序：clear;clc;x=0:20;y=[9.018.967.967.978.029.0510.1311.1812.2613.2813.3212.6111.2910.229.157.907.958.869.8110.8010.93];d=y;plot(x,y,'k.','markersize',15)holdon%%%计算差商fork=1:2fori=21:-1:(k+1)d(i)=(d(i)-d(i-1))/(x(i)-x(i-k));endend%%%设定d的边界条件fori=2:20d(i)=6*d(i+1);endd(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);%%%追fork=2:21l(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);fork=20:-1:1m(k)=(yy(k)-r(k)*m(k+1))/u(k);end%%%绘制曲线k=1;nn=100;xx=linspace(0,20,nn);l=0;forj=1:nnfori=2:20ifxx(j)<=x(i)k=i;break;elsek=i+1;endendh=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;sp(j)=-m(k-1)*(x(k)-xx(j))^2/(2*h)+m(k)*(xx(j)-x(k-1))^2/(2*h)+(y(k)-y(k-1))/h-(m(k)-m(k-1))*h/6;l(j+1)=(1+sp(j)^2)^0.5*(20/nn)+l(j);%求解光缆长度end%%%绘图plot(xx,s,'r-','linewidth',1.5)disp(['¹光缆长度为ª',num2str(l(nn+1)),'Ã×'])曲线图如图2-1所示，计算光缆长度如图2-2所示。图2-1光缆插值曲线图图2-1光缆计算长度显示3.假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温变化的规律;试计算这一天的平均气温,并试估计误差。时刻0123456789101112平均气温15141414141516182020232528时刻131415161718192021222324平均气温313431292725242220181716实现思想：此题中所给数据点数目较多，采用拉格朗日插值法或者牛顿插值法需要很高次的多项式，计算困难，误差大；采用样条插值计算量虽然不大，但是存放参数Mi的量很大，且没有一个统一的数学公式来表示，也不是很方便。所以可考虑用最小二乘法进行拟合。计算过程中，分别使用二次函数、三次函数以及四次函数，计算其相应的系数，估算误差并作图比较各个函数之间的区别。算法结构：（参考课本P123）1.1[形成矩阵Qk]1.2[变换Gk-1到Gk]2.[求解三角方程]3.[计算误差]源代码：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];

m=length(x);

n=input('请输入函数的次数

');

plot(x,y,'k.',x,y,'-')

grid;

hold

on;

n=n+1;

G=zeros(m,n+1);

G(:,n+1)=y';

c=zeros(1,n);%建立c来存放σ

q=0;

f=0;

b=zeros(1,m);%建立b用来存放β

形成矩阵G

j=1:n

i=1:m

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

建立矩阵Qk

k=1:n

i=k:m

c(k)=G(i,k)^2+c(k);

c(k)=-sign(G(k,k))*(c(k)^0.5);

w(k)=G(k,k)-c(k);%建立w来存放ω

j=k+1:m

w(j)=G(j,k);

b(k)=c(k)*w(k);

变换矩阵Gk-1到Gk

G(k,k)=c(k);

j=k+1:n+1

q=0;

i=k:m

q=w(i)*G(i,j)+q;

s=q/b(k);

i=k:m

G(i,j)=s*w(i)+G(i,j);

求解三角方程

a(n)=G(n,n+1)/G(n,n);

i=n-1:(-1):1

j=i+1:n

f=G(i,j)*a(j)+f;

a(i)=(G(i,n+1)-f)/G(i,i);

（i）存放各级系数

f=0;

a

回代过程

p=zeros(1,m);

j=1:m

i=1:n

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

plot(x,p,'r*',x,p,'-');

E2=0;%用E2来存放误差

误差求解

i=n+1:m

E2=G(i,n+1)^2+E2;

E2=E2^0.5;

disp('误差为');

disp(E2);

t=0;

i=1:m

t=t+p(i);

%%%平均温度

disp(['平均温度为',num2str(t),'℃'])实验结果：二次函数拟合，结果如下图所示图3-1二次函数拟合结果三次函数拟合，结果如下图所示图3-2三次函数拟合结果四次函数拟合，结果如下图所示图3-3四次函数拟合结果结果对比：将二次函数、三次函数和四次函数拟合结果绘制在同一个坐标内，如图3-4所示。其计算误差结果见表3-1所示。图3-4拟合结果对比分析

4.设计算法，求出非线性方程的所有实根，并使误差不超过。算法思想：本题可采用牛顿法迭代求解,令,得带格式为根据函数图像可以找出根的大致分布区间,带入不同的初值即可解出不同的根.源代码:functiony=f2(x)y=6*x.^5-45*x.^2+20;%定义原函数functiony=f3(x)y=30*x^4-90*x;%定义原函数倒数i=-5:0.1:5;y=f2(i);plot(i,y)holdonplot(i,0,'-')%画出原函数图像%%Newton法求根x1=input('输入初值');e=10^(-4);%误差设定Nmax=1000;%迭代最大次数限定forn=1:Nmaxf0=f2(x1);ifabs(f2(x1))<efprintf('输出的f(x)已经足够小');x=x1;breakelseF0=f3(x1);x=x1-f0/F0;ifabs(x-x1)<ebreakelsex1=x;endendendfprintf('输出方程的根x=%2f',x)计算结果:函数图像如图4-1所示。计算结果分别见图4-2所示。图4-1函数图像图4-2计算结果根据带入不同的初值，可以求出不同的根，有图4-2可以看出，原函数的根大约有三个，分别是-0.654542、0.681174、1.870799。

5.线性方程组求解。（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('文件标识符为');idver=dec2hex(head(2));%读取版本号fprintf('文件版本号为');vern=head(3);%读取阶数fprintf('矩阵A的阶数');nq=head(4);%上带宽fprintf('矩阵A的上带宽');qp=head(5);%下带宽fprintf('矩阵A的下带宽');pdist=4*num;fseek(file,dist,'bof');%把句柄值转向第六个元素开头处[A,count]=fread(file,inf,'float');%读取二进制文件，获取系数矩阵fclose(file);%关闭二进制头文件%%对非压缩带状矩阵进行求解ifver=='102',a=zeros(n,n);fori=1:n,forj=1:n,a(i,j)=A((i-1)*n+j);%求系数矩阵a(i,j)endendb=zeros(n,1);fori=1:n,b(i)=A(n*n+i);endfork=1:n-1,%列主元高斯消去法m=k;fori=k+1:n,%寻找主元ifabs(a(m,k))<abs(a(i,k))m=i;endendifa(m,k)==0%遇到条件终止disp('错误！')returnendforj=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;endfori=k+1:n,%计算l(i,k)并将其放到a(i,k)中a(i,k)=a(i,k)/a(k,k);forj=k+1:na(i,j)=a(i,j)-a(i,k)*a(k,j);endb(i)=b(i)-a(i,k)*b(k);endendx=zeros(n,1);%回代过程x(n)=b(n)/a(n,n);fork=n-1:-1:1,x(k)=(b(k)-sum(a(k,k+1:n)*x(k+1:n)))/a(k,k);endend%%对压缩带状矩阵进行求解ifver=='202',%高斯消去法m=p+q+1;a=zeros(n,m);fori=1:1:nforj=1:1:ma(i,j)=A((i-1)*m+j);endendb=zeros(n,1);fori=1:1:nb(i)=A(n*m+i);%求b(i)endfork=1:1:(n-1)%开始消去ifa(k,(p+1))==0disp('错误！');break;endst1=n;if(k+p)<nst1=k+p;endfori=(k+1):1:st1a(i,(k+p-i+1))=a(i,(k+p-i+1))/a(k,(p+1));forj=(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);endb(i)=b(i)-a(i,k+p-i+1)*b(k);endendx=zeros(n,1);%回代x(n)=b(n)/a(n,p+1);sum=0;fork=(n-1):-1:1sum=b(k);st2=n;if(k+q)<nst2=k+q;endforj=(k+1):1:st2sum=sum-a(k,j+p-k+1)*x(j);endx(k)=sum/a(k,p+1);sum=0;endenddisp('方程组的的解为：')%输出解disp(x’)求解结果对数据文件dat51求解，结果如下：文件标识符为id=F1E1D1A0文件版本号为ver