西安交大计算方法b大作业Word文档下载推荐.docx
- 文档编号:21091516
- 上传时间:2023-01-27
- 格式:DOCX
- 页数:21
- 大小:353.40KB
西安交大计算方法b大作业Word文档下载推荐.docx
《西安交大计算方法b大作业Word文档下载推荐.docx》由会员分享,可在线阅读,更多相关《西安交大计算方法b大作业Word文档下载推荐.docx(21页珍藏版)》请在冰豆网上搜索。
已探测到一组等分点位置的深度数据(单位:
米)如下表所示:
分点
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)估算所需光缆长度的近似值,并作出铺设河底光缆的曲线图;
算法思想:
由于题中所给点数为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:
fori=21:
(k+1)
d(i)=(d(i)-d(i-1))/(x(i)-x(i-k));
end
%%%设定d的边界条件
fori=2:
d(i)=6*d(i+1);
end
d
(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:
21
l(k)=a(k)/u(k-1);
u(k)=b(k)-l(k)*r(k-1);
yy(k)=d(k)-l(k)*yy(k-1);
%%%赶
m(21)=yy(21)/u(21);
fork=20:
m(k)=(yy(k)-r(k)*m(k+1))/u(k);
%%%绘制曲线
k=1;
nn=100;
xx=linspace(0,20,nn);
l=0;
forj=1:
nn
ifxx(j)<
=x(i)
k=i;
break;
else
k=i+1;
h=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);
%求解光缆长度
%%%绘图
plot(xx,s,'
r-'
linewidth'
1.5)
disp(['
¹
光缆长度为ª
'
num2str(l(nn+1)),'
Ã
×
])
曲线图如图2-1所示,计算光缆长度如图2-2所示。
图2-1光缆插值曲线图
图2-1光缆计算长度显示
3.假定某天的气温变化记录如下表所示,试用数据拟合的方法找出这一天的气温变化的规律;
试计算这一天的平均气温,并试估计误差。
时刻
平均气温
23
25
28
22
24
31
34
29
27
此题中所给数据点数目较多,采用拉格朗日插值法或者牛顿插值法需要很高次的多项式,计算困难,误差大;
采用样条插值计算量虽然不大,但是存放参数Mi的量很大,且没有一个统一的数学公式来表示,也不是很方便。
所以可考虑用最小二乘法进行拟合。
计算过程中,分别使用二次函数、三次函数以及四次函数,计算其相应的系数,估算误差并作图比较各个函数之间的区别。
(参考课本P123)
1.1[形成矩阵Qk]
1.2[变换Gk-1到Gk]
2.[求解三角方程]
3.[计算误差]
源代码:
x=0:
24;
y=[15
14
15
16
18
20
23
25
28
31
34
29
27
24
22
17
16];
m=length(x);
n=input('
请输入函数的次数
plot(x,y,'
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
for
j=1:
n
i=1:
m
G(i,j)=x(1,i)^(j-1);
end
%%%建立矩阵Qk
k=1:
i=k:
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:
w(j)=G(j,k);
b(k)=c(k)*w(k);
%%%变换矩阵Gk-1到Gk
G(k,k)=c(k);
n+1
q=w(i)*G(i,j)+q;
s=q/b(k);
G(i,j)=s*w(i)+G(i,j);
%%%求解三角方程
Rx=h1
a(n)=G(n,n+1)/G(n,n);
i=n-1:
(-1):
j=i+1:
f=G(i,j)*a(j)+f;
a(i)=(G(i,n+1)-f)/G(i,i);
%a(i)存放各级系数
a
%%%回代过程
p=zeros(1,m);
p(j)=p(j)+a(i)*x(j)^(i-1);
plot(x,p,'
r*'
x,p,'
E2=0;
%用E2来存放误差
%%%误差求解
i=n+1:
E2=G(i,n+1)^2+E2;
E2=E2^0.5;
disp('
误差为'
disp(E2);
t=0;
t=t+p(i);
t=t/m;
%%%平均温度
平均温度为'
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)
holdon
plot(i,0,'
)%画出原函数图像
%%Newton法求根
x1=input('
输入初值'
e=10^(-4);
%误差设定
Nmax=1000;
%迭代最大次数限定
forn=1:
Nmax
f0=f2(x1);
ifabs(f2(x1))<
e
fprintf('
输出的f(x)已经足够小'
x=x1;
break
F0=f3(x1);
x=x1-f0/F0;
ifabs(x-x1)<
x1=x;
fprintf('
输出方程的根x=%2f'
x)
计算结果:
函数图像如图4-1所示。
计算结果分别见图4-2所示。
图4-1函数图像
图4-2计算结果
根据带入不同的初值,可以求出不同的根,有图4-2可以看出,原函数的根大约有三个,分别是-0.654542、0.681174、1.870799。
5.线性方程组求解。
(1)编写程序实现大规模方程组的高斯消去法程序,并对所附的方程组进行求解。
所附方程组的类型为对角占优的带状方程组。
(2)针对本专业中所碰到的实际问题,提炼一个使用方程组进行求解的例子,并对求解过程进行分析、求解。
算法思想:
高斯消去法是利用现行方程组初等变换中的一种变换,将一个不为零的数乘到一个方程后加到另一个方程,使方程组变成同解的上三角方程组,然后再自下而上对上三角方程组求解。
1.读取二进制文件,存入计算矩阵
2.对矩阵进行初等变换,然后求解(详见计算方法教程第2版高斯消去法以及列主元高斯消去法算法)
%%读取系数矩阵
[f,p]=uigetfile('
*.dat'
选择数据文件'
%读取数据文件
num=5;
%输入系数矩阵文件头的个数
name=strcat(p,f);
file=fopen(name,'
r'
head=fread(file,num,'
uint'
%读取二进制头文件
id=dec2hex(head
(1));
%读取标识符
文件标识符为'
id
ver=dec2hex(head
(2));
%读取版本号
文件版本号为'
ver
n=head(3);
%读取阶数
矩阵A的阶数'
q=head(4);
%上带宽
矩阵A的上带宽'
q
p=head(5);
%下带宽
矩阵A的下带宽'
p
dist=4*num;
fseek(file,dist,'
bof'
%把句柄值转向第六个元素开头处
[A,count]=fread(file,inf,'
float'
%读取二进制文件,获取系数矩阵
fclose(file);
%关闭二进制头文件
%%对非压缩带状矩阵进行求解
ifver=='
102'
a=zeros(n,n);
fori=1:
n,
a(i,j)=A((i-1)*n+j);
%求系数矩阵a(i,j)
b=zeros(n,1);
b(i)=A(n*n+i);
n-1,%列主元高斯消去法
m=k;
fori=k+1:
n,%寻找主元
ifabs(a(m,k))<
abs(a(i,k))
m=i;
ifa(m,k)==0%遇到条件终止
disp('
错误!
return
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;
n,%计算l(i,k)并将其放到a(i,k)中
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);
x=zeros(n,1);
%回代过程
x(n)=b(n)/a(n,n);
fork=n-1:
1,
x(k)=(b(k)-sum(a(k,k+1:
n)*x(k+1:
n)))/a(k,k);
%%对压缩带状矩阵进行求解
202'
%高斯消去法
m=p+q+1;
a=zeros(n,m);
1:
n
a(i,j)=A((i-1)*m+j);
b(i)=A(n*m+i);
%求b(i)
(n-1)%开始消去
ifa(k,(p+1))==0
st1=n;
if(k+p)<
st1=k+p;
fori=(k+1):
st1
a(i,(k+p-i+1))=a(i,(k+p-i+1))/a(k,(p+1));
forj=(k+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);
b(i)=b(i)-a(i,k+p-i+1)*b(k);
%回代
x(n)=b(n)/a(n,p+1);
sum=0;
fork=(n-1):
sum=b(k);
st2=n;
if(k+q)<
st2=k+q;
st2
sum=sum-a(k,j+p-k+1)*x(j);
x(k)=sum/a(k,p+1);
方程组的的解为:
)%输出解
disp(x’)
求解结果
对数据文件dat51求解,结果如下:
文件标识符为
id=F1E1D1A0
文件版本号为
ver=102
矩阵A的阶数
n=15
矩阵A的上带宽
q=3
矩阵A的下带宽
p=3
1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000
对数据文件dat52求解,结果如下:
ver=202
n=20
q=5
p=5
1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000
对数据文件dat53求解,结果如下:
n=2160
q=5
p=5
方程组的解截图如图5-1所示(由于矩阵阶数较大,计算结果未能截取完整)。
图5-1数据文件dat53求解结果
n=43240
q=4
p=4
方程组的解截图如图5-2所示(由于矩阵阶数较大,计算结果未能截取完整)。
图5-2数据文件dat54求解结果
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 西安 交大 计算方法 作业