计算方法实验指导书.doc

计算方法实验指导书实验一：拉格朗日插值算法目的与要求：熟悉拉格朗日插值多项式构造类型：必做命名（源程序.cpp及工作区.dsw）：lagrange问题：0.40.550.800.910.410750.578150.888111.026521.17520用四次拉格朗日多项式求的函数近似值运行及结果显示： #include #include #define N 4void main()double xN+1=0.4,0.55,0.8,0.9,1;double yN+1=0.41075,0.57815,0.88811,1.02652,1.17520;double varx=0.5;int checkvalid(double x,int n);double Lagrange(double x,double y,double varx,int n);if(checkvalid(x,N)=1)printf(nn插值结果：P(%f)=%fn,varx,Lagrange(x,y,varx,N);else printf(输入的插值节点的x值必须互异！);getch();int checkvalid(double x,int n)/检查插值节点是否互异，有重节点返回-1；否则返回1double Lagrange(double x,double y,double varx,int n)/参数说明：x，y存放插值节点值，varx变量x的值，n插值多项式的次数容错测试：对double xN+1=0.4,0.4,0.8,0.9,1;将能提示错误。实验二：牛顿插值算法目的与要求：熟悉牛顿插值多项式构造，注意与拉氏插值的不同特点类型：必做命名（源程序.cpp及工作区.dsw）：newton_cz问题：0.40.550.800.910.410750.578150.888111.026521.17520用牛顿插值多项式求 运行及结果显示： #include #include #define N 4static double xN+1=0.4,0.55,0.8,0.9,1;static double yN+1=0.41075,0.57815,0.88811,1.02652,1.17520;void main()double varx,fN+1N+1;int checkvalid(double x,int n);void chashang(double x,double y,double fN+1);double compvalue(double tN+1,double x,double varx);varx=0.5;if(checkvalid(x,N)=1)chashang(x,y,f);printf(nn牛顿插值结果：P(%f)=%fn,varx,compvalue(f,x,varx);else printf(输入的插值节点的x值必须互异！);getch();int checkvalid(double x,int n) /检查插值节点是否互异，有重节点返回-1；否则返回1void chashang(double x,double y,double fN+1)/参数说明：x，y存放插值节点值，varx变量x的值，fN+1返回各阶差商值要求：在函数chashang中完成对各阶差商的计算及以下三角格式打印.显示量，本题可采用一维数组仅保存差商表对角线元素，则算法可进一步优化，有兴趣的同学可尝试设计算法。double compvalue(double fN+1,double x,double varx) ./ 参数说明：fN+1传入各阶差商值，varx变量x的值，x传入插值节点值/ 返回结果：插值多项式的值，利用嵌套乘法实现求值运算实验三：自适应梯形求积公式目的与要求：通过实际计算体会自适应梯形公式的原理命名（源程序.cpp及工作区.dsw）：autotrap注意事项：对于中间计算结果不得采用数组来保存，否则空间浪费太严重；类型：必做问题：计算的近似值，使得误差（EPS）不超过运行及结果显示： #include #include #include #define EPS 1e-6void main()double s;double f(double);double AutoTrap(double (*)(double),double,double,double);s=AutoTrap(f,0.0,1.0,EPS);getch();double f(double x)return 4/(1+x*x);double AutoTrap(double (*f)(double),double a,double b,double eps)/参数说明：f为求积函数指针，a,b分别为积分上下限，eps为绝对误差限； 精度控制：实验四：龙贝格求积算法目的与要求：理解并掌握龙贝格求积算法的原理命名（源程序.cpp及工作区.dsw）：romberg类型：选做问题：求的近似值，要求误差不超过运行及结果显示：#include #include #include #define EPS 1e-6double f(double x)return 4/(1+x*x);double Romberg(double a,double b,double (*f)(double),double eps) /参数说明：f为求积函数指针，a,b分别为积分上下限，eps为绝对误差限； void main()printf(the Result is R=%f,Romberg(0,1,f,EPS);实验五：埃特金加速迭代算法目的与要求：掌握埃特金加速算法原理，比较其相对于不定点迭代法加速效果命名（源程序.cpp及工作区.dsw）：aitken类型：必做问题：求方程在，附近的根（精度）注意事项：对于迭代结果，不得采用数组来保存，否则空间浪费太严重；运行及结果显示：#include math.h#include stdio.h#define MaxK 100 /*最大迭代次数*/#define EPS 0.5e-3double g(double x)return x*x*x-1;/迭代函数int aitken(double (*g)(double),double &x,double eps)/参数说明：g传入迭代函数指针， x引用方式传入参数，可返回根值；eps为精度返回值：迭代次数；-1迭代次数越界void main()double x=1.5;int k;k=aitken(g,x,EPS);if(k=-1)printf(迭代次数越界!n);else if(k=-2) printf(除零错误!n);else printf(n经k=%d次迭代，所得方程根为：x=%fn,k,x);getch();精度控制： 实验六：牛顿切线法目的与要求：理解并掌握牛顿切线法解非线性方程求根的原理命名（源程序.cpp及工作区.dsw）：newton_qxf类型：必做问题：分别以,为初始值求方程在1.5附近的根（精度）注意事项：对于迭代结果不得采用数组来保存，否则空间浪费太严重；运行及结果显示：#include math.h#include stdio.h#define MaxK 1000/最大迭代次数#define EPS 0.5e-3double f(double x)return x*x*x-x-1;/函数double f1(double x)return 3*x*x-1;/导函数void main()double x=1.5;int newton(double(*)(double),double(*)(double),double &x0,double);if(newton(f,f1,x,EPS)=1)printf(the root is x=%f n,x);else printf(the method is fail!);int newton(double (*f)(double),double(*f1)(double),double &x0,double eps)/参数说明：f传入求积函数指针，f1传入导函数指针，x0引用方式传入参数，可返回根值；eps为两次迭代值的误差控制精度；函数返回值为：1正常结果；-1为导数出现0；-2迭代次数越界精度控制：实验六：牛顿下山法目的与要求：掌握牛顿下山法，理解其扩大普通牛顿法根收敛区间的原理命名（源程序.cpp及工作区.dsw）：newton_downhill类型：选做问题：求方程在，附近的根（精度）注意事项：对于迭代结果不得采用数组来保存，否则空间浪费太严重；运行及结果显示： #include #include #include #include #define Et 1e-3 /下山因子下界#define E1 1e-3 /根的误差限#define E2 1e-3 /残量精度void main()int b;double f(double);double f1(double);void errormess(int);double x0;int Newton(double (*)(double),double (*)(double),double &x0);x0=0.6;b=Newton(f,f1,x0);/f,f1为函数 if(b=1)printf(nThe root of is x=%fn,x0);else errormess(b);int Newton(double (*f)(double),double (*f1)(double),double &x0)/参数说明：f传入求积函数指针，f1传入导函数指针，x0引用方式传入参数，可返回根值；函数返回值为：1正常结果；-1为导数出现0；-2为下山因子越界；double f(double x)return x*x*x-x-1;double f1(double x)return 3*x*x-1;void errormess(int b)char *mess;switch(b)case-1:mess=f(x)的导数为0！;break;case-2:mess=下山因子已越界，下山处理失败;break;default:mess=其它类型错误！;printf(ther method has failed ! because %s,mess);实验七：正割法算法（弦截法）目的与要求：掌握正割法算法原理 命名（源程序.cpp及工作区.dsw）：ZhengGe类型：必做问题：求方程在附近的根（精度），取注意事项：对于迭代结果不得采用数组来保存，否则空间浪费太严重；运行及结果显示：#include math.h#include stdio.h#include conio.h#define MaxK 1000#define EPS 0.5e-8double f(double x)return x*x*x-x-1;void main()double x0=1,x1=2;int ZhengGe(double (*f)(double),double x0,double &x1,double eps);if(ZhengGe(f,x0,x1,EPS)=1)printf(the root is x=%f n,x1);else printf(the method is fail!);getch();int ZhengGe(double (*f)(double),double x0,double &x1,double eps)/参数说明：f传入函数指针，x0传入x0，x1引用方式传入参数x1并返回根值；eps为精度返回值：迭代次数；-1除0操作，-2迭代次数越界；1方法成功精度控制：实验九：高斯列选主元算法命名（源程序.cpp及工作区.dsw）：colpivot类型：必做问题：求解方程组并计算系数矩阵行列式值运行及结果显示：#include #include stdio.h#include conio.h#define N 3static double aaNN=1,2,-1,1,-1,5,4,1,-2;static double bbN+1=3,0,2;void main()int i,j;double aN+1N+1,bN+1,xN+1,det;double gaussl(double aN+1,double b,double x);for(i=1;i=N;i+)for(j=1;j=N;j+) aij=aai-1j-1;bi=bbi-1;det=gaussl(a,b,x);if(det!=0)printf(n方程组的解为:);for(i=1;i=N;i+)printf( x%d=%f,i,xi);printf(nn系数矩阵的行列式值%f,det);else printf(nn系数矩阵奇异阵，高斯方法失败 !:);getch();double gaussl(double aN+1,double b,double x)/参数说明：a传入增广矩阵（有效元素为），x负责传入迭代初值并返回解向量（有效值为）；返回值：系数矩阵行列式值；void disp(double aN+1,double b)/用于显示选主元及消元运算中间增广矩阵结果，由gaussl()负责调用。实验十：高斯全主元消去算法命名（源程序.cpp及工作区.dsw）：fullpivot类型：选做问题：利用全主元消去法求解方程组运行及结果显示：#include #include stdio.h#include conio.h#define N 3double xN+1;static double aaNN=0.01,2,-0.5,-1,-0.5,2,5,-4,0.5;static double bbN=-5,5,9;void main()int i,j;double aN+1N+1,bN+1,xN+1,det;int tN+1;/引入列交换保存x向量的各分量的位置顺序double gaussl(double aN+1,double b,double x,int t);for(i=1;i=N;i+)for(j=1;j=N;j+) aij=aai-1j-1;bi=bbi-1;for(i=1;i=1;i-)printf( x%d=%f,ti,xi);if(i1)printf( -);printf(nn系数矩阵的行列式值%f,det);else printf(nn系数矩阵奇异阵，高斯方法失败 !:);getch();double gaussl(double aN+1,double b,double x,int t)/参数说明：a传入增广矩阵（有效元素为），x负责传入迭代初值并返回解向量（有效值为）；返回值：系数矩阵行列式值；void disp(double aN+1,double b,int x)显示选主元及消元运算中间增广矩阵结果，由gaussl()负责调用。实验十一：LU分解算法命名（源程序.cpp及工作区.dsw）：LU类型：必做问题：利用LU法求解方程组运行及结果显示：#include math.h#include stdio.h#define N 3static aaNN=3,1,6,2,1,3,1,1,1;static bbN=2,7,4;void main()int i,j;double aN+1N+1,bN+1;void solve(double aN+1,double bN+1);int LU(double aN+1);for(i=1;i=N;i+)for(j=1;j=N;j+)aij=aai-1j-1;bi=bbi-1;if(LU(a)=1)printf(矩阵L如下n);for(i=1;i=N;i+)for(j=1;j=i;j+)if(i=j)printf(%12d ,1);else printf(%12f,aij); printf(n);printf(n矩阵U如下n);for(i=1;i=N;i+)for(j=1;j=N;j+) if(i=j)printf(%12f,aij);else printf(%12s,); printf(n);solve(a,b);for(i=1;i=N;i+)printf(x%d=%f,i,bi);printf(n);int LU(double aN+1)/对N阶矩A阵进行LU分解，结果L、U存放在A的相应位置/返回值：-1存在元素为0，1正常结束void solve(double aN+1,double bN+1)/利用分解的LU求x/回代求解，L与U元素均在A矩阵相应位置；b存放常数列，返回时存放方程组的解实验十二：Guass-Sediel迭代算法命名（源程序.cpp及工作区.dsw）：GS类型：必做问题：利用GS法求解方程组（精度为）运行及结果显示：#include iostream.h#include math.h#include stdio.h#define N 3/方程的阶数#define MaxK 100 /最大迭代次数#define EPS 0.5e-4/精度static double aaNN=10,-1,-2,-1,10,-2,-1,-1,5;static double bbN=7.2,8.3,4.2;void main()int i,j; double xN+1; double aN+1N+1,bN+1;int GaussSeidel(double aN+1,double b,double x);for(i=1;i=N;i+) for(j=1;j=N;j+) aij=aai-1j-1; bi=bbi-1; xi=0; if(GaussSeidel(a,b,x)=1) printf(the roots is:); for(i=1;i=N;i+) printf( x%d=%f ,i,xi);printf(n); else printf(nthe G-S method failed!n); int GaussSeidel(double aN+1,double b,double x)/参数说明：a传入系数矩阵（有效元素为），b为方程组右边常数列，x传入迭代初值并返回解向量（有效值为）；返回值：1正常结束；-1存在元素为0；-2为迭代次数越界；实验十三：SOR迭代算法命名（源程序.cpp及工作区.dsw）：sor类型：选做问题：利用SOR法求解方程组（精度为）运行及结果显示（）：运行及结果显示（）：#include math.h#include stdio.h#include conio.h#define N 4#define MaxK 1000#define EPS 0.5e-5static double aaNN=-4,1,1,1,1,-4,1,1,1,1,-4,1,1,1,1,-4;static double bbN=1,1,1,1;void main()int i,j; double xN+1; double aN+1N+1,bN+1;int Sor(double aN+1,double b,double w,double x);for(i=1;i=N;i+) for(j=1;j=N;j+) aij=aai-1j-1; bi=bbi-1; xi=0; if(Sor(a,b,1.3,x)=1) printf(the root