fortran数值计算基础.doc

数值计算基础目录实验一 直接法解线性方程组的1实验二 插值方法4实验三 数值积分6实验四 常微分方程的数值解8实验五 迭代法解线性方程组与非线性方程10实验一 直接法解线性方程组一、实验目的掌握全选主元消去法与高斯-塞德尔法解线性方程组。二、实验内容分别写出Guass列选主元消去法与追赶法的算法，编写程序上机调试出结果，要求所编程序适用于任何一解线性方程组问题，即能解决这一类问题，而不是某一个问题。实验中以下列数据验证程序的正确性。1、用Guass列选主元消去法求解方程组2、用追赶法求解方程组三、实验仪器设备与材料主流微型计算机四、实验原理1、Guass列选主元消去法对于AX =B 1）、消元过程：将（A|B）进行变换为，其中是上三角矩阵。即：k从1到n-1a、 列选主元 选取第k列中绝对值最大元素作为主元。b、 换行 c、 归一化 d、 消元 2）、回代过程：由解出。 2、追赶法线性方程组为：做LU分解为： 分解公式：则回代公式：五、实验步骤1、理解并掌握全选主元消去法与高斯-塞德尔迭代法公式；2、画出全选主元消去法与高斯-塞德尔迭代法的流程图3、使用C语言编写出相应的程序并调试验证通过六、实验报告要求1、统一使用武汉科技大学实验报告本书写，实验报告的内容要求有：实验目的、实验内容、程序流程图、源程序、运行结果及实验小结六个部分。2、源程序需打印后粘贴在实验报告册内；3、运行结果以屏幕截图形式保存并打印后粘贴在实验报告册内。七、实验注意事项注意如何定义数据结构以保存矩阵和解以降低算法的复杂性。八、思考题若使用全主元消去法，在编程中应如何记录保存对于未知数的调换。实验二 插值方法一、实验目的掌握拉格郎日插值法与牛顿插值法构造插值多项式。二、实验内容分别写出拉格郎日插值法与牛顿插值法的算法，编写程序上机调试出结果，要求所编程序适用于任何一组插值节点，即能解决这一类问题，而不是某一个问题。实验中以下列数据验证程序的正确性。已知下列函数表xi0.561600.562800.564010.56521yi 0.827410.826590.825770.82495求x=0.5635时的函数值。三、实验仪器设备与材料主流微型计算机四、实验原理已知n个插值节点的函数值，则可由拉格郎日插值公式与牛顿插值公式构造出插值多项式，从而由该插值多项式求出所要求点的函数值。拉格郎日插值公式与牛顿插值公式如下：1、Lagrange插值公式2、Newton插值公式五、实验步骤1、理解并掌握拉格郎日插值法与牛顿插值法的公式；2、画出拉格郎日插值法与牛顿插值法算法的流程图；3、使用C语言编写出相应的程序并调试验证通过。六、实验报告要求1、统一使用武汉科技大学实验报告本书写，实验报告的内容要求有：实验目的、实验内容、程序流程图、源程序、运行结果及实验小结六个部分。2、源程序需打印后粘贴在实验报告册内；3、运行结果以屏幕截图形式保存并打印后粘贴在实验报告册内。七、实验注意事项Newton插值法在编程时应注意定义何种数据结构以保存差商。八、思考题比较Lagrange插值法与Newton插值法的异同。实验三 数值积分一、实验目的掌握复化梯形法与龙贝格法计算定积分。二、实验内容分别写出变步长梯形法与Romberge法计算定积分的算法，编写程序上机调试出结果，要求所编程序适用于任何类型的定积分，即能解决这一类问题，而不是某一个问题。实验中以下列数据验证程序的正确性。求。三、实验仪器设备与材料主流微型计算机四、实验原理通过变步长梯形法与龙贝格法，我们只要知道已知n个求积节点的函数值，则可由相应的公式求出该函数的积分值，从而不需要求该函数的原函数。变步长梯形法与龙贝格法公式如下： 1、变步长梯形法用来控制精度2、龙贝格法用来控制精度五、实验步骤1、理解并掌握变步长梯形法与龙贝格法的公式；2、画出变步长梯形法与龙贝格法的流程图3、使用C语言编写出相应的程序并调试验证通过六、实验报告要求1、统一使用武汉科技大学实验报告本书写，实验报告的内容要求有：实验目的、实验内容、程序流程图、源程序、运行结果及实验小结六个部分。2、源程序需打印后粘贴在实验报告册内；3、运行结果以屏幕截图形式保存并打印后粘贴在实验报告册内。七、实验注意事项在积分中，被积函数在x=0点函数值为1，对该点在程序设计中应注意对其的定义。八、思考题使用复化梯形法与复化Simpson法来计算该问题有何缺点？实验四 常微分方程的数值解一、实验目的掌握改进欧拉法与四阶龙格-库塔求解一阶常微分方程的初值问题。二、实验内容分别写出改进欧拉法与四阶龙格-库塔求解的算法，编写程序上机调试出结果，要求所编程序适用于任何一阶常微分方程的数值解问题，即能解决这一类问题，而不是某一个问题。实验中以下列数据验证程序的正确性。求步长h=0.25。三、实验仪器设备与材料主流微型计算机四、实验原理常微分方程的数值解主要采用“步进式”，即求解过程顺着节点排列次序一步一步向前推进，在单步法中改进欧拉法和四阶龙格-库塔法公式如下： 1、改进欧拉法2、四阶龙格-库塔法五、实验步骤1、理解并掌握改进欧拉法与四阶龙格-库塔法的公式；2、画出改进欧拉法与四阶龙格-库塔法的流程图3、使用C语言编写出相应的程序并调试验证通过六、实验报告要求1、统一使用武汉科技大学实验报告本书写，实验报告的内容要求有：实验目的、实验内容、程序流程图、源程序、运行结果及实验小结六个部分。2、源程序需打印后粘贴在实验报告册内；3、运行结果以屏幕截图形式保存并打印后粘贴在实验报告册内。七、实验注意事项 的精确解为，通过调整步长，观察结果的精度的变化八、思考题如何对四阶龙格-库塔法进行改进，以保证结果的精度。实验五 迭代法解线性方程组与非线性方程一、实验目的掌握高斯-塞德尔迭代法求解线性方程组与牛顿迭代法求方程根。二、实验内容分别写出高斯-塞德尔迭代法与牛顿迭代法的算法，编写程序上机调试出结果，要求所编程序适用于任何一个方程的求根，即能解决这一类问题，而不是某一个问题。实验中以下列数据验证程序的正确性。1、高斯-塞德尔迭代法求解线性方程组2、用牛顿迭代法求方程的近似根，牛顿法的初始值为1。三、实验仪器设备与材料主流微型计算机四、实验原理二分法通过将含根区间逐步二分，从而将根的区间缩小到容许误差范围。牛顿通过迭代的方法逐步趋进于精确解，该两种方法的公式如下： 1、高斯-塞德尔迭代法1）判断线性方程组是否主对角占优 2）直接分离xi，即建立高斯-塞德尔迭代格式为：3）取初值迭代求解至所要求的精度为止。2、牛顿法五、实验步骤1、理解并掌握二分法与牛顿法的公式；2、画出二分法与牛顿法的流程图3、使用C语言编写出相应的程序并调试验证通过六、实验报告要求1、统一使用武汉科技大学实验报告本书写，实验报告的内容要求有：实验目的、实验内容、程序流程图、源程序、运行结果及实验小结六个部分。2、源程序需打印后粘贴在实验报告册内；3、运行结果以屏幕截图形式保存并打印后粘贴在实验报告册内。七、实验注意事项对于二分法应注意二分后如何判断根的区间，对于牛顿法注意如何确定迭代过程的结束八、思考题若使用牛顿法是发散的，如何对牛顿法进行改进以保证其收敛性。前三个实验的程序代码（C/C+）和运行结果截图Gauss全选主元解方程组的源程序及运行结果#include #include #include using namespace std;class Matrixpublic:Matrix();Matrix();void SetMatrix(const int n,const double esp1);/构造线性方程组相应的矩阵，n为方程的未知数数目，esp1为要求的精度void Max(const int r);/全选主元void ChangeRC(const int r);/根据主元变换矩阵的行或列void Eliminate(const int r);/处理消元工作void Result()const;/计算方程的解void Calculate();int GetRank()const;/返回矩阵的行数double GetX(const i)const;/确定方程组的第i个解(1=i=N)private:/指针a和b分别用于存储方程组的未知数系数和方程=右边的常数，esp存/储精度double *a,*b,esp;/指针flag用于记录方程组解的顺序int *flag;/以下的结构体用于在全选主元中记录最大主元的位置struct coordinateint row,column;location;int N;/方程组未知数的数目;int main()int n;double esp1;Matrix matrix;docoutn;if(n0) n=0;/如何控制非法字符的输入？while(n=0);docoutesp1;if(esp10) esp1=0;/输入不合法的精度就把精度置0while(esp1=0); cout输入线性方程组的增广矩阵:n;matrix.SetMatrix(n,esp1);/设置矩阵内的数据matrix.Calculate();/计算方程组的解/输出方程组的解coutnn方程组的解如下:n;for(int i=1;i=matrix.GetRank();+i)coutXi:matrix.GetX(i)endl;return 0;Matrix:Matrix()/将Matrix类的数据成员初始化a=NULL;b=NULL;flag=NULL;location.row=0;location.column=0;esp=0;N=0;Matrix:Matrix()/释放指针a、b和flag指向的内存空间deletea;deleteb;deleteflag;void Matrix:SetMatrix(const int n,const double esp1)N=n;esp=esp1;a=new doubleN*N;b=new doubleN;flag=new intN;/判断是否成功分配存储区if(a=NULL|b=NULL|flag=NULL)cout分配存储区失败！n;exit(EXIT_FAILURE); /读取线性方程组的增广矩阵for(int i=0;iN;+i)for(int j=0;j*(a+i*N+j);cin*(b+i); /flag中存储的值对应相应的x值，当方程的解由于列变换交换后，flag中/的值也相应交换，最后用于恢复解的顺序for( i=0;iN;+i) *(flag+i)=i;void Matrix:Max(const int r)double max=0;for(int i=r;iN;+i)for(int j=r;jN;+j)if(maxfabs(*(a+i*N+j)max=fabs(*(a+i*N+j);/设定最大主元的行、列location.row=i;location.column=j; /最大主元小于输入的精度时，认为方程组无解，退出程序if(max=esp)cout方程组无解！n;exit(EXIT_FAILURE);void Matrix:ChangeRC(const int r)double temp;/如果最大主元所在的行不在当前行，则进行行变换if(location.row!=r)for(int i=r;iN;+i)temp=*(a+r*N+i);*(a+r*N+i)=*(a+location.row*N+i);*(a+location.row*N+i)=temp;temp=br;br=blocation.row;blocation.row=temp;/若最大主元所在的列不在当前的r列，则进行列变换if(location.column!=r)for(int i=r;iN;+i)temp=*(a+i*N+r);*(a+i*N+r)=*(a+i*N+location.column);*(a+i*N+location.column)=temp;/交换flag中的元素来标记方程解的位置变化int temp1;temp1=*(flag+r);*(flag+r)=*(flag+location.column);*(flag+location.column)=temp1;void Matrix:Eliminate(const int r)if(fabs(*(a+N*r+r)=esp)cout方程组无解！n;exit(EXIT_FAILURE);for(int i=r+1;iN;+i) for(int j=r+1;jN;+j) (*(a+i*N+j)-=(*(a+i*N+r)*(*(a+r*N+j)/(*(a+r*N+r); (*(b+i)-=(*(b+r)*(*(a+i*N+r)/(*(a+r*N+r);void Matrix:Result()constif(fabs(*(a+N*(N-1)+N-1)=esp)cout=0;-i)temp=0;for(int j=i+1;jN;+j) temp+=(*(a+i*N+j)*(*(b+j);*(b+i)=(*(b+i)-temp)/(*(a+i*N+i); /根据flag中的数据用冒泡排序法恢复方程组解的次序for(i=0;iN-1;+i)for(int j=0;j*(flag+j+1)int temp1; /交换解的顺序temp=*(b+j);*(b+j)=*(b+j+1);*(b+j+1)=temp; /交换用于标记的元素的顺序temp1=*(flag+j);*(flag+j)=*(flag+j+1);*(flag+j+1)=temp1;void Matrix:Calculate()/根据矩阵行数重复进行寻找最大主元、变换行或列、消元for(int i=0;iGetRank()-1;+i)Max(i); ChangeRC(i);Eliminate(i);Result();int Matrix:GetRank()constreturn N;/返回矩阵的行数double Matrix:GetX(const int i)constreturn *(b+i-1);运行结果追赶法求解方程组的算法：1 输入方程组的维数n，将主对角元素b(i)(i=0:n-1)，主对角元素左边的元素a(i)(i=0:n-2)，主对角元素右边的元素c(i)(i=0:n-2)，右端项的元素f(i)(i=0:n-1)2 对方程组的系数矩阵作Crout分解, (0)=b(0),对于i=0:n-2,c(i):=(i):= c(i)/b(i), b(i+1):=(i+1):= b(i+1)-a(i)* (i)3.解方程组Ly=f b(0):=y(0):=f(0)/ (0):=f(0)/b(0) 对于i=1：n-1， b(i):=y(i):=f(i)-a(i-1)*y(i-1)/b(i)4.解方程组Ux=y a(n-1):=x(n-1):=b(n-1)对于i=n-2:0,a(i)=x(i):=b(i)-c(i)*a(i+1);5输出方程组的解a(0:n-1)用追赶法求解方程组的源程序及运行结果#include#includeusing namespace std;class MatrixThrpublic:MatrixThr();MatrixThr();void SetMatrixThr(const int n);/设置三对角矩阵的数据void Result();/计算三对角矩阵的解double GetX(const int i)const;/取得第i个解,i从1开始int GetN() const;/返回未知数的数目private:int N;/N为未知数的数目 /b为矩阵主对角线的元素首地址，a为主对角线左边一斜条元素的首地址，/c为主对角线右边一斜条元素首地址，f为方程组的常数首地址double *a,*b,*c,*f;int main()MatrixThr matrix;int n;docoutn;while(n3);cout请依次输入三对角方程组每行的数据(0元素除外):n;matrix.SetMatrixThr(n); /计算方程组的解 matrix.Result();/输出方程组的解cout方程的解如下:n; for(int i=1;i=matrix.GetN();+i) coutXi:matrix.GetX(i)endl;return 0;MatrixThr:MatrixThr()/初始化相关数据N=0;a=NULL;b=NULL;c=NULL;f=NULL;MatrixThr:MatrixThr() /释放分配的内存空间delete a;delete b;delete c;delete f;void MatrixThr:SetMatrixThr(const int n) /根据输入的未知数个数设置矩阵的数据N=n;a=new doubleN;b=new doubleN;c=new doubleN;f=new doubleN; /若内存分配失败，退出程序if(a=NULL|b=NULL|c=NULL|f=NULL) cout*b*c*f;for(int i=1;i*(a+i-1)*(b+i)*(c+i)*(f+i);cin*(a+i-1)*(b+i)*(f+i);void MatrixThr:Result() /对系数矩阵A作Crout分解for(int i=0;iN-1;+i) /将U中的存于指针b中，L中的存于指针c中*(c+i)/=(*(b+i);*(b+i+1)-=(*(a+i)*(*(c+i); /解方程组Ly=f，求得的y值存于指针b中*b=(*f)/(*b);for(i=1;i=0;-i) *(a+i)=*(b+i)-(*(c+i)*(*(a+i+1);int MatrixThr:GetN()constreturn N;double MatrixThr:GetX(const int i)constreturn *(a+i-1);运行结果Lagrange插值法的源程序：#include#includeusing namespace std;class Lagrangepublic:Lagrange();Lagrange();void SetLagrange(const int n);/根据用户的输入设置Lagrange类中的插值点数据bool Exist(const double x,const int i);/检测是否输入了与前i个插值结点横坐标相同的点int GetN()const;/获取插值结点的数目void Calculate(const double a);/计算横坐标a对应的函数值double GetResult()const;/返回计算的函数值private:int N;/插值结点的数目double *x,*y,zx,zy;/x、y分别用于存储插值点的数据，zx、zy表示所求的坐标点;int main()int n=2;Lagrange L;doif(n2) cout2endl;coutn;while(n2);cout请输入各插值点横、纵坐标的数据n;L.SetLagrange(n);double a;couta;L.Calculate(a);cout横坐标a-对应的函数值为:L.GetResult()endl;return 0;Lagrange:Lagrange() /初始化相关数据N=0;x=y=NULL;zx=zy=0;Lagrange:Lagrange() /释放分配给指针的内存空间delete x;delete y;bool Lagrange:Exist(const double a,const int i) /遍历以输入的插值点，看是否重复插入横坐标相同的插值点for(int j=0;ji;+j) if(a=*(x+j) return true;return false;void Lagrange:SetLagrange(const int n)N=n;x=new doubleN;y=new doubleN;/判断是否成功为指针分配了内存空间if(x=NULL|y=NULL) cout分配存储空间失败！endl;exit(EXIT_FAILURE); /输入插值点for(int i=0;i*(x+i)*(y+i);/如果不是输入第一个坐标值，则会对输入的横坐标进行合法性检测while(i!=0&Exist(*(x+i),i)=true)cout输入了重复的横坐标-请重新输入第i+1*(x+i)*(y+i);void Lagrange:Calculate(const double a)zx=a;for(int i=0;iGetN();+i) /计算插值基函数double temp=1;for(int j=0;ji;+j) temp*=(zx-*(x+j)/(*(x+i)-*(x+j);for(+j;jGetN();+j) temp*=(zx-*(x+j)/(*(x+i)-*(x+j);zy+=(*(y+i)*temp);int Lagrange:GetN()constreturn N;/返回插值点的数目double Lagrange:GetResult() constreturn zy;/返回求得的函数值Lagrange插值法的运行结果Newton插值法的源程序：#include using namespace std;class Matrixpublic:Matrix();Matrix();void SetMatrix(const int n);/根据用户输入的插值点的数据设置计算结果的矩阵int GetN()const;/返回插值点的数目void Calculate();/计算差商double GetResult(double x)const;/根据输入的横坐标求函数值，返回运算结果private:double *a,*f;/a和f分别用于存储插值点的横坐标和相应的函数值int N;/记录插值点的数目;int main()Matrix matrix;int n;do coutn;while(n2);cout输入插值点的数据：endl;matrix.SetMatrix(n);matrix.Calculate(); double x;coutx;cout所求的函数值为：matrix.GetResult(x)endl;return 1;Matrix:Matrix() /初始化数据a=f=NULL;N=0;Matrix:Matrix()/释放指针a、f指向的内存空间delete a;delete f;void Matrix:SetMatrix(int n)N=n;/为插值点创建动态数组a=new doubleN;f=new doubleN;/输入插值点的数据for(int i=0;i*(a+i)*(f+i);int Matrix:GetN() constreturn N;void Matrix:Calculate() /将差商存储在一个一维数组内for(int i=0;ii;-j)*(f+j)=(*(f+j)-*(f+j-1)/(*(a+j)-*(a+j-i-1);double Matrix:GetResult(double x) const /利用差商和插值点的横坐标及第一个插值点的纵坐标计算函数值double result=*f;/指针f指向第一个插值点的纵坐标for(int i=1;iGetN();+i) double temp=1;for(int j=0;ji;+j)temp*=(x-*(a+j);/计算(x-x0)*(x-x1)*(x-x(i-1)result+=(*(f+i)*temp;/指针(f+i)指向第i号差商return result;Newton插值法的运行结果变步长梯形法求定积分的源程序：#include#includeusing namespace std;double function(const double x);/求被积函数的值并返回/accumulate()为求定积分的函数，a、b分别为积分的上下限,默认精度为0.00001double accumulate(const double a,const double b,const double eps=0.00001);int main()double a,b,eps;/a,b分别为定积分的上限和下限,h为步长,eps为要求的精度a=0;b=1;eps=0.00001;cout(sinX/X)在0到1上的积分为:accumulate(a,b)eps)n*=2;h/=2;/步长折半T1=T2; /利用T1计算T2double temp=0;for(int i=1;i=n;+i) temp+=function(a+(i-1.0/2)*h);T2=T1/2+h*temp/2;return T2;/返回步长h最合适时定积分的近似值变步长梯形法求定积分的运行结果Romberg算法求定积分的源程序：#include #include using namespace std;/求被积函数的值并返回double fun(const double x);/Romberge()为求定积分的函数，a、b分别为积分的上下限,默认精度为0.00001double Romberg(const double a,const double b,const double eps=0.0001);/函数Tm()为T-数表的计算公式double Tm(const double T1,const double T2,const int m);int main()double a,b,eps;/a,b分别为定积分的上限和下限,h为步长,eps为要求的精度a=0;b=1;eps=0.00001;cout所求积分为:Romberg(a,b,eps)endl;return 0;double fun(const double x)if(x=0) return 1;/x为0时函数值为1return(sin(x)/x);/返回被积函数的值double Tm(const double T1,const double T2,const int m) /根据T1、T2和m计算出下一个数并返回return(pow(4,m)*T1-T2)/(pow(4,m)-1);double Romberg(const double a,const double b,const double eps)double T10;/用于存储T-数表某一行的数据int l=1;T0=(b-a)*(fun(a)+fun(b)/2;double T1,T2;int flag;/flag用于标记十分求得符合精度的值doflag=0;int n=pow(2,l-1);/n表示区间等分数double h=(b-a)/n;/h为步长double temp=0;for(int i=0;in;+i) temp+=fun(a+(i+1.0/2)*h);T1=T0/2+h*temp/2;/利用变步长梯形公式计算T-数表第（l+1）行的第一个数据/将T-数表第l+1行的前l-1个数据存入数组for(int m=1;meps)/比较T-数表第（l+1）行和第i行的最后一个数据的差值 /将T-数表第l+1行的最后两个个数据存入数组Tl-1=T1;Tl=T2;+l;flag=1;/表示没有求得符合精度的值while(flag);return T2;/返回符合精度的解Romberg算法求定积分运行结果数值计算基础考试样卷一、单项选择题（每小题3分，共15分）1、数值x的近似值x*=0.1215102，若满足( )，则称x有4位有效数字. (A) 103 (B) 104 (C) 105 (D) 1062、若为矩阵A的k阶主子矩阵，则矩阵A满足（ ）时，则存在唯一单位下三角阵和上三角阵，使。(A) (B) 某个(C) (D) 3、通过四个互异节点的插值多项式P(x)，只要满足( ), 则P(x)是不超过一次多项式。(A) 初始值y0=0 (B) 所有一阶均差为0 (C) 所有二阶均差为0 (D) 所有三阶均差为04、牛顿切线法求解方程f(x)=0的近似根，若初始值x0满足( )，则解的迭代数列一定收敛。(A)0 (C)0 (D)05、改进欧拉法的平均形式公式是( )(A) (B) (C) (D)二、填空题（每小题3分，共15分）1、sin1有2位有效数字的近似值0.84的相对误差限是 . 2、设f(x)可导，求方程x=f(x) 根的牛顿迭代格式是 .3、设，则 . 4、在区间上的插值型求积公式系数满足 .5、二阶龙格库塔法的局部截断误差是 .三、解答题（每小题10分，共50分）1、用列主元消去法解线性方程组2、用牛顿法求的近似值，取初始值，进行二次迭代。3、已知有y=f(x)的函数表如下x123y137求其代数插值多项式并给出其余项。4、给出数值积分公式：确定A、B使得该数值积分公式的代数精度尽可能的高，并确定其代数精度为多少？5、用欧拉法解初值问题，要求保留4位有效数字。四、综合题（每小题10分，共20分）1、试利用数值积分的方法推导求解初值问题的梯形公式为 ，并证明该方法是二阶方法。2、设l0(x)是以n+1个互异点x0,x1,x2,xn 为节点的拉格朗日插值基函数 试利用牛顿插值法证明： 数值计算基础考试样卷参考答案一、单项选择题（每小题3分，共15分）1、D 2、D 3、C 4、B 5、D二、填空题（每小题3分，共15分）1、2、3、 64、 b-a5、 O(h3)三、解答题（每小题10分，共50分）1、解： 8分回代得 2分2、解： 7分 3分3、解法一： 待定系数法设，则 （3分） （3分）即 （1分）法二：Lagrange插值法法三：Newton插值法xiyi一阶差商二阶差商1123237441 （3分） （4分）余项为 （3