西安交通大学计算方法B大作业资料docWord格式文档下载.docx
- 文档编号:20992616
- 上传时间:2023-01-26
- 格式:DOCX
- 页数:20
- 大小:168.61KB
西安交通大学计算方法B大作业资料docWord格式文档下载.docx
《西安交通大学计算方法B大作业资料docWord格式文档下载.docx》由会员分享,可在线阅读,更多相关《西安交通大学计算方法B大作业资料docWord格式文档下载.docx(20页珍藏版)》请在冰豆网上搜索。
t=0
n值加至n7
s=3.1415926536
②t=0
n值加至n22
s=3.14159265358979311599796346854
从上述的算法思想中可以看出,运算中不仅要满足误差要求,还要尽可能地减少计算量,此外还要考虑舍入误差的影响,这时就要对所运算数据的性质进行分析,设置合适的算法,从而提高运算的精度。
而逆序算法能够很好的满足上述要求。
题目二
2.1题目内容
某通信公司在一次施工中,需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。
在铺设光缆之前需要对沟底的地形进行初步探测,从而估计所需光缆的长度,为工程预算提供依据。
已探测到一组等分点位置的深度数据(单位:
米)如下表所示:
分点
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
15
16
17
18
19
20
9.15
7.90
7.95
8.86
9.81
10.80
10.93
(1)请用合适的曲线拟合所测数据点;
(2)估算所需光缆长度的近似值,并作出铺设河底光缆的曲线图;
2.2算法思想
利用曲线拟合数据点,即利用数据点拟合差值多项式,我们可以利用Newton法进行拟合,也可以用复化Simpson求积公式、三次样条插值来拟合,但三次样条插值使用方程组计算增大了计算量,同时还要附加边界条件,分段三次样条插值对图形的控制能力还不够灵活。
因此这里用Newton形式的差值多项式进行拟合。
首先计算出各差商,然后计算出Newton差值多项式的每一项,最后将所有项相加,即可计算出Newton差值多项式,然后利用所得的差值多项式一次算出多个点的函数值。
MATLAB的plot函数进行绘图。
计算长度近似值,只需将每隔两点之间的距离算出,然后一次相加,所得的折线长度,即为长度的近似值。
2.3Matlab源程序
Untitled2
clear
clc
x=0:
1:
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];
%输入给定的数据点
xi=0:
[Nx,Ni]=Newton(x,y,xi);
%调用函数,建立Newton差值多项式
plot(xi,Ni);
%绘制拟合的曲线图
long=0;
%为长度赋初值
fori=1:
20%将每一段折线相加算出长度的近似值
long=long+sqrt(1+((y(i)-y(i+1))^2));
disp('
需要的光缆长度为'
)%显示需要的光缆长度
disp(long)
Newton插值法
function[Nx,N0]=Newton(X,Y,x0)
n=size(X);
%插值点个数
y=Y;
Nx=Y
(1);
N=1;
n-1%计算Newton插值多项式
forj=i+1:
n
yi(j)=(y(j)-y(i))/(X(j)-X(i));
end
m(i)=yi(i+1);
N=N*(x-X(i));
Nx=Nx+N*m(i);
y=yi;
N0=subs(Nx,'
x'
x0);
2.4计算结果及总结
针对上述Matlab程序,铺设海底光缆的曲线图如下图所示:
结果如下:
Nl=26.4844,即为所求近似计算光缆长度。
本题利用Newton法进行拟合,既简单又使用,运行也可以用复化Simpson求积公式、三次样条插值来拟合,但三次样条插值使用方程组计算增大了计算量,同时还要附加边界条件,分段三次样条插值对图形的控制能力还不够灵活。
题目三
3.1题目内容
假定某天的气温变化记录如下表所示,试用数据拟合的方法找出这一天的气温变化的规律;
试计算这一天的平均气温,并试估计误差。
时刻
平均气温
23
25
28
21
22
24
31
34
29
27
3.2算法思想
在本题中,数据点的数目较多。
当数据点的数目很多时,用“多项式插值”方法做数据近似要用较高次的多项式,这不仅给计算带来困难,更主要的缺点是误差很大。
用“插值样条函数”做数据近似,虽然有很好的数值性质,且计算量也不大,但存放参数的量很大,且没有一个统一的数学公式来表示,也带来了一些不便。
另一方面,在有的实际问题中,用插值方法并不合适。
当数据点的数目很大时,要求通过所有数据点,可能会失去原数据所表示的规律。
如果数据点是由测量而来的,必然带有误差,插值法要求准确通过这些不准确的数据点是不合适的。
在这种情况下,不用插值标准而用其他近似标准更加合理。
通常情况下,是选取使最小,这就是最小二乘近似问题。
在本题中,采用“最小二乘法”找出这一天的气温变化的规律,使用二次函数、三次函数、四次函数以及指数型函数,计算相应的系数,估算误差,并作图比较各种函数之间的区别。
给定数据点和一组函数。
求数(假定),使函数
满足
为计算方便,只需考虑
最小二乘近似问题的解应使上式成立;
在
中,它是ak的二次函数,因此其取极小值的必要条件为
,即
经过推导知,若记
则最小二乘问题的法方程为
解此法方程,
就可以得到各项系数,数据拟合之后便可知道一天中温度变化趋势。
3.3Matlab源程序
A=0:
24;
B=[15,14,14,14,14,15,16,18,20,20,23,25,28,31,32,41,29,27,25,24,22,20,18,17,16];
%输入给定的数据
G=zeros(25,7);
%给矩阵G赋初值
fori=1:
25;
forj=1:
7;
G(i,j)=A(i)^(j-1);
x1=(G'
*G)\(G'
*B'
%按法方程法算出各项系数
a=zeros(1,7);
forn=1:
a(1,n)=x1(8-n,1);
end%各项系倒序排列
x=0:
0.1:
plot(x,polyval(a,x),'
b'
)%绘制曲线
holdon
fx=polyval(a,x)
sum=0
240%计算曲线下面积
sum=sum+(fx(i)+fx(i+1))*0.1/2
average=sum/24%计算平均温度
e=0
10:
241%计算误差
e=e+(fx(i)-B((i+9)/10)).^2
disp('
平均温度为'
)%输出平均温度
disp(average)
误差为'
)%输出误差值
disp(e)
3.4计算结果及总结
铺设海底光缆的曲线图如下图所示:
这一天平均温度为average=21.7729
估算误差e=118.7770
本算法利用最小二乘法,这里用法方程法进行求解,相比于正交化方法,此方法运算量较大,但对于本题来说,数据量并不大,因此在电脑上运行也不会很慢,这里生成的矩阵G由7列,因此有六阶算法的精度。
题目四
4.1题目内容
设计算法,求出非线性方程
的所有实根,并使误差不超过
。
4.2算法思想
采用Newton法进行求解,牛顿法的迭代函数
该式的导数
则迭代函数
迭代的初值为-1,0.5,1.5
4.3Matlab源程序
主程序
x1=-1;
%初值为-1时
whileabs(6*x1^5-45*x1^2+20)>
0.0001%误差大于0.0001时继续循环
x1=x1-(6*x1^5-45*x1^2+20)/(30*x1^4-90*x1);
%Newton迭代法
x2=0.5;
%初值为0.5时
whileabs(6*x2^5-45*x2^2+20)>
0.0001
x2=x2-(6*x2^5-45*x2^2+20)/(30*x2^4-90*x2);
x3=1.5;
%初值为1.5时
whileabs(6*x3^5-45*x3^2+20)>
0.0001
x3=x3-(6*x3^5-45*x3^2+20)/(30*x3^4-90*x3);
x1=%.6f,x2=%.6f,x3=%.6f'
x1,x2,x3)%输出三个根
绘图程序
x=-2:
y=6*x.^5-45*x.^2+20
plot(x,y)
xlabel('
X轴'
ylabel('
Y轴'
%坐标轴表示对象标签
gridon;
%显示网格线
axison;
%显示坐标轴
4.4计算结果及总结
x1=-0.654542x2=0.681174x3=1.870799
所得曲线
所给的方程形式较为简单,导数较为好求所以用牛顿法非常方便,先作图知道根的大致范围后再取初值可以有效的减小迭代次数。
题目五
5.1题目内容
线性方程组求解。
(1)编写程序实现大规模方程组的高斯消去法程序,并对所附的方程组进行求解。
所附方程组的类型为对角占优的带状方程组。
(2)针对本专业中所碰到的实际问题,提炼一个使用方程组进行求解的例子,并对求解过程进行分析、求解。
5.2算法思想
消去法的中心是“降维”,即求解n元方程组的问题转化为先解n-1元方程组,一旦这个n-1元方程组的解取得,则剩余的一个未知量自然可以求得。
这样逐步减少未知量个数的方法,是求解多元方程组的一个重要思想。
Gauss消去过程中,适当交换方程的顺序对保证消去过程能够顺利进行及计算解的精确度都是必要的。
消去过程中产生的数称为第k步消去的主元。
交换方程的原则是使中,绝对值最大的一个换到(k,k)位置而成为第k步消去的主元,带有这种交换的Gauss消去法为列主元Gauss消去法。
5.3Matlab源程序
5.3.1非压缩带状对角方程组
(1)dat51.dat、dat52.dat主程序:
%非压缩格式dat51.dat;
dat52.dat
formatshort
%读出数据
fid=fopen('
E:
\2015上机题目\实验五\dat52.dat'
'
r'
D=fread(fid,3,'
long'
0);
n=D(3)%矩阵维数
offset=20;
fseek(fid,offset,'
bof'
A=zeros(n);
B=zeros(1,n);
C=fread(fid,n*(n+1),'
float'
0)
forj=1:
A(i,j)=C((i-1)*n+j);
A;
B(i)=C(n*n+i);
B;
x=GAUSSPP(A,B)
(2)列主元高斯消去法Matlab函数程序:
%列主元高斯消去法
functionx=GAUSSPP(A,B)
n=length(B);
%--------消去过程------------
%找出列主元
fork=1:
n-1
MAX=abs(A(k,k));
Rmax=k;
fori=k+1:
ifabs(A(i,k))>
MAX
MAX=abs(A(i,k));
Rmax=i;
ifA(Rmax,k)==0
break;
%交换当前行和主元所在行
temp=A(k,j);
A(k,j)=A(Rmax,j);
A(Rmax,j)=temp;
temp=B(k);
B(k)=B(Rmax);
B(Rmax)=temp;
A(i,k)=A(i,k)/A(k,k);
forj=k+1:
A(i,j)=A(i,j)-A(i,k)*A(k,j);
B(i)=B(i)-A(i,k)*B(k);
A
%--------回代过程----------------
x=zeros(1,n);
x(n)=B(n)/A(n,n);
fork=n-1:
S=B(k);
S=S-A(k,j)*x(j);
x(k)=S/A(k,k);
5.3.2压缩带状对角方程组
(1)dat53.dat主程序:
%压缩格式dat143.dat
\2015上机题目\实验五\dat53.dat'
D=fread(fid,5,'
n=D(3);
%矩阵维数
p=D(5);
%上带宽
q=D(4);
%下带宽
A=fread(fid,[p+q+1,n],'
B=fread(fid,[n,1],'
C=A'
;
D=B'
%采用求解压缩型的列主元Gauss消去法
x=GAUSSPP1(C,D,n,p,q);
(2)针对压缩格式的带状矩阵,尤其是大型矩阵的列主元Gauss消去法子程序GAUSSPP1如下:
%针对压缩格式的带状方程组的列主元Gauss消去法
functionx=GAUSSPP1(A,B,n,p,q)
%-------消去过程---------
pmax=k+p;
ifpmax>
n
pmax=n;
pmax
A(i,k-i+p+1)=A(i,k-i+p+1)/A(k,p+1);
qmax=k+q;
ifqmax>
qmax=n;
qmax
A(i,j-i+p+1)=A(i,j-i+p+1)-A(i,k-i+p+1)*A(k,j-k+p+1);
B(i)=B(i)-B(k)*A(i,k-i+p+1);
x(n)=B(n)/A(n,p+1);
sum=0;
sum=sum+A(k,j-k+p+1)*x(j);
x(k)=(B(k)-sum)/A(k,p+1);
5.4实验结果及分析
5.4.1Matlab运行结果
方程组一:
文件标识符为id=F1E1D1A0
文件版本号为ver=102
矩阵A的阶数n=15
矩阵A的上带宽q=6
矩阵A的下带宽p=6
方程组的的解为:
1.0000
方程组二:
矩阵A的阶数n=2050
矩阵A的下带宽p=5
1.9600
方程组三:
文件版本号为ver=202
矩阵A的阶数n=43512
矩阵A的上带宽q=4
矩阵A的下带宽p=3
3.1413
5.4.2总结分析
对于任意一个方程组,所有根都相等,因此Matlab求解结果如下表所示:
方程
阶数
根
dat141.dat
1.0000
dat142.dat
2050
1.9600
dat143.dat
43512
3.1413
由上面结果可以看出,所编写程序可以解决正常系数矩阵方程组,同时能解决带状系数矩阵方程组,并且在解决该类系数矩阵时具有很高的效率,说明了列主元消去法在这一类方程组的有效性。
如果要继续提高计算的精度,可以采用全主元消去法。
但是此时的工作量要大的多。
5.5本专业算例
本专业算例:
一底边绝热的正方形导热物体,内无热源,其余三边温度如下,求解该正方形内节点1、2、3、4的温度。
解:
导热物体无内热源,物性为常数的控制方程可以写成
考虑到空间步长相同时,离散方程可以写成:
式中aE=1,aW=1,aN=1,aS=1,aP=4
从而有:
所有的节点方程可整理为:
用列主元高斯消去法对上述矩阵进行求解,其程序如下:
A=[-4110-70;
1-401-50;
10-31-15;
011-3-10];
[m,n]=size(A);
(m-1)
max=i;
fork=i+1:
m
ifabs(A(k,i))>
abs(A(max,i))
max=k;
end
temp=A(i,i:
m);
A(i,i:
m)=A(max,i:
A(max,i:
m)=temp;
forj=(i+1):
m
A(j,:
)=A(j,:
)-A(i,:
)*A(j,i)/A(i,i);
end
回代'
)
x(m)=A(m,n)/A(m,m);
fork=(m-1):
1
x(k)=(A(k,n)-A(k,k+1:
m)*x(k+1:
m)'
)/A(k,k);
x
运算结果为T1=28.7368,T2=24.2632,T3=20.6842,T4=18.3158。
运算结果和正确结果吻合。
学习感悟
学习是一件忙碌的事情,但主旋律是充实和快乐。
在学期末,首先要感谢的是老师这学期对我们的教学。
您板书的教学方式,能让我们更快的理解知识,面对面交流,会让我们对知识更深一步的掌握。
编程的过程很乏味,刚入学的我也许还没有特别深入的掌握一门编程语言,所以我去借阅书籍,咨询师长,通过这样的过程,自身的编程水平也得到了一定的提升,同时也进一步掌握了我们所学的知识,这个过程让人快乐。
值此元旦来临之际,祝老师节日快乐。
报告中如有不足之处,还希望老师批评指正。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 西安交通大学 计算方法 作业 资料 doc