西安交大计算方法b大作业课件.docx
- 文档编号:14376843
- 上传时间:2023-06-22
- 格式:DOCX
- 页数:28
- 大小:683.84KB
西安交大计算方法b大作业课件.docx
《西安交大计算方法b大作业课件.docx》由会员分享,可在线阅读,更多相关《西安交大计算方法b大作业课件.docx(28页珍藏版)》请在冰点文库上搜索。
西安交大计算方法b大作业课件
《计算方法 B》上机
实验报告
学院:
机械工程学院
班级:
姓名:
学号:
2015 年 12 月 22 日
1
∞
n=0 16
n
⎛ ⎫
⎝ ⎭
(1)若保留 11 个有效数字,给出计算结果,并评价计算的算法;
(2)若要保留 30 个有效数字,则又将如何进行计算。
实现思想:
以上问题出现了近似数相减的问题,为了减小误差,可分别求得减数之和
以及被减数之和,最后将两者相减。
另外,减数与被减数求和均为同号计算,
按照绝对值递增顺序相加可减小舍入误差。
此题中对有效数字有要求,因而计
算时首先需要根据有效数字位数计算得出迭代次数,以保证计算值的精度。
源程序:
m=input('输入有效数字个数m=');
s0=1;s1=0;s2=0;n=0;
%判断迭代次数
while s0>=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
%分别求解各项并求和
for k=n-1:
-1:
0
a1=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;
end
S=vpa(s1-s2,m)
2
实验结果:
11位有效数字计算结果如图1所示;30为有效数字计算结果如图2所
示。
图 1.11 位有效数字计算结果图 2.30 为有效数字计算结果
3
分点
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
15
16
17
18
19
20
深度
9.15
7.90
7.95
8.86
9.81
10.80
10.93
1. 某通信公司在一次施工中,需要在水面宽度为 20 米的河沟底部
沿直线走向铺设一条沟底光缆。
在铺设光缆之前需要对沟底的
地形进行初步探测,从而估计所需光缆的长度,为工程预算提
供依据。
已探测到一组等分点位置的深度数据(单位:
米)如下表
所示:
(1)请用合适的曲线拟合所测数据点;
(2)估算所需光缆长度的近似值,并作出铺设河底光缆的曲线图;
算法思想:
由于题中所给点数为 20,若采用高次多项式插值将产生很大的误差,
所以拉格朗日或牛顿并不适用。
题中光缆为柔性,可光滑铺设于水底,鉴于此
特性,采用三次样条插值插值法较为合适。
算法结构:
三次样条算法结构见《计算方法教程》P110;
光缆长度计算公式:
l =
20
⎰
0
1 + ( f '( x ) )2dx =
19 20
∑ ⎰k ⎰
0
1 + ( f '( x ) )2dx
源程序:
clear;
clc;
x=0:
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];
d=y;
plot(x,y,'k.','markersize',15)
hold on
%%%计算差商
for k=1:
2
for i=21:
-1:
(k+1)
d(i)=(d(i)-d(i-1))/(x(i)-x(i-k));
end
end
%%%设定d的边界条件
4
for i=2:
20
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);
%%%追
for k=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);
end
%%%赶
m(21)=yy(21)/u(21);
for k=20:
-1:
1
m(k)=(yy(k)-r(k)*m(k+1))/u(k);
end
%%%绘制曲线
k=1;
nn=100;
xx=linspace(0,20,nn);
l=0;
for j=1:
nn
for i=2:
20
if xx(j)<=x(i)
k=i;
break;
else
k=i+1;
end
end
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)-
5
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 光缆计算长度显示
6
时刻
0
1
2
3
4
5
6
7
8
9
10
11
12
平均气温
15
14
14
14
14
15
16
18
20
20
23
25
28
时刻
13
14
15
16
17
18
19
20
21
22
23
24
平均气温
31
34
31
29
27
25
24
22
20
18
17
16
3.假定某天的气温变化记录如下表所示,试用数据拟合的方法找出
这一天的气温变化的规律;试计算这一天的平均气温,并试估计误差。
实现思想:
此题中所给数据点数目较多,采用拉格朗日插值法或者牛顿插值法
需要很高次的多项式,计算困难,误差大;采用样条插值计算量虽然不大,但
是存放参数 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';
7
c=zeros(1,n);%建立 c 来存放 σ
q=0;
f=0;
b=zeros(1,m);%建 立 b 用来存放 β %%%形成矩阵 G
for j=1:
n
for i=1:
m
G(i,j)=x(1,i)^(j-1);
end
end %%%建立矩阵 Qk
for k=1:
n
for i=k:
m
c(k)=G(i,k)^2+c(k);
end
c(k)=-sign(G(k,k))*(c(k)^0.5);
w(k)=G(k,k)-c(k);%建立 w 来存放 ω
for j=k+1:
m
w(j)=G(j,k);
end
b(k)=c(k)*w(k);
%%%变换矩阵 Gk-1 到 G k
G(k,k)=c(k);
for j=k+1:
n+1
q=0;
for i=k:
m
q=w(i)*G(i,j)+q;
end
s=q/b(k);
for i=k:
m
G(i,j)=s*w(i)+G(i,j);
end
end
end %%%求解三角方程 Rx=h1
a(n)=G(n,n+1)/G(n,n);
for i=n-1:
(-1):
1
for j=i+1:
n
f=G(i,j)*a(j)+f;
end
a(i)=(G(i,n+1)-f)/G(i,i); %a(i)存放各级系数
f=0;
end
a
%%%回代过程
p=zeros(1,m);
for j=1:
m
8
for i=1:
n
p(j)=p(j)+a(i)*x(j)^(i-1);
end
end
plot(x,p,'r*',x,p,'-');
E2=0;%用 E2 来存放误差
%%%误差求解
for i=n+1:
m
E2=G(i,n+1)^2+E2;
end
E2=E2^0.5;
disp('误差为 ');
disp(E2);
t=0;
for i=1:
m
t=t+p(i);
end
t=t/m;
%%%平均温度
disp(['平均温度为 ',num2str(t),'℃'])
实验结果:
二次函数拟合,结果如下图所示
图 3-1 二次函数拟合结果
9
三次函数拟合,结果如下图所示
图 3-2 三次函数拟合结果
四次函数拟合,结果如下图所示
10
图 3-3 四次函数拟合结果
结果对比:
将二次函数、三次函数和四次函数拟合结果绘制在同一个坐标内,如
图 3-4 所示。
其计算误差结果见表 3-1 所示。
图 3-4 拟合结果对比分析
11
4.设计算法,求出非线性方程 6x5 - 45x2 + 20 = 0 的所有实根,并
使误差不超过10-4 。
算法思想:
本题可采用牛顿法迭代求解,令 f ( x) = 6 x5- 45 x2+ 20 ,得带格式为
xk+1 = xk -
f ( xk )
f '( xk )
根据函数图像可以找出根的大致分布区间,带入不同的初值即可解出不同的根.
源代码:
function y=f2(x)
y=6*x.^5-45*x.^2+20;%定义原函数
function y=f3(x)
y=30*x^4-90*x;%定义原函数倒数
i=-5:
0.1:
5;
y=f2(i);
plot(i,y)
hold on
plot(i,0,'-')%画出原函数图像
%%Newton 法求根
x1=input('输入初值');
e=10^(-4);%误差设定
Nmax=1000;%迭代最大次数限定
for n=1:
Nmax
f0=f2(x1);
if abs(f2(x1)) fprintf('输出的f(x)已经足够小'); x=x1; break else F0=f3(x1); x=x1-f0/F0; if abs(x-x1) break else x1=x; end end end 12 fprintf('输出方程的根x=%2f',x) 计算结果: 函数图像如图4-1所示。 计算结果分别见图4-2所示。 图4-1 函数图像 图 4-2 计算结果 根据带入不同的初值,可以求出不同的根,有图 4-2 可以看出,原函数的 根大约有三个,分别是-0.654542、0.681174、1.870799。 13 5.线性方程组求解。 (1)编写程序实现大规模方程组的高斯消去法程序,并对所附的方 程组进行求解。 所附方程组的类型为对角占优的带状方程组。 (2)针对本专业中所碰到的实际问题,提炼一个使用方程组进行求 解的例子,并对求解过程进行分析、求解。 算法思想: 高斯消去法是利用现行方程组初等变换中的一种变换,将一个不为零的数乘 到一个方程后加到另一个方程,使方程组变成同解的上三角方程组,然后再自下 而上对上三角方程组求解。 算法结构: 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('文件标识符为'); id ver=dec2hex(head (2)); %读取版本号 fprintf('文件版本号为'); ver n=head(3); %读取阶数 fprintf('矩阵 A 的阶数'); n q=head(4); %上带宽 fprintf('矩阵 A 的上带宽'); 14 q p=head(5); %下带宽 fprintf('矩阵 A 的下带宽'); p dist=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-1, %列主元高斯消去法 m=k; for i=k+1: n, %寻找主元 if abs(a(m,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)-a(i,k)*a(k,j); end 15 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); end end %% 对压缩带状矩阵进行求解 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); end end b=zeros(n,1); for i=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) 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); 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); 16 st2=n; if (k+q) 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; end end disp('方程组的的解为: ') %输出解 disp(x’) 求解结果 对数据文件 dat51 求解,结果如下: 文件标识符为 id =F1E1D1A0 文件版本号为 ver =102 矩阵 A 的阶数 n =15 矩阵 A 的上带宽 q =3 矩阵 A 的下带宽 p =3 方程组的的解为: 1.00001.00001.00001.00001.00001.00001.0000 1.00001.00001.00001.00001.00001.00001.0000 1.0000 对数据文件 dat52 求解,结果如下: 文件标识符为 id = F1E1D1A0 文件版本号为 ver =202 矩阵 A 的阶数 n =20 矩阵 A 的上带宽 q =5 矩阵 A 的下带宽 p =5 17 方程组的的解为: 1.00001.00001.00001.00001.00001.00001.0000 1.00001.00001.00001.00001.00001.00001.0000 1.00001.00001.00001.00001.00001.0000 对数据文件 dat53 求解,结果如下: 文件标识符为 id =F1E1D1A0 文件版本号为 ver =102 矩阵 A 的阶数 n =2160 矩阵 A 的上带宽 q =5 矩阵 A 的下带宽 p =5 方程组的解截图如图 5-1 所示(由于矩阵阶数较大,计算结果未能截取完整)。 图 5-1 数据文件 dat53 求解结果 对数据文件 dat53 求解,结果如下: 文件标识符为 id =F1E1D1A0 文件版本号为 ver =202 矩阵 A 的阶数 n =43240 矩阵 A 的上带宽 18 q =4 矩阵 A 的下带宽 p =4 方程组的解截图如图 5-2 所示(由于矩阵阶数较大,计算结果未能截取完整)。 图 5-2 数据文件 dat54 求解结果 19
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 西安 交大 计算方法 作业 课件