fortran数值计算基础

1、 数值计算基础 目录 实验一 直接法解线性方程组的 . 2 实验二 插值方法 . 11 实验三 数值积分 . 5 实验四 常微分方程的数值解 . 7 实验五 迭代法解线性方程组与非线性方程 . 9 . 实验一 直接法解线性方程组 一、实验目的 掌握全选主元消去法与高斯-塞德尔法解线性方程组。 二、实验内容 分别写出Guass列选主元消去法与追赶法的算法，编写程序上机调试出结果，要求所编程序适用于任何一解线性方程组问题，即能解决这一类问题，而不是某一个问题。实验中以下列数据验证程序的正确性。 1、用Guass列选主元消去法求解方程组 2.52.3?5.1x3.7?1?8.5?.x35.39.61

2、?2?5.3.5x8.11.7?4?3 2、用追赶法求解方程组 x?10?20000?1?x00001?2?2?x0001?203?x001?200?4?x00001?2?5 三、实验仪器设备与材料 主流微型计算机 四、实验原理 1、Guass列选主元消去法 对于AX =B )A|B(A）（ 进行变换为1）、消元过程：将是上三角矩阵。即：A|B，其中aa?b1aa?ba? 1121n1111n12?b?aa?a10ab?22122n222n? ?b0ab0?aaa?nn2nnnnnn1n-1 从k1到 、 列选主元a maxa作为主元。 选取第k列中绝对值最大元素 ikk?i?nb、 换行 .

3、 a?a,j?k?1,?,nijkj bb?ikc、 归一化 a/a?a,j?k?1,?,nkjkjkk ba?b/kkkk d、 消元a?aa?a,i?k?1,?,n;j?k?1,?,nijijkjik n,1,?b,i?k?b?ab?iikikx,x,?,x)|(AB。解出）2、回代过程：由 11nn?b/a?xnnnnn ?12,?1,?,b?ax?x,k?nkkkjj?k?1j2、追赶法 线性方程组为： acfx?1111?xfcba?22222?fcbax?33333? ?fcbax1n?1n?11n?n?1?n?fbxa?nnnn做LU分解为： ?1?11?1?222?33?L?,

4、R ?1?1?n?1?nn 分解公式： . ?a(i?2,3,?,n)?ii?b,b3,?,n)?(i?2, ?i1i?i11i?c?i?(i?1,2,?,n?1)? i?i 则f?Ly?f?LUx?Ax?f ?y?Ux? 回代公式：f?1?y? 1?1 ?y?f?1iii?)n?,2,3,(y?i? i?iy?x?nn ?)1,?,n?1,n?2?xy?(xi?1iiii? 五、实验步骤 塞德尔迭代法公式；1、理解并掌握全选主元消去法与高斯- 塞德尔迭代法的流程图2、画出全选主元消去法与高斯- 语言编写出相应的程序并调试验证通过3、使用C 六、实验报告要求、统一使用武汉科技大学实验报告本书写

5、，实验报告的内容要求有：实验目的、1 实验内容、程序流程图、源程序、运行结果及实验小结六个部分。 、源程序需打印后粘贴在实验报告册内；2 、运行结果以屏幕截图形式保存并打印后粘贴在实验报告册内。3 七、实验注意事项 注意如何定义数据结构以保存矩阵和解以降低算法的复杂性。 八、思考题 若使用全主元消去法，在编程中应如何记录保存对于未知数的调换。 . 实验二 插值方法 一、实验目的 掌握拉格郎日插值法与牛顿插值法构造插值多项式。 二、实验内容 分别写出拉格郎日插值法与牛顿插值法的算法，编写程序上机调试出结果，要求所编程序适用于任何一组插值节点，即能解决这一类问题，而不是某一个问题。实验中以下列数据

6、验证程序的正确性。 已知下列函数表 x i0.56160 0.56280 0.56401 0.56521 y i0.82741 0.82659 0.82577 0.82495 时的函数值。求x=0.5635三、实验仪器设备与材料 主流微型计算机 四、实验原理 已知n个插值节点的函数值，则可由拉格郎日插值公式与牛顿插值公式构造出插值多项式，从而由该插值多项式求出所要求点的函数值。拉格郎日插值公式与牛顿插值公式如下： 1、Lagrange插值公式 n?)xly?(y)y?.?l(x)lL(x)?l(x)y?(xkn011nkn00?kx?xn(x?x)(x?x)?(x?x)(x?x)?(x?x)?

7、jn?01k?11k?)?l(x k(x?x)(x?x)?(x?x)(x?x)?(x?x)x?x0j?jkn1kkkkk?kk?011j?k2、Newton插值公式 N(x)?f(x)?fx,x(x?x)?fx,x,x(x?x)(x?x)110n001020 )x(x?)(x?x)?x?fx,x,?(x?x1nn0?011五、实验步骤 1、理解并掌握拉格郎日插值法与牛顿插值法的公式； 2、画出拉格郎日插值法与牛顿插值法算法的流程图； 3、使用C语言编写出相应的程序并调试验证通过。 六、实验报告要求 1、统一使用武汉科技大学实验报告本书写，实验报告的内容要求有：实验目的、 . 实验内容、程序流程

8、图、源程序、运行结果及实验小结六个部分。 2、源程序需打印后粘贴在实验报告册内； 3、运行结果以屏幕截图形式保存并打印后粘贴在实验报告册内。 七、实验注意事项 Newton插值法在编程时应注意定义何种数据结构以保存差商。 八、思考题 比较Lagrange插值法与Newton插值法的异同。 . 实验三 数值积分 一、实验目的 掌握复化梯形法与龙贝格法计算定积分。 二、实验内容 分别写出变步长梯形法与Romberge法计算定积分的算法，编写程序上机调试出结果，要求所编程序适用于任何类型的定积分，即能解决这一类问题，而不是某一个问题。实验中以下列数据验证程序的正确性。 sinx1 ?000010dx

9、,.?。求 x0三、实验仪器设备与材料 主流微型计算机 四、实验原理 通过变步长梯形法与龙贝格法，我们只要知道已知n个求积节点的函数值，则可由相应的公式求出该函数的积分值，从而不需要求该函数的原函数。变步长梯形法与龙贝格法公式如下： 1、变步长梯形法 n?1h?)f(x?xf()?T 1iin?20i?1n?h?)(b)f(x?f?f(a)?2 i2 1?in?1h1?f(?x)?TT 2/in2?1n22 0?i? ?TT?来控制精度 用n2n2、龙贝格法 n?1h1?f(xT?T?) 21nn2/i?22 0i?41T?TS? nnn233 116SS?C? n2nn1515 164C?R

10、?C nnn26363 . ? ?R?R 用来控制精度n2n 五、实验步骤 、理解并掌握变步长梯形法与龙贝格法的公式；1 、画出变步长梯形法与龙贝格法的流程图2 语言编写出相应的程序并调试验证通过3、使用C 六、实验报告要求、统一使用武汉科技大学实验报告本书写，实验报告的内容要求有：实验目的、1 实验内容、程序流程图、源程序、运行结果及实验小结六个部分。 2、源程序需打印后粘贴在实验报告册内； 3、运行结果以屏幕截图形式保存并打印后粘贴在实验报告册内。 七、实验注意事项xsin1?dx，对该点在程序设计中应注意对其1积分中，被积函数在在x=0点函数值为 x0 的定义。 八、思考题 法来计算该问

11、题有何缺点？使用复化梯形法与复化Simpson . 实验四 常微分方程的数值解 一、实验目的 掌握改进欧拉法与四阶龙格-库塔求解一阶常微分方程的初值问题。 二、实验内容 分别写出改进欧拉法与四阶龙格-库塔求解的算法，编写程序上机调试出结果，要求所编程序适用于任何一阶常微分方程的数值解问题，即能解决这一类问题，而不是某一个问题。实验中以下列数据验证程序的正确性。 2?xy?y求步长h=0.25。 ?)?x?5y(0)?2(0?三、实验仪器设备与材料 主流微型计算机 四、实验原理 常微分方程的数值解主要采用“步进式”，即求解过程顺着节点排列次序一步一步向前推进，在单步法中改进欧拉法和四阶龙格-库塔

12、法公式如下： 1、改进欧拉法 y?y?hf(x,y) n?1nnnh y?y?f(x,y)?f(x,y) 11n?1nnnn?n2 2、四阶龙格-库塔法 h?y?y?(k?2k?2k?k)? n1?3n2416?k?f(x,y)?n1n?hh?k?f(x?,y?k) ? 2n1n22?hh?k?f(x?,y?k)? 3n2n22?k?f(x?h,y?hk)?34nn五、实验步骤 1、理解并掌握改进欧拉法与四阶龙格-库塔法的公式； 2、画出改进欧拉法与四阶龙格-库塔法的流程图 3、使用C语言编写出相应的程序并调试验证通过 六、实验报告要求 . 1、统一使用武汉科技大学实验报告本书写，实验报告的内

13、容要求有：实验目的、实验内容、程序流程图、源程序、运行结果及实验小结六个部分。 2、源程序需打印后粘贴在实验报告册内； 3、运行结果以屏幕截图形式保存并打印后粘贴在实验报告册内。 七、实验注意事项 2?xy?y? 2)x1?y?2/(，通过调整步长，观察结果的的精确解为?)x?5(0)?20?y(? 精度的变化八、思考题 如何对四阶龙格-库塔法进行改进，以保证结果的精度。 . 实验五 迭代法解线性方程组与非线性方程 一、实验目的 掌握高斯-塞德尔迭代法求解线性方程组与牛顿迭代法求方程根。 二、实验内容 分别写出高斯-塞德尔迭代法与牛顿迭代法的算法，编写程序上机调试出结果，要求所编程序适用于任何

14、一个方程的求根，即能解决这一类问题，而不是某一个问题。实验中以下列数据验证程序的正确性。 1、高斯-塞德尔迭代法求解线性方程组 x47?221?1?x7?23915?2? ?x15?2?2113?x013132?4?3?0.0000101x?x?，牛顿法的初始值为2、用牛顿迭代法求方程的近似根，1。 三、实验仪器设备与材料 主流微型计算机 四、实验原理 二分法通过将含根区间逐步二分，从而将根的区间缩小到容许误差范围。牛顿通过迭代的方法逐步趋进于精确解，该两种方法的公式如下： 1、高斯-塞德尔迭代法 1）判断线性方程组是否主对角占优 n? ?a,i?1,2,?,an iiij1?jij? 2）直

15、接分离xi，即 n?n,?1,2d?)bx/a,i(x? iiijiij1j?建立高斯-塞德尔迭代格式为： i?1n?()k)k(?1(k?1)/a,i?1?a,xx2,?(?d?ax,n iiijijiijj1?i?j?j1 ）取初值迭代求解至所要求的精度为止。32、牛顿法 . f(x)k?xx? k?k1?)xf(k 五、实验步骤 1、理解并掌握二分法与牛顿法的公式； 2、画出二分法与牛顿法的流程图 C语言编写出相应的程序并调试验证通过3、使用 六、实验报告要求、统一使用武汉科技大学实验报告本书写，实验报告的内容要求有：实验目的、1 实验内容、程序流程图、源程序、运行结果及实验小结六个部分

16、。 2、源程序需打印后粘贴在实验报告册内； 3、运行结果以屏幕截图形式保存并打印后粘贴在实验报告册内。 七、实验注意事项对于二分法应注意二分后如何判断根的区间，对于牛顿法注意如何确定迭代过程的结 束 八、思考题 若使用牛顿法是发散的，如何对牛顿法进行改进以保证其收敛性。 ）和运行结果截图前三个实验的程序代码（C/C+ 全选主元解方程组的源程序及运行结果Gauss#include #include #include using namespace std; class Matrix public: Matrix(); Matrix(); n构造线性方程组相应的矩阵， void SetMatrix

17、(const int n,const double esp1);/ 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=in; if(nesp1; if(esp10) e

18、sp1=0;/输入不合法的精度就把精度置0 while(esp1=0); 潣瑵?输入线性方程组的增广矩阵:n; matrix.SetMatrix(n,esp1);/设置矩阵内的数据 matrix.Calculate();/计算方程组的解 /输出方程组的解 . coutn方程组的解如下: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; locat

19、ion.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) 潣瑵?分配存储区失败！n; exit(EXIT_FAILURE); /读取线性方程组的增广矩阵 . f

20、or(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; locatio

21、n.column=j; /最大主元小于输入的精度时，认为方程组无解，退出程序 if(max=esp) 潣瑵?方程组无解！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;

22、 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 Matri

23、x:Eliminate(const int r) if(fabs(*(a+N*r+r)=esp) 潣瑵?方程组无解！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()const . if(fabs(*(a+N*(N-1)+N-1)=0;-i) temp=0; for(int j=i+1;jN;+

24、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() /根据矩阵行

25、数重复进行寻找最大主元、变换行或列、消元 for(int i=0;iGetRank()-1;+i) Max(i); ChangeRC(i); Eliminate(i); Result(); . int Matrix:GetRank()const return N;/返回矩阵的行数 double Matrix:GetX(const int i)const return *(b+i-1); 运行结果 追赶法求解方程组的算法：，主对角元素左边的元素b(i)(i=0:n-1)，输入方程组的维数n，将主对角元素 1f(i)(i=0:n-1) c(i)(i=0:n-2)，右端项的元素a(i)(i=0:n-

26、2)，主对角元素右边的元素i=0:n-2, 对于, (0)=b(0),Crout2 对方程组的系数矩阵作分解(i) (i):= c(i)/b(i), b(i+1):=(i+1):= b(i+1)-a(i)* c(i):= Ly=f 3.解方程组(0):=f(0)/b(0) b(0):=y(0):=f(0)/ b(i):=y(i):=f(i)-a(i-1)*y(i-1)/b(i) ， 对于i=1：n-1 Ux=y 4.解方程组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) 用追赶法求

27、解方程组的源程序及运行结果 #include #include using namespace std; class MatrixThr public: 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为主对角线左边一斜条元素的首地址

28、， /c为主对角线右边一斜条元素首地址，f为方程组的常数首地址 double *a,*b,*c,*f; ; int main() MatrixThr matrix; int n; do 潣瑵?输入三对角方程组的变量的数目N:; cinn; while(n3); 潣瑵?请依次输入三对角方程组每行的数据(0元素除外):n; matrix.SetMatrixThr(n); /计算方程组的解 matrix.Result();/输出方程组的解 潣瑵?方程的解如下:n; for(int i=1;i=matrix.GetN();+i) coutXi:matrix.GetX(i)*b*c*f; for(int

29、 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()const r

30、eturn N; double MatrixThr:GetX(const int i)const return *(a+i-1); 运行结果 Lagrange插值法的源程序：#include . #include using namespace std; class Lagrange public: Lagrange(); Lagrange(); void SetLagrange(const int n);/根据用户的输入设置Lagrange类中的插值点数据 bool Exist(const double x,const int i);/检测是否输入了与前i个插值结点横坐标相同的点 int G

31、etN()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; do 晩渨监?潣瑵?用于模拟曲线的插值点数目N2n; while(na; L.Calculate(a); 潣瑵?横坐标a-对应的函数值为:L.GetResult()endl; ret

32、urn 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

33、doubleN; y=new doubleN; /判断是否成功为指针分配了内存空间 if(x=NULL|y=NULL) 潣瑵?分配存储空间失败！endl; exit(EXIT_FAILURE); /输入插值点 for(int i=0;i*(x+i)*(y+i); /如果不是输入第一个坐标值，则会对输入的横坐标进行合法性检测 while(i!=0&Exist(*(x+i),i)=true) 潣瑵?输入了重复的横坐标-请重新输入第?椼?个插入结点的数据n; cin*(x+i)*(y+i); . void Lagrange:Calculate(const double a) zx=a; for(in

34、t 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()const return N;/返回插值点的数目 double Lagrange:GetResult() const return zy;/返回求得的函数值 Lagrange插值法的运行结果 . Newton插值法的源程序： #in

35、clude using namespace std; class Matrix public: 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 潣瑵?输入插值点的数目N：; cinn; while(n2); 潣瑵?输入插值点的数据：x