计算方法上机报告.docxVIP

计算方法上机作业计算方法上机作业共轭梯度法方法描述共轭梯度法是把求解线性方程组的问题转化为求解一个与之等价的二次泛函极小值问题，即： 其中，泛函共轭梯度法在形式上迭代法的特征只要确定了初始向量，然后确定确定最佳步长和搜索方向就可以得出迭代格式。给定初始向量后，由于负梯度方向是下降最快的方向，所以第一次迭代的搜索方向，得；第二次迭代时，选取的搜索方向不再是 ，而是使得与是关于矩阵K的共轭向量，即 ，由此可以求得参数，然后沿方向进行搜索得。最佳步长的计算是在给定迭代点以及搜索方向后，要求选取非负实数使得达到最小，即对求导，令，得：所以由此确定了每一步的最佳步长。实际上的计算还可以化简，在此不再推导，化简后的表达式在下列算法中。程序流程程序分为两个部分，主程序读取或产生矩阵和右端向量A,b以及初始值x0和误差限和最大迭代次数tol、Imax，然后调用方法部分进行计算得出x，方法部分放在iter_solver模块中，调用函数接口为call CG_method(A,b,x,x0,tol,Imax)。程序流程如下： 子程序CG_method(A,b,x,x0,tol,Imax)的流程如下：给定最大迭代次数以及误差限：Imax， tol给定初值：CG: Do i=1,Imax IF()EXIT END DO CG 程序编制实际上利用Fortran语言编写了一个迭代法求解线性方程组的模块module_iter_solver，其中包含Jacobi迭代法，Gauss-Seidel迭代法，Successive over-relaxation迭代法，以及Conjugate gradient method。每个子程序都有验证，其中有两个主程序验证共轭梯度法，并在Linux系统下Gfortran编译器以及Windows系统下InterFortran编译器通过验证。调用共轭梯度法的接口为CALL CG_method(A,b,x,x0,tol,Imax)，子程序CG_method(A,b,x,x0,tol,Imax)中输入值A,b,x0,tol,Imax分别为线性方程组左端矩阵、右端向量、迭代初始值、误差限、最大迭代次数，返回值为解x；实际子程序中将迭代中间值记录到一个专门的文件中。程序验证编写了两个主程序验证Module iter_solver中的CG_method(A,b,x,x0,tol,Imax)，分别为在文件中读取A,b,x0,tol,Imax然后调用CG_method(A,b,x,x0,tol,Imax)和程序中产生A,b,x0,tol,Imax然后调用CG_method(A,b,x,x0,tol,Imax)。验证一主程序从in.dat文件中读取tol,Imax ,N,A,b其中N为A的维数，输入文件格式如下结果为中间结果为从中间结果可以看出只迭代两次就求出真解。验证二主程序内部产生A,b,x0,tol,Imax然后调用CG_method(A,b,x,x0,tol,Imax)；产生的100阶的A,b分别为给定的x0=0. ;tol=1.0E-8 ;Imax=10000 ；结果为：中间结果为：从中间结果可以看出迭代50次得出真解。由以上结果可以验证module_iter_solver模块中的CG_method是可靠的，实际上其它子程序也进行了相关验证，这里没有给出验证结果。所以module_iter_solver是通用的，可以用其中的各种迭代解法来求解线性方程组。其它子程序的调用方法参考程序说明。三次样条插值方法描述三次样条插值函数定义设在区间[a,b]上给定n+1个节点,在节点处的函数值，若函数满足以下三条：在每个子区间上，是三次多项式；；在区间[a,b]上，的二阶导函数连续，则称函数为函数在区间[a,b]上的三次样条插值函数。由定义可知由这4n-2个条件加上两个边界条件可求出4n个待定参数。三次样条插值函数设，则其中对式两端积分两次得到由插值条件得将带入可以得出在区间上的表达式为：关于参数的方程组对求导得：由得该式可以化简为其中第一种边界条件为已知，取用矩阵表示为其中第二种边界条件为已知即，将其带入可以得到其中用矩阵表示为：其中第三种边界条件为周期型边界条件此时有，i=n时式也成立，即其中用矩阵表示为其中程序流程程序分为三个部分，包括求解线性方程组的equation_solver模块，三次样条插值模块cub_spline，和主程序main；主程序读取或产生插值点(x,y)、区间数n、边界条件boundary以及要计算的函数点值得个数，然后调用cub_spline模块中的cal_function_parameter(M,x,y,n,boundary)计算出三次样条函数中的参数M，然后调用cub_spline模块中的cal_functi