西安交大计算方法b大作业.pdf

1 计算方法 B 上机 实验报告 学院 机械工程学院 班级 姓名 学号 2015 年 12 月 22 日 2 1 计算以下和式 0 14211 8184858616 n n S nnnn 要求 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 3 实验结果 11位有效数字计算结果如图1所示 30为有效数字计算结果如图2所 示 图 1 11 位有效数字计算结果图 2 30 为有效数字计算结果 4 1 某通信公司在一次施工中 需要在水面宽度为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 15 16 17 18 19 20 深度9 15 7 90 7 95 8 86 9 81 10 80 10 93 1 请用合适的曲线拟合所测数据点 2 估算所需光缆长度的近似值 并作出铺设河底光缆的曲线图 算法思想 由于题中所给点数为20 若采用高次多项式插值将产生很大的误差 所以拉格朗日或牛顿并不适用 题中光缆为柔性 可光滑铺设于水底 鉴于此特 性 采用三次样条插值插值法较为合适 算法结构 三次样条算法结构见 计算方法教程 P110 光缆长度计算公式 2020 19 1 2 2 0 00 l 1 f 1 f k k k xdxxdx 源程序 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的边界条件 for i 2 20 5 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 m k h 2 6 xmao h 6 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 1 光缆长度为 a num2str l nn 1 曲线图如图 2 1 所示 计算光缆长度如图2 2 所示 图 2 1 光缆插值曲线图 图 2 1 光缆计算长度显示 7 3 假定某天的气温变化记录如下表所示 试用数据拟合的方法找出这 一天的气温变化的规律 试计算这一天的平均气温 并试估计误差 时刻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 实现思想 此题中所给数据点数目较多 采用拉格朗日插值法或者牛顿插值法需 要很高次的多项式 计算困难 误差大 采用样条插值计算量虽然不大 但是存 放参数 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 8 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 for i 1 n p j p j a i x j i 1 end end 9 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 二次函数拟合结果 10 三次函数拟合 结果如下图所示 图 3 2 三次函数拟合结果 四次函数拟合 结果如下图所示 图 3 3 四次函数拟合结果 11 结果对比 将二次函数 三次函数和四次函数拟合结果绘制在同一个坐标内 如图 3 4 所示 其计算误差结果见表3 1 所示 图 3 4 拟合结果对比分析 12 4 设计算法 求出非线性方程 52 645200 xx的所有实根 并使 误差不超过 4 10 算法思想 本题可采用牛顿法迭代求解 令 52 f x 6 x 45 x 20 得带格式为 k k 1k k f x x x f x 根据函数图像可以找出根的大致分布区间 带入不同的初值即可解出不同的根 源代码 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 e fprintf 输出的 f x 已经足够小 x x1 break else F0 f3 x1 x x1 f0 F0 if abs x x1 e break else x1 x end end end 13 fprintf 输出方程的根x 2f x 计算结果 函数图像如图 4 1 所示 计算结果分别见图 4 2 所示 图4 1 函数图像 图 4 2 计算结果 根据带入不同的初值 可以求出不同的根 有图 4 2 可以看出 原函数的根 大约有三个 分别是 0 654542 0 681174 1 870799 14 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 的上带宽 15 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 abs a i 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 16 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 n 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 17 st2 n if k q n 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