西安交通大学计算方法上机作业.docx
- 文档编号:26872864
- 上传时间:2023-06-23
- 格式:DOCX
- 页数:26
- 大小:168.43KB
西安交通大学计算方法上机作业.docx
《西安交通大学计算方法上机作业.docx》由会员分享,可在线阅读,更多相关《西安交通大学计算方法上机作业.docx(26页珍藏版)》请在冰豆网上搜索。
西安交通大学计算方法上机作业
计算方法上机作业
1.对以下和式计算:
,要求:
(1)若只需保留11个有效数字,该如何进行计算;
(2)若要保留30个有效数字,则又将如何进行计算;
(1)解题思想和算法实现:
根据保留有效位数的要求,可以由公式
得出计算精度要求。
只需要很少内存,时间复杂度和d呈线性,不需要高浮点支持。
先根据while语句求出符合精度要求的n值的大小,然后利用for语句对这n项进行求和,输出计算结果及n值大小即可。
(2)matlab源程序:
保留11位有效数字时;
clear
clc
formatlong
n=0;
sum=1/(16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6));
whilesum>=5*10^(-11);
n=n+1;
sum=1/(16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6));
end
fori=0:
n-1;
sum=sum+1/(16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6));
end
vpa(sum,11)
n
保留30位有效数字时;
clear
clc
formatlong
n=0;
sum=1/(16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6));
whilesum>=5*10^(-30);
n=n+1;
sum=1/(16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6));
end
fori=0:
n-1;
sum=sum+1/(16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6));
end
vpa(sum,30)
n
(3)实验结果分析
图1.1保留11位有效数字的n值及计算结果图图1.2保留30位有效数字的n值及计算结果图
由计算结果可知,通过合理的误差控制,分别通过7次和22次循环,可以实现题目所要求的精确度。
2.某通信公司在一次施工中,需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。
在铺设光缆之前需要对沟底的地形进行初步探测,从而估计所需光缆的长度,为工程预算提供依据。
已探测到一组等分点位置的深度数据(单位:
米)如下表所示:
分点
0
1
2
3
4
5
6
深度
9.01
8.96
7.96
7.97
8.02
9.05
10.13
分点
7
8
9
10
11
12
13
深度
11.18
12.26
13.28
13.32
12.61
11.29
10.22
分点
14
0
深度
9.15
7.90
7.95
8.86
9.81
10.80
10.93
(1)请用合适的曲线拟合所测数据点;
(2)预测所需光缆长度的近似值,并作出铺设河底光缆的曲线图;
(1)解题思想和算法原理
给定区间[a,b]一个分划
⊿:
a=x0 若函数S(x)满足下列条件: 1)S(x)在每个区间[xi,xj]上是不高于3次的多项式。 2)S(x)及其2阶导数在[a,b]上连续。 则称S(x)使关于分划⊿的三次样条函数。 (2)matlab源程序: clc,clear x=0: 1: 20; y=[9.018.967.967.978.029.0510.1311.1812.2613.2813.3212.6111.2910.229.157.97.958.869.8110.8010.93]; n=length(x); l (1)=0;m(n)=0; h=diff(x); df=diff(y)./diff(x); d (1)=0;d(n)=0; forj=2: n-1 l(j)=h(j)/(h(j-1)+h(j)); m(j)=h(j-1)/(h(j-1)+h(j)); d(j)=6*(df(j)-df(j-1))/(h(j-1)+h(j)); end m=m(2: end); u=diag(m,-1);r=diag(l,1);a=diag(2*ones(1,n)); A=u+r+a; M=inv(A)*d'; symsg forj=1: n-1 s(j)=M(j)*(x(j+1)-g)^3/(6*h(j))+M(j+1)*((g-x(j))^3/(6*h(j)))+(y(j)-M(j)*h(j)^2/6)*(x(j+1)-g)/h(j)+(y(j+1)-M(j+1)*h(j)^2/6)*(g-x(j))/h(j); end s r=0; forj=1: n-1 df=diff(s(j),g); warningoffall; q=int(sqrt(1+df.^2),g,j-1,j); r=r+q; end L=vpa(r,8); disp('thelengthofthelabelisL='); disp(L); forj=1: n-1 S(j,: )=sym2poly(s(j)); end forj=1: n-1 x1=x(j): 0.1: x(j+1); y1=polyval(S(j,: ),x1); ifj==1 y2=y1; else fori=1: 11 k=(j-1)*10+i; y2(k)=y1(i); end end end x2=x (1): 0.1: x(n); plot(x,y,'o') grid holdon plot(x2,y2,'r') (3)实验结果分析 图2.1铺设河底电缆长度 图2.2铺设河底光缆的曲线图 由三次样条插值得出的函数曲线的长度和即铺设河底电缆的长度为26.498514。 为了提高插值精度,用三次样条插值可以增加插值节点的办法来满足要求,而且在给定节点数的条件下,三次样条插值的精度要优于多项式插值以及线性分段插值,虽然舍弃了降低误差这个优点,但是其曲线的光滑性要好一些。 3.假定某天的气温变化记录如下表所示,试用数据拟合的方法找出这一天的气温变化的规律;试计算这一天的平均气温,并试估计误差。 时刻 0 11 12 平均气温 5 16 8 时刻 8 19 20 21 22 23 24 平均气温 3 24 22 20 18 17 16 (1)解题思想和数学原理: 对于具体实验时,通常不是先给出函数的解析式,再进行实验,而是通过实验的观察和测量给出离散的一些点,再来求出具体的函数解析式。 又因为测量误差的存在,实际真实的解析式曲线并不一定通过测量给出的所有点。 最小二乘法,形成法方程是求解这一问题的很好的方法,本实验运用这一方法实现对给定数据的拟合。 对于给定的测量数据(xi,fi)(i=1,2,…,n),设函数分布为 特别的,取 为多项式 (j=0,1,…,m) 则根据最小二乘法原理,可以构造泛函 令 (k=0,1,…,m) 则可以得到法方程 求该解方程组,则可以得到解 ,因此可得到数据的最小二乘解 (2)matlab源程序: x=[0123456789101112131415161718192021222324];%给出题目数据(时间) y=[15141414141516182020232528313231292725242220181716];%给出题目数据(温度) plot(x,y,'m*')%画出各个离散数据点 holdon forn=2: 4;%2、3、4代表拟合函数最高项次数 alltemp=25;%alltemp代表数据点总共有25个 A=zeros(n+1,n+1);%定义初始正规方程组的系数矩阵A C=ones(n+1,1);%定义初始正规方程组的系数向量C D=zeros(n+1,1);%定义初始正规方程组的向量D fori=1: n+1 forj=1: n+1 fork1=1: alltemp A(i,j)=A(i,j)+(x(k1).^(i-1+j-1)); end end fork2=1: alltemp D(i,1)=D(i,1)+(x(k2).^(i-1)).*(y(k2)); end end%以上为计算出正规方程组矩阵A、D的所有元素的程序 tol=1.0e-12; maxit=1000; C=bicg(A,D,tol,maxit);%使用bicg迭代算出正规方程组的系数向量C p=0;%误差分量 E=0;%误差总量 ifn==2 b=0: 24; f=C (1)+C (2).*b+C(3).*(b.^2); p=y(b+1)-f; forv=1: 25 E=E+(p(v)).^2; end plot(b,f,'r-') end%以上是对2阶拟合函数的图形处理与误差计算 ifn==3 b=0: 24; f=C (1)+C (2).*b+C(3).*(b.^2)+C(4).*(b.^3); p=y(b+1)-f; forv=1: 25 E=E+(p(v)).^2; end plot(b,f,'g-') end%以上是对3阶拟合函数的图形处理与误差计算 ifn==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; forv=1: 25 E=E+(p(v)).^2; end plot(b,f,'b-') end%以上是对4阶拟合函数的图形处理与误差计算 C,E end n=2;%重新对n赋值,进行指数函数拟合 A=zeros(n+1,n+1);%重新对A矩阵赋初值 C=zeros(n+1,1);%重新对C向量赋初值 D=zeros(n+1,1);%重新对D向量赋初值 fori=1: n+1 forj=1: n+1 fork=1: alltemp A(i,j)=A(i,j)+(x(k).^(i-1+j-1)); end end forl=1: alltemp D(i,1)=D(i,1)+((x(l).^(i-1)).*(log(y(l)))); end end%计算出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; forv=1: 25 E=E+(p(v)).^2; end plot(b,f,'c-')%对指数拟合函数进行图形处理和误差计算 b=-C(3); c=C (2)/(2*b); a=exp(b*(c^2)+C (1));%算出题设要求的指数拟合函数的各个系数 a,b,c,E gridon legend('测量数据','二次函数','三次函数','四次函数','指数拟合','Location','NorthWest') holdoff%holdon与holdoff联合使用,表示将各个曲线画在同一个图中 图3.1二次曲线拟合系数与2范数误差 图3.2三次曲线拟合系数与2范数误差 图3.3四次曲线拟合系数与2范数误差 图3.4指数曲线拟合系数与2范数误差 图3.5数据原始点与拟合曲线对比图 (3)实验结果分析: 根据国家有关规定,用的是02、08、14、20时四个观测时次的数据做平均,最有代表性。 从图中可以看出并不是多项式次数越高越好,随着次数的增高,曲线所呈现出的给定点处和实际的吻合度越好,但对于其他地方的吻合度降低了。 4.设计算法,求出非线性方程 的所有根,并使误差不超过 。 解: (1)解题思想和算法实现: 对于一个非线性方程的数值解法很多。 在此介绍两种最常见的方法: 二分法和Newton法。 首先要研究函数的形态,确定根的数量和大致区间的位置。 对于二分法,其数学实质就是说对于给定的待求解的方程f(x),其在[a,b]上连续,f(a)f(b)<0,且f(x)在[a,b]内仅有一个实根x*,取区间中点c,若,则c恰为其根,否则根据f(a)f(c)<0是否成立判断根在区间[a,c]和[c,b]中的哪一个,从而得出新区间,仍称为[a,b]。 重复运行计算,直至满足精度为止。 这就是二分法的计算思想。 Newton法通常预先要给出一个猜测初值x0,然后根据其迭代公式 产生逼近解x*的迭代数列{xk},这就是Newton法的思想。 当x0接近x*时收敛很快,但是当x0选择不好时,可能会发散,因此初值的选取很重要。 另外,若将该迭代公式改进为 其中r为要求的方程的根的重数,这就是改进的Newton法,当求解已知重数的方程的根时,在同种条件下其收敛速度要比Newton法快的多。 本题中用matlab画曲线 如下: 图4.1y=6*x.^5-45*x.^2+20的曲线 由曲线可以看出,该方程有三个实根,根所在区间依次为: [-1,-0.5][0.5,1][1.5,2] 采用Newton迭代法,依据迭代式: ,k=0,1,2,…(8-1) 该方程的迭代式为: (8-2) 分别选用迭代初值 、 、 进行迭代,分别求得满足精度的三个实根。 (2)matlab源程序: clc;clear x=-2: 0.1: 2; y=6*x.^5-45*x.^2+20; plot(x,y,'g')%画相应的曲线 grid title('y=6*x.^5-45*x.^2+20曲线') formatlong roots([600-45020])%根的真实解 clear x0=input('请输入迭代初值x0='); formatlong i=1; e=10^-4;%精度 x=x0-(6*x0.^5-45*x0.^2+20)/(30*x0.^4-90*x0); whileabs(x-x0)>e i=i+1; x0=x; x=x0-(6*x0.^5-45*x0.^2+20)/(30*x0.^4-90*x0);%迭代 end display('方程的根') x display('迭代的次数') i (3)实验结果分析: 图4.2运行结果 对于Newton迭代法,三个初值x0都使得迭代收敛,这是非常重要的。 考虑Newton法迭代的收敛性条件: (1)存在一个区间,满足。 由曲线和所选的三个区间可知这一条件满足。 (2)f(x)是[a,b]上的单调函数,即对一切不变号。 经验证所选的三个区间满足这一条件。 (3)f(x)的凹向在[a,b]上不变,即在[a,b]上不改变符号。 经验证所选的三个区间满足这一条件。 (4) -1 1 1.5 >0 >0 >0 另外,可以看出Newton迭代法收敛速度也很快,且很快达到很高的精度,源于它一般是超线性收敛的。 5.编写程序实现大规模方程组的列主元高斯消去法程序,并对所附的方程组进行求解。 针对本专业中所碰到的实际问题,提炼一个使用方程组进行求解的例子,并对求解过程进行分析、求解。 解: (1)算法原理 由于一般线性方程在使用Gauss消去法求解时,从求解的过程中可以看到,若 =0,则必须进行行交换,才能使消去过程进行下去。 有的时候即使 0,但是其绝对值非常小,由于机器舍入误差的影响,消去过程也会出现不稳定得现象,导致结果不正确。 因此有必要进行列主元技术,以最大可能的消除这种现象。 这一技术要寻找行r,使得 并将第r行和第k行的元素进行交换,以使得当前的 的数值比0要大的多。 这种列主元的消去法的主要步骤如下: 1.消元过程 对k=1,2,…,n-1,进行如下步骤。 1)选主元,记 很小,这说明方程的系数矩阵严重病态,给出警告,提示结果可能不对。 2)交换增广阵A的r,k两行的元素。 (j=k,…,n+1) 3)计算消元 (i=k+1,…,n;j=k+1,……,n+1) 2.回代过程 对k=n,n-1,…,1,进行如下计算 至此,完成了整个方程组的求解。 (2)matlab源程序: %非压缩,dat51.dat、dat53.dat clear;clc fp=fopen('dat53.dat','rb'); id=fread(fp,1,'int32'); ver=fread(fp,1,'int32'); N=fread(fp,1,'int32'); q=fread(fp,1,'int32'); p=fread(fp,1,'int32'); fori=1: N A(i,: )=fread(fp,N,'float'); end fori=1: N d(i)=fread(fp,1,'float'); end %正向消去过程 fori=1: N-q fork=1: p ll=A(i+k,i)/A(i,i); forj=i: i+q A(i+k,j)=A(i+k,j)-ll*A(i,j); end d(i+k)=d(i+k)-ll*d(i); end end fori=N-q+1: N fork=1: N-i ll=A(i+k,i)/A(i,i); forj=i: N A(i+k,j)=A(i+k,j)-ll*A(i,j); end d(i+k)=d(i+k)-ll*d(i); end end %回代过程 x(N)=d(N)/A(N,N); fori=N-1: -1: 1 S=0; ifi+q>N cv=N;%cv-criticalvalue elsecv=i+q; end forj=i+1: cv S=S+A(i,j)*x(j); end x(i)=(d(i)-S)/A(i,i); end x %压缩,dat52.dat、dat54.dat clear;clc fp=fopen('dat54.dat','rb'); id=fread(fp,1,'int32'); ver=fread(fp,1,'int32'); N=fread(fp,1,'int32'); q=fread(fp,1,'int32'); p=fread(fp,1,'int32'); fori=1: N A(i,: )=fread(fp,p+q+1,'float'); end fori=1: N d(i)=fread(fp,1,'float'); end %正向消去过程 fori=1: N ifi+p cv=p;%cv-criticalvalue else cv=N-i; end fork=1: cv ll=A(i+k,p+1-k)/A(i,p+1); forj=p+1: p+q+1 A(i+k,j-k)=A(i+k,j-k)-ll*A(i,j); end d(i+k)=d(i+k)-ll*d(i); end end %回代过程 x(N)=d(N)/A(N,p+1); fori=N-1: -1: 1; S=0; ii=i+1; jj=p+2; while(ii<=N)&&(jj<=p+q+1) S=S+A(i,jj)*x(ii); ii=ii+1; jj=jj+1; end x(i)=(d(i)-S)/A(i,p+1); end x (3)实验结果: 非压缩矩阵求解结果(部分) 压缩矩阵求解结果(部分) (4)分析心得: 采用Gauss消去法时,如果在消元时对角线上的元素始终较大,那么本方法不需要进行列主元计算,计算结果一般就可以达到要求,否则必须进行列主元这一步,以减少机器误差带来的影响,使方法得出的结果正确。 (5)实例 化学反应方程式严格遵守质量守恒定律,书写化学反应方程式写出反应物和生成物后,往往左右两边各原子数目不相等,不满足质量守恒定律,这就需要通过计算配平来解决。 对于化学反应方程式X1KMnO4+x2MnSO4+x3H2O=x4MnO2+x5K2SO4+x6H2SO4,可按以下方式配平: 上述化学反应式中包含5种不同的原子(钾、锰、氧、硫、氢),按照原子守恒有: K: x1=2x5 Mn: x1+x2=x4 O: 4x1+4x2+x3=2x4+4x5+4x6 S: x2=x5+x6 H: 3x2=2x6 上述方程组中有5个方程,6个未知量,因此有无数解。 可先令x6=1,于是形成方程组。 使用Gauss消去法求解此方程组得: x=[1.0000,1.5000,1.0000,2.5000,0.5000,1.0000];由于化学方程式通常取最简的正整数,因此将x乘2得配平的化学方程式: 2KMnO4+3MnSO4+2H2O=5MnO2+K2SO4+2H2SO4 程序如下: clear;clc; a=[1000-20; 110-100; 441-2-4-4; 0100-1-1; 00200-2; 00000-1]; b=[00000-1]'; n=length(b); fork=1: n-1 ifa(k,k)==0 disp('error'); return; end fori=k+1: n a(i,k)=a(i,k)/a(k,k); forj=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); fork=n-1: -1: 1 S=b(k); forj=k+1: n S=S-a(k,j)*x(j); end x(k)=S/a(k,k); end x
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 西安交通大学 计算方法 上机 作业
![提示](https://static.bdocx.com/images/bang_tan.gif)