计算方法第一次上机报告.doc

计算方法实验报告一题目：班级： 计1402 学号：1403140126 姓名： 应成龙 日期： 2016.3.17 程序名： 一、问题和要求：1. 土木工程和环境工程师在设计一条排水渠道时必须考虑渠道的各种参数（如宽度，深度，渠道内壁光滑度）及水流速度、流量、水深等物理量之间的关系。假设修一条横断面为矩形的水渠，其宽度为B，假定水流是定常的，也就是说水流速度不随时间而变化。为了不同的工业目的（比如说要把污染物稀释到一定的浓度以下，或者为某工厂输入一定量的水），需要指定流量Q和B，求出水的深度。这样，就需要求解其中Q 是水的流量（），U是流速（），H是水的深度（）。 一个具体的案例是.求出渠道中水的深度H。2. 分析用迭代法求解方程组的收敛性，并求出使 | x(k + 1) x(k) |1 0.0001 的近似解及相应的迭代次数。考虑用（1） 雅可比迭代法；（2）赛德尔迭代法；（3）超松驰迭代法（w 取1.2，1.3，1.9，0.9）。3. 对下面的9阶对称矩阵分别做LU分解和LDLT分解二、设计总体方案及分析：2.1 问题分析第1题考查的是非线性方程组的迭代法（我用的是牛顿迭代法）第2题考查的是线性代数方程组的迭代法，分别是（1）雅可比迭代法，（2）赛德尔迭代法， （3）超松驰迭代法第3题考查的是线性代数方程组的矩阵三角分解法，分别是LU分解和LDLT分解2.2 概要设计2.3 详细设计（二合一）第1题用牛顿迭代法， 先判断f(x)是否收敛，从x=0开始，利用公式不断迭代，直至Xk+1-Xk的值小于0.000001，才停止。并输出所有的迭代值第2-1题用雅可比迭代法时，由于公式中，Xj迭代的次数是k次，所以需要用一个额外的数组来保存之前第k次迭代结果，我用了x6N= 0 ，来存放那个数组，通过控制每个X k+1次迭代与X k次迭代差=0.001控制循环次数，当每个差都同时=0.001，循环停止，迭代结束第2-2题用的是赛德尔迭代法，由于已经进行过的迭代，直接覆盖原值，所以赛德尔迭代法程序设计比雅可比迭代要简单，不用在设计数组存放X的k次迭代值，而且只需在在雅可比迭代的基础上稍作修改就可以，而且循环控制方法与雅可比迭代相同第2-3题用的是超松弛迭代法，超松弛法在赛德尔迭代法基础上，加入w，只需稍作修改就好，循环控制与上面几题相同，这3题互相关联，只要写出一个程序，其余2个就一定会写了。在难度上雅可比迭代稍难一些。在程序中只要修改宏定义w就行了。j = k, , n第3-1题用的是LU分解i = k+1, , n这题我用了三个嵌套循环，在最内层循环有两个并列循环，先给U赋值，再给L赋值，只有U赋值了，L才能赋值，最外层两个循环实现a,l,u的遍历，最终实现全部赋值。第3-2题用的是LDLT分解实际用到的是这3个公式这一题我通过修改我自己写的LU分解，而且只用L和D两个公式行不通，一直出不来结果，借鉴了同学的代码，也是通过3个循环嵌套，实际用的是这3个公式，方法同LU分解一样。2.4 调试分析第1题结果我就算出来两个0，和0.206213，估计是0这个数太接近实际结果了，所以迭代较快。第2-1这题我是这里面用时比较长的一题，一开始是不知道怎么写，然后从公式入手，还得保存X的k次迭代，通过数组的方法来保存，最后终于写出来了。我还把迭代次数及到该迭代完成后X的值在同一行表示出来第2-2,2-3，都是在2-1的基础上稍作修改就很快算出了结果。第3-1也是从公式入手，通过观察公式，发现规律，有了前面的经验之后很快就能写出来。第3-2是由于公式选择的问题，导致一直编译出错，最后在借鉴了同学的代码之后完成。2.5 测试结果第1题第2-1第2-2题第2-3题第3-1第3-2三、源程序及注释第1题#include#include#include#define n 0.03#define B 20#define S 0.0002#define Q 5using namespace std;bool shoulian(double H)if (0.8 / B*pow(n*Q / pow(S, 0.5) / (B + 2 * H), 0.6) 1)return true;else return false;double diedai(double H)return (pow(n*Q / pow(S, 0.5)*pow(B + 2 * H, 2 / 3), 0.6) / B);int main()double H;for (double h = 0; h =0.00000001)cout H t diedai(H) t endl;H = diedai(H);cout H t diedai(H) t endl;第2-1题#include#include#includeusing namespace std;#define N 1000int main()int counter = 0;int a6 = 4, -1, 0, -1, 0, 0, -1, 4, -1, 0, -1, 0, 0, -1, 4, -1, 0, -1, -1, 0, -1, 4, -1, 0, 0, -1, 0, -1, 4, -1, 0, 0, -1, 0, -1, 4 ;double x6N= 0 ;int b6 = 0,5,0,6,-2,6;for (int n = 0; n N; n+)for (int i = 0; i 6; i+)double s = 0;for (int j = 0; j 6; j+)if (j != i)s -= aij * xjn;s += bi;s /= aii;xin + 1 = s;if (fabs(x0n + 1 - x0n) = 0.0001&fabs(x1n + 1 - x1n) = 0.0001&fabs(x2n + 1 - x2n) = 0.0001&fabs(x3n + 1 - x3n) = 0.0001&fabs(x4n + 1 - x4n) = 0.0001&fabs(x5n + 1 - x5n) = 0.0001)counter = n + 1;break;for (int i = 0; i=counter; i+)cout i t x0i t x1i t x2i t x3i t x4i t x5i endl;第2-2题#include#include#includeusing namespace std;#define N 1000int main()int counter = 0;int a6 = 4, -1, 0, -1, 0, 0,-1, 4, -1, 0, -1, 0,0, -1, 4, -1, 0, -1,-1, 0, -1, 4, -1, 0,0, -1, 0, -1, 4, -1,0, 0, -1, 0, -1, 4 ;double x6 = 0 ;double x16 = 0 ;int b6 = 0, 5, 0, 6, -2, 6 ;for (int n = 0; n N; n+)cout n t x0 t x1 t x2 t x3 t x4 t x5 t endl;if (n != 0)if (fabs(x0 - x10) = 0.0001&fabs(x1 - x11) = 0.0001&fabs(x2 - x12) = 0.0001&fabs(x3 - x13) = 0.0001&fabs(x4 - x14) = 0.0001&fabs(x5 - x15) = 0.0001)counter = n + 1;break;for (int i = 0; i 6; i+)x1i = xi;double s = 0;for (int j = 0; j 6; j+)if (j != i)s -= aij * xj;s += bi;s /= aii;xi = s;第2-3题#include#include#includeusing namespace std;#define N 1000#define w 1.2 /w取1.2,1.3,1.9,0.9int main()int counter = 0;int a6 = 4, -1, 0, -1, 0, 0,-1, 4, -1, 0, -1, 0,0, -1, 4, -1, 0, -1,-1, 0, -1, 4, -1, 0,0, -1, 0, -1, 4, -1,0, 0, -1, 0, -1, 4 ;double x6 = 0 ;double x16 = 0 ;int b6 = 0, 5, 0, 6, -2, 6 ;for (int n = 0; n N; n+)cout n t x0 t x1 t x2 t x3 t x4 t x5 t endl;if (n != 0)if (fabs(x0 - x10) = 0.0001&fabs(x1 - x11) = 0.0001&fabs(x2 - x12) = 0.0001&fabs(x3 - x13) = 0.0001&fabs(x4 - x14) = 0.0001&fabs(x5 - x15) = 0.0001)counter = n + 1;break;for (int i = 0; i 6; i+)x1i = xi;double s = 0;for (int j = 0; j 6; j+)if (j != i)s -= aij * xj;s += bi;s /= aii;s *= w;s += (1 - w)*xi;xi = s;第3-1题#include#include#includeusing namespace std;int main()double a99 = 4,-1,0,-1,0,0,0,0,0, -1,4,-1,0,-1,0,0,0,0, 0,-1,4,0,0,-1,0,0,0, -1,0,0,4,-1,0,-1,0,0, 0,-1,0,-1,4,-1,0,-1,0, 0,0,-1,0,-1,4,-1,0,-1, 0,0,0,-1,0-1,4,-1,0, 0,0,0,0,-1,0-1,4,-1, 0,0,0,0,0,-1,0,-1,4;double l99 = 0 ;double u99 = 0 ;for (int i = 0; i 9; i+)u0i = a0i;lii = 1;if (i != 0)li0 = ai0 / a00;for(int r = 1; r 9; r+)int i = r;double s2 = 0;for (; i 9; i+)double s1 = 0;for (int k = 0; k r; k+)s1 -= lrk * uki;s1 += ari;uri = s1;for (int k = 0; k r; k+)s2 -= lr + 1k * ukr;s2 += ar + 1r;s2 /= urr;lr + 1r = s2;cout l99 endl;for (int i = 0; i 9; i+)for (int j = 0; j 9; j+)cout lij t;cout endl;cout u99 endl;for (int i = 0; i 9; i+)for (int j = 0; j 9; j+)cout uij t;cout endl;第3-2题#include#include#include#define n 9using namespace std;int main()double ann = 4, -1, 0, -1, 0, 0, 0, 0, 0,-1, 4, -1, 0, -1, 0, 0, 0, 0,0, -1, 4, 0, 0, -1, 0, 0, 0,-1, 0, 0, 4, -1, 0, -1, 0, 0,0, -1, 0, -1, 4, -1, 0, -1, 0,0, 0, -1, 0, -1, 4, -1, 0, -1,0, 0, 0, -1, 0, -1, 4, -1, 0,0, 0, 0, 0, -1, 0, -1, 4, -1,0, 0, 0, 0, 0, -1, 0, -1, 4 ;double tnn = 0 ;double lnn = 0 ;double dnn = 0 ;int i, j, k;for ( i = 0;i n; i+) dii = aii;for (k = 0; k i - 1; k+)dii -= tik * lik;for (j = 0; j i; j+)tij = aij;for (k = 0; k i - 1; k+)tij -= tik * ljk;lij = tij / dii;cout l99 endl;for (int i = 0; i 9; i+)for (int j = 0; j 9; j+)cout setw(10)lij ;cout endl;cout d99 endl;f