数值计算方法

1 题目 造倒数表 并例求 18 的倒数 精度为 0 0005 2 算法原理 2 1 牛顿迭代法 牛顿迭代法是通过非线性方程线性化得到迭代序列的一种方法 对于非线性方程 f x 0 若已知根 x 的一个近似值 xk 将 f x 在 xk 处展 成一阶泰勒公式后忽略高次项可得 f x f x k f xk x xk 右端是直线方程 用这个直线方程来近似非线性方程 f x 将非线性方程 f x 0 的根 x 代入 f x 0 即 f x k f xk x xk 0 xk f xk 解出x f xk 将右端取为 xk 1 则 xk 1 是比 xk 更接近于 x 的近似值 即 f xk xk 1 xk f xk 这 就是牛顿迭代公式 相应的迭代函数是 f x x x f x 2 2 牛顿迭代法的应用 11 计算 是求cx 1 0 的解 解出 x 即得到 取 cc 有牛顿迭代公式 精品文档 1欢迎下载1欢迎下载1欢迎下载1欢迎下载1欢迎下载1欢迎下载 cxk 1 1 xk 1 xk c c 这样就失去了迭代的 意义 达不到迭代的效果 1 f x cx 1 f x c 故重新构造方程 cx2 x 0 也是该式的解 故取 f x cx2 x c f x 2cx 1 则有牛顿迭代公式 xk 1 xk cxk2 xk cxk2 k 0 1 2cxk 12ck 1 11 的值在 之间 取初值 x0 0 1 2010 3 流程图 4 输出结果 5 结果分析 当k 3 时 得 5 位有效数字 0 05 564 此时 x3 x4 0 00 000 x 时必收敛 但是当 x0 0 时迭代结果发散 较小尚不确定 6 心得体会 起初对题目的理解并不是很透彻 另外对构建牛顿迭代公式理论依据不是特 别充分 比如说为什么在原有直接得到的式子两边各乘一个 x 只是试出来的 在编程方面不够成熟 当然也加深了对牛顿迭代法的理解和应用的具体实现 实验二 例 3 4 1 题目 用列主元消去法求解方程组 12x1 3x2 3x3 15 18x1 3x2 x3 15 x1 x2 x3 6 并求出系数矩阵 A 的行列式的值 det A 2 算法原理 2 1 顺序高斯消去法 顺序高斯消去法是利用线性方程组初等变换中的一种变换 即用一个不为零 的数乘一个方程后加至另一个方程 使方程组变成同解的上三角方程组 然后再 自下而上对上三角方程组求解 这样 顺序高斯消去法可分成 消去 和 回代 两个过程 在用顺序高斯消去法时 在消元之前检查方程组的系数矩阵的顺序主子式 当阶数较高时是很难做到的 若线性方程组的系数具有某种性质时 如常遇到的 对角占优方程组 自然能够用高斯消去法求解 2 2 列选主元消去法 线性方程组只要系数矩阵非奇异 就存在惟一解 但是按顺序消元过程中可 能出现主元素akk k 0 这时尽管系数矩阵非奇异 消元过程无法再进行 或 者 即使akk k 0 但如果其绝对值很小 用它作除数也会导致其他元素的数量级急 剧 精品文档 3欢迎下载3欢迎下载3欢迎下载3欢迎下载3欢迎下载3欢迎下载 增大和使舍入误差扩大 将严重影响计算的精度 为避免在校园过程确定乘数时的所用除数是零或绝对值小的数 即零主元或 小主元 在每一次消元之前 要增加一个选主元的过程 将绝对值大的元素交换 到主对角线的位置上来 列选主元是当高斯消元到第 k 步时 从 k 列的akk 以下 包括akk 的各元 素中选出绝对值 大的 然后通过行交换将其交换到akk 的位置上 交换系数矩阵 中的两行 包括常数项 只相当于两个方程的位置交换了 因此 列选主元不 影响求解的结果 列选主元消去法常用来求行列式 设有矩阵 a11 L a1n A MM an1 L ann 用列主元消去法将其化为上三角形矩阵 对角线上元素为a11 1 a22 2 L a33 3 于是行列式 det A 1 ma a11 1 22 2 Lann n 其中 m 为所进行的行交换次数 这是实际中求行列式值的可靠方法 精品文档 4欢迎下载4欢迎下载4欢迎下载4欢迎下载4欢迎下载4欢迎下载 3 流程图 精品文档 5欢迎下载5欢迎下载5欢迎下载5欢迎下载5欢迎下载5欢迎下载 4 输出结果 5 结果分析 采用计算机运算在计算大数据时有明显的优点 另外也需要考虑到存储 高斯消去法的使用条件是akk k 0 k 1 2 L n 而列选主元法可以保证这一 条件 并且可以避免在消元过程确定乘数时所用除数是绝对值小的数 相对全选 主元的运算量小 一般也可以满足精度要求 6 心得体会 此次上机不仅需要对原理了解透彻 而且要求的编程能力较强 在定义和思 路上没问题 只是在编程软件的使用上遇到些不稳定的问题 如头文件的使用 在存储空间上得到了新的认识 另外发现了当代码多时流程框图的好处 编程是 一件很需要耐心的事 自己还有很大进步空间 实验三 例 3 10 1 题目 用杜里特尔分解法求矩阵 A的逆矩阵 A 1 11 1 A 12 2 2 11 2 算法原理 2 1 杜里特尔分解法设线性方程组Ax b 对系数矩阵 A进行除不交换两行 位置得初等行变换相当于用初等矩阵M1 左乘 A 在对方程组第一次消元后 A 1 和b 1 分别化为 A 2 和 精品文档 6欢迎下载6欢迎下载6欢迎下载6欢迎下载6欢迎下载6欢迎下载 b 2 即 M1 A 1 A 2 M 1b 1 b 2 1 m211 其中M1 m311 MO mn11 第 k 次消元时 A k 和b k 分别化为 A k 1 和b k 1 即 M A k k A k 1 M b k b k 1 k 1 其中 M1 O 1 mk 1 k M mnk O O 1 消元过程是对k 1 n 1 进行的 因此有 Mn 1LM M A2 1 1 A n M LM M b2 1 b n n 11 将上三角形矩阵 A n 记为 U 于是有 精品文档 7欢迎下载7欢迎下载7欢迎下载7欢迎下载7欢迎下载7欢迎下载 A M M1 12 1LMn 11U LU 其中 1 m211 mm1 L M M1 12 1LMn 11 3132 MMm43 O MMMO mn1mn2mn3 L 为单位下三角形矩形 1 这样高斯消去法的实质是将系数矩阵 A 分解为两个三角形矩阵 L 和 U 相 乘 即A LU 在上述矩阵描述中遇到了下三角形矩阵运算 主对角线以上元素全为零的 方阵称为下三角形矩阵 下三角形矩阵的乘积仍是下三角形矩阵 若下三角形 矩阵可逆 其逆矩阵仍是下三角形矩阵 而且下三角形矩阵的乘积和逆矩阵很 容易求得 把 A 分解成一个单位下三角阵和一个上三角阵 U 的乘积成为杜里特尔分解 这种分解是惟一的 2 2 高斯 约当法 高斯消去法有消元和回代两个过程 当对消元过程稍加改变便可以使方程组 化为对角形方程组 Dx b 的形式 其中矩阵D为对角形矩阵 即 a11 1 2 D a22 O ann n 当高斯 约当消去法消元的每一步都先用主元去除其所在行的各元素 包括 常数项 时 方程组便可化成 1 x1 b1 n 精品文档 8欢迎下载8欢迎下载8欢迎下载8欢迎下载8欢迎下载8欢迎下载 1 O xM2 bM2 n 1 xn n bn 这是等号右端即为方程组的解 高斯 约当消去法每一步都用主元去除其所在行 的各元素 包括常数项 这个个过程成为归一化 这时方程组的系数阵转化为 单位阵 为减小误差 高斯 约当消去法还常用列选主元技术 精品文档 9欢迎下载9欢迎下载9欢迎下载9欢迎下载9欢迎下载9欢迎下载 精品文档 10欢迎下载10欢迎下载10欢迎下载10欢迎下载10欢迎下载10欢迎下载 4 输出结果 5 结果分析 采用杜里特尔分解法求解方程组时 由于把对系数矩阵的计算和对右端项的 计算分开 这使计算线性方程组系非常方便 只需进行一次矩形三角分解 然后 再解多个三个方程组 且多解一个方程组仅需要增加大约n2 次乘除法运算 采用高斯约当法仅需要进行消元归一 而不需要回代 为编程实现提供了便 利 6 心得体会 步骤很重要 审题 确定算法 解题步骤 流程图 程序 代入简单值进行 验证 在编程时先在代码输入区打好框架 并且尽量在每一命令后注释 方便检 查错误和日后复习 定义和变量存储很灵活 如我把单位向量直接赋给了 A 矩 阵变量中 还有根据 终的目的直接简化计算 另外赋值前 确定存储空间并且 要定义初值为零 实验四 例 4 6 1 题目已知 f x 的 观测数据 x1234 f x0 5 63 构造插值多项式 精品文档 11欢迎下载11欢迎下载11欢迎下载11欢迎下载11欢迎下载11欢迎下载 2 算法原理 首先构造基函数l xk i n0 xxk xxii 可以证明基函数满足下列条件 i k 0i k l xk i 1 i k 对于给定n 1 个节点 n次拉格朗日插值多项式由下式给出 n L x l x yk k k 0 其中l xk i n0 xxk xxii i k 由于l xk 是一个关于 x的n次多项式 所以L x 为关于 x的不高于n 次的代数多项式 当 x xi 时 L x i yi 满足插值条件 3 流程图 4 输出结果 精品文档 12欢迎下载12欢迎下载12欢迎下载12欢迎下载12欢迎下载12欢迎下载 5 结果分析 由于所知的拉格朗日计算机算法只能实现计算某一特定值的近似函数值 而 不知如何导出表达式 故例求 x 2 5 处的函数值以说明表达式以得出 只是在 计算机程序中 并且也能达到拉格朗日插值法使用的目的 6 心得体会 编程不够细心 程序没问题 却因为不知道是输入文件错了而检查了好长时 间 但同时也加深了对拉格朗日基函数性质的认识和理解 实验五 习题 5 2 1 题目 给出平面函数 z x y ax by c的数据 i12345 xi0 10 20 40 60 9 yi0 20 30 50 70 8 zi0 580 630 730 830 92 按 小二乘原理确定a b c 2 算法原理 2 1 小二乘原理 设已知某物理过程 y f x 的一组观测数据 xi f x i i 1 2 L m 要求在某特定函数类 x 中寻找一个函数 x 作为 y f x 的近似函 数 使得二者在 xi 上的误差或称残差 i xi f x i i 1 2 L m 按某种度量标准为 小 这就是拟合问题 精品文档 13欢迎下载13欢迎下载13欢迎下载13欢迎下载13欢迎下载13欢迎下载 要求残差 i 按某种度量标准为 小 即要求由残差 i 构成的残差向量 0 1 L m T 的某种范数 为 小 例如 要求 或 即 mm 1 i xi f x i 0i 0 max i max xi f x ii 为 小 这本来都是很自然的 可是计算不太方便 所以通常要求 11 2 m i2 2 m xi f x i 2 2 i 0 i 0 或者 mm 22 i2 xi f x i 2 i 0i 0 为 小 这种要求误差平方和 小的拟合称为曲线拟合的 小二乘法 2 2 多变量数据拟合 对于给定的一组数据 xi yi i 1 2 m 寻求做n次多项式 n y a xkk k 0 使性能指标 J a a 0 1 L an Error Error yi Error Error a xk ik 2 为 小 i 1k 0 由于性能指标J可以被看做关于ak k 0 1 n的多元函数 故上述拟 合多项式的构造问题可转化为多元函数的极值问题 令 精品文档 14欢迎下载14欢迎下载14欢迎下载14欢迎下载14欢迎下载14欢迎下载 J 0 ak xinxin 1xin 2 L xi2n an xin yi 对多变量 或称多元 线性模型 y a0 a x1 1 a x22 L a xn n 进行了 m 次观测 y1 a0 a x1 11 a x221 L a xn1n y2 a0 a x1 12 a x222 L a xn2n M 从而有正则方程组 m xixi2 3 xixi2 xi M MM L L L xin a0 yi n 1 a1 xi yi xi M M M 精品文档 15欢迎下载15欢迎下载15欢迎下载15欢迎下载15欢迎下载15欢迎下载 yn a0 a x1 1n a x22n L a xnnn 这个称为回归方程组 写成矩阵形式 y A y1 1x11 其中 y y 2 A A 1x12 y M MM 1x1m ym x21 L xn1 a0 x22 L xn2 a1 M M M M x2m L xnm an 当m n 时 要确定一组a ii 0 1 L n 使之精确地满足 m 个方程 这是 超定方程组的问题 只能在 小平方误差的基础上确定 i 定义残差向量 1 2 L m T 则 y Ay A 其中 y y y1 2 L ym T 代替输出向量 取性能指标 J T 使之 小 以此确定出 由 J T y A A T y A A y yT TA Ay yTA A TA A A AT 利用向量和矩阵的运算公式 有 A A A AT A AT y 此即为正则方程组 当 A A A AT 非奇异时 可求得 A A A AT 1A AT y 精品文档 16欢迎下载16欢迎下载16欢迎下载16欢迎下载16欢迎下载16欢迎下载 3 流程图 4 输出结果 5 结果分析 曲线拟合的 小二乘法是反映所给数据点的总的趋势 并不是严格的通过每 个数据点 这样就避免了大量数据插值时需要高次多项式 同时又去掉了数据 所含的测量误差 6 心得体会 整理思路越来越熟练了 所以执行各个步骤也相对简单了很多 另外对原 理也加深了认识 精品文档 17欢迎下载17欢迎下载17欢迎下载17欢迎下载17欢迎下载17欢迎下载 附录 1 造倒数表 1 源程序 include stdio h include math h define N 30 void main int i float x N c FILE fp1 fp2 fp1 fopen input1 txt r fp2 fopen output1 txt w fscanf fp1 f fscanf fp1 f 初值 fprintf fp2 倒数表 n for i 0 i N i 牛 顿迭代法 x i 1 x i x i c 2 c x i 1 fprintf fp2 k d tx d 5f n i i x i if fabs x i 1 x i 输入文件 input1 txt 18 0 1 3 输出文件 output1 txt 倒数表 k 0 x 0 0 10000 k 1x 1 0 06923 k 2x 2 0 05781 k 3x 3 0 05564 k 4x 4 0 05564 计算结果 1 18 000000 0 056 2 例 3 4 1 源程序 include iostrea m include cmath using namespace std define N 10 void main int i j k l n float b N a N N t d det 1 0 精品文档 19欢迎下载19欢迎下载19欢迎下载19欢迎下载19欢迎下载19欢迎下载 数据输入 FILE fp1 fp2 fp1 fopen input2 txt r fp2 fopen output2 txt w fscanf fp1 d for i 0 i n i for j 0 j n j fscanf fp1 f for i 0 i n i fscanf fp1 f 数据输入 高斯消去 消元 列选主元函数 for k 0 k n 1 k 从第一次消 元到第 N 1 次消元 d a k k l k for i k 1 ifabs d d a i k i l if i n 判断是否奇异 不奇异进行行交换 精品文档 20欢迎下载20欢迎下载20欢迎下载20欢迎下载20欢迎下载20欢迎下载 if d 0 fprintf fp2 奇异 如果所有行的 首列 都为 0 为奇异 else if l k 如果第 k 行的 首列 并不是 大 det det 1 for j k j n j 交换系数矩阵中的两行 t a l j a l j a k j a k j t t b l b l b k b k t 交换右端常向量中的两行 列选主元函数 for i k 1 i n i 第 k 1 次消元要得到 N k 1 个乘数 a i k a i k a k k for j k 1 j 0 i 从倒数第二项开始依次回代 N 1 次 t 0 for j i 1 j n j t t a i j b j b i b i t a i i 高斯消去 数据输出 for i 0 i n i 输出方 程组的解 fprintf fp2 x d 4f n i 1 b i for i 0 i输入文件 input2 txt 3 12 3 3 精品文档 22欢迎下载22欢迎下载22欢迎下载22欢迎下载22欢迎下载22欢迎下载 18 3 1 1 1 1 15 15 6 3 输出文件 output2 txt x 1 1 0000 x 2 2 0000 x 3 3 0000 detA 66 0000 3 例 3 10 1 源程序 include iostrea m include cmath using namespace std define N 30 void main int i j r k n float a N N 0 s 数据输入 FILE fp1 fp2 fp1 fopen input3 txt r fp2 fopen output3 txt 精品文档 23欢迎下载23欢迎下载23欢迎下载23欢迎下载23欢迎下载23欢迎下载 w fscanf fp1 d for i 0 i n i for j 0 j n j fscanf fp1 f for i 0 i n i 输入单位阵 j i n a i j 1 数据输入 LU 分解 for i 1 i n i r 0 时 1 区间不变化 2 区间变化 a i 0 a i 0 a 0 0 for i 1 i n i 第二行列至末行列进行变化时 for j i j 2 n j 第 2 r 1 1 区间的 变化行 s 0 0 for k 0 k i k 对应行列元素的乘积进行求和 s a i k a k j a i j a i j s for j i 1 j n j 2 r 1 区间的变化列 s 0 0 for k 0 kLUY LU 分解 高斯约当法解 Ux Y 提取 UY 减少计算 for i 1 i n i for j 0 j i j a i j 0 提取 UY 减少计算 消元 for j 0 j 2 n j 首 行归一化 a 0 j a 0 j a 0 0 a 0 0 1 第一列其余行已为零 for i 1 i n i n 1 次归一消元 for j i 1 j 2 n j 第 i 1 行的各列进 行归一化 a i j a i j a i i a i i 1 for r 0 r i r 对第一行至第 i 1 行的 i 列进行置零 精品文档 25欢迎下载25欢迎下载25欢迎下载25欢迎下载25欢迎下载25欢迎下载 for j r 2 j 2 n j r 行的各列与第 r 1 行的对应列进行减法运算 a r j a r j a i j a r r 1 第 r 1 行归一后 乘数即为要置零处的值 a r r 1 0 乘数用 完之后置零即可 消元 高斯约当法解 Ux Y 数据输出 fprintf fp2 A 的逆矩阵为 n n for i 0 i n i fprintf fp2 for j 0 j输入文件 input3 txt 精品文档 26欢迎下载26欢迎下载26欢迎下载26欢迎下载26欢迎下载26欢迎下载 3 1 1 1 1 2 2 2 1 1 3 输出文件 output3 txt A 的逆矩阵为 2 0000 1 0000 0 0000 1 5000 0 5000 0 5000 2 5000 1 5000 0 5000 4 例 4 6 1 源程序 include stdio h include math h define N 50 int n i j k float xx 0 0 yy 0 0 t x N y N c N A N main FILE fp1 fp2 fp1 fopen input4 txt r fp2 fopen output4 txt w fscanf fp1 d n 4 for i 0 i n i 精品文档 27欢迎下载27欢迎下载27欢迎下载27欢迎下载27欢迎下载27欢迎下载 fscanf fp1 f f 1 0 2 5 3 6 4 3 fscanf fp1 f 要计算的值 for k 0 k n k 拉格 朗日插值 t 1 0 for j 0 j k j t xx x j t x k x j for j k 1 j输入文件 input4 txt 4 1 0 2 5 3 6 4 3 2 5 3 输出文件 output4 txt x 2 5000000 处的函数值为 y 6 3750000 5 习题 5 2 精品文档 28欢迎下载28欢迎下载28欢迎下载28欢迎下载28欢迎下载28欢迎下载 1 源程序 include stdio h include math h define N 30 void main int i n k j l float x N y N z N 定义输入变量 float AT 3 N 定义 A 的转置 float X 3 3 Y 3 定义中间变量 ATA 和 ATy float s t d