1、 计算方法（ A）大作业 姓名： 班级： 专业： 学号： 1 共轭梯度法 一、 算法原理 共轭梯度法是把求解线性方程组的问题转化为求解一个与之等价的二次函数极小化的问题，因此从任意给定的初始点出发，沿一组关于矩阵 A 的共轭方向进行线性搜索，在无舍入无差的假定下，最多迭代 n(其中 n 为矩阵 A 的阶数）次就可求得二次函数的极小点，也就求得了线性方程组 Ax=B 的解。 下述定理给出了求系数矩阵 A 是对称正定矩阵的线性方程组 Ax=b 的解与求二次函数 () = 12 极小点的等价性。 定理 3.4.1 设 A 是 n 阶对称正定矩阵，则 是方程组 Ax=b 的解的充分必要条件是 是二次函
2、数 () = 12 的极小点，即 = () = () 证明：充分性 .设 是 f(x)的极小点，则由无约束最优化问题最优解的必要条件知 () = 即 是方程组 Ax=b 的解。 必要性 . 若 是方程组 Ax=b 的解，即 A=b，注意到 A 是对称正定矩阵，故x 有 ()() = 12 12 + = 12( 2 +) + = 12( 2() +)( ) = 12( )( ) 0 即 是 f(x)的极小点，进而由 A 是正定矩阵知， 是 f(x)的最小点。证毕。 共轭梯度法在形式上具有迭代法的特征，在给定初始值情况下，根据迭代公式： (+1) = () +() 产生的迭代序列 (1),(2),
3、(3) 在无舍入误差假定下，最多经过 n 次迭代，就可求得 f(x) 的最小值，也就是方程 Ax=b 的解。 共轭梯度法中关键的两点是，确定迭代格式中的搜索向量 ()和最佳步长 （ 0） 。 实际上，搜索方向 d()是关于矩阵 A 的共轭向量，在迭代中逐步构造之。步长 的确定原则是给定迭代点 ()和搜索方向 d()后，要求选取非负实数 ，使得 f() +d()达到最小，即选择 ，满足 () +() = 0 () +() 2 下面首先导出最佳步长 k 的计算公式。 假 设 迭 代 点 () 和 搜 索 方 向 d() 已 经 给 定 ， 可 以 通 过 一 元 函 数 () = () +() 的
4、极小化 () =0() +() 来求得，根据多元复合函数的求导法则得 : () = () +()() 令 0 = () = () +()() = () +()() = () +()() = () +()() 其中 () = ()。由此解得最佳步长 = ()()()() (3.4.4) 下面确定搜索方向 d()。 给定初始向量 (0)后，由于负梯度方向是函数下降最快的方向，故第 1 次迭代取搜索方向 (0) = (0) = (0) = (0) 。令 (1) = (0) +0(0) 其中 0 = (0)(0)(0)(0)。第 2次迭代时，从 (1) 出发的搜索方向不再取 (1)，而是 选取 (1)
5、 = (1) +0(0)，使得 (1)与 (0)是关于矩阵 A 的共轭向量，即要求 (1)满足 (1),(0) = 0，由此可求得参数 0: 0 = (1)(0)(0)(0) 然后从 (1)出发，沿 (1)进行搜索得到 (2) = (1) +1(1) 其中 1由式 （ 3.4.4） 确定。 一般地， 设已经求出 (+1) = () +()，计算 (+1) = (+1)。令 (+1) = (+1) +()，选取 ，使得 (+1)和 ()是关于 A 的共轭向量，即要求(+1)满足 (+1),() = 0。由此可得 : = (+1)()()() 从而确定了搜索方向 (+1)。 3 综上所述，共轭梯度
6、法的计算公式为 (0) = (0) = (0), = ()()()(), (+1) = () +(), (+1) = (+1), = (+1)()()() , (+1) = (+1) +(). 利用定理 3.4.3，可以将 和 的表达式简化为 = ()()()() = (+1)(+1)()() 定理 3.4.3 设分别是共轭梯度法中产生的非零残向量和搜索方向，则 （ 1） (),(1) = 0 （ 2） (),() = (),() （ 3） (0), (1), ()(k n1)是正交向量组， (0), (1), ()(k n1)是关于 A的共轭向量组。 关于共轭向量法的误差有如下定理： 定理
7、3.4.5 设是用共轭梯度法求得的迭代序列，则有误差估计式 () 2( 1 +1)(0) 其中范数 = , = 2(). 若 () = 0，则 ()就是线性方程组的解，故在计算过程中将 (+1) 作为迭代终止的条件。 共轭梯度法的 计算过程如下： （ 1） 给定初始近似向量 (0) 以及精度 0； （ 2） 计算 (0) = (0)，取 (0) = (0)； （ 3） for k=0 to n-1 do （ i） = ()()()()； （ ii） (+1) = () +()； （ iii） (+1) = (+1)； （ iv） 若 (+1) 或 k+1=n，则输出近似解 (+1)，停止；否则
8、，转（ v）； 4 （ v） = (+1)22()22 ； （ vi） (+1) = (+1) +()； end do 计算经验表明，对于不是十分病态的问题，共轭梯度法收敛较快，迭代次数远小于矩阵的阶数 n。对于病态问题，只要进行足够多次的迭代（迭代次数大约为矩阵阶数 n 的 3-5 倍）后，一般也能得到满意的结果，因而共轭梯度法是求解高阶稀疏线性方程组的一个有效、常用的方法 。 二、 程序框图 5 三、程序使用说明及案例计算结果 程序使用说明 本共轭梯度法求解线性方程的程序需要用户给定对称正定矩阵 A 的阶数，误差以及初始向量，然后程序自动调用定义的函数 Gongetidu(A,b,x0,e
9、psa)实现求解线性方程组 Ax=b。其中 A， b 的阶数，初始向量 x0 和 epsa 都是可变的。 案例计算结果 可以发现，当 n=100， 200， 300 时，方程组 Ax=b 的解 x 为元素均为 1 的 n 阶向量。 6 四、 Matlab 程序（ M-文件） clear A b x0;clc; %程序运行前首先清除 A， b 和 X0 的数值，以免影响计算 n=input(请输入对称正定矩阵 A 的阶数 n=); epsa=input(请输入误差 =); A=zeros(0,0); %给矩阵 A 预先配置内存空间，减少耗时 for k=1:(n-1) %读取题目 3.2 的矩阵
10、 A A(k,k)=-2; A(k,k+1)=1; A(k+1,k)=1; end A(n,n)=-2; A; b(1,1)=-1; %读取题目 3.2 的矩阵 b b(n,1)=-1; b; x0(1,1)=input(假定初始向量各元素相同，均为： ); %给题目 3.2 的迭代过程赋初值 for kk=2:n x0(kk,1)=x0(1,1); end x=Gongetidu(A,b,x0,epsa); %调用共轭梯度法求线性方程组的函数 function x=Gongetidu(A,b,x0,epsa) n=size(A,1); x=x0; r=b-A*x; d=r; for k=0:
11、(n-1) alpha=(r*r)/(d*A*d); x=x+alpha*d; r2=b-A*x; if (norm(r2)=eps R(2k)=C(2(k+1),a,b,n)+1/(43-1)*(C(2(k+1),a,b,n)-C(2k,a,b,n); eps1=abs(R(2k)-R(2(k-1); k=k+1; end y=R(2(k-1); end 22 四阶龙格 -库塔法 七、 算法原理 利用泰勒展开可以导出龙格 -库塔法（简称 R-K 方法）。 m 级 的龙格 -库塔法的 一般形式为 +1 = +11 +22 + 1 = (,) 2 = ( +2, +211) 3 = ( +3,
12、+311 +322) = ( +, +11 +22 +,11) 其中 ,均为常数，有待定系数法确定。确定的原则是将局部截断误差 =(+1)+1在处 泰勒展开，适当选取 h 的系数，使得局部截断误差 的阶数尽可能高。 下面以 m=2 级的 R-K 法为例来说明 R-K 法的导出思想。 2 级 R-K 法形式如下： +1 = +11 +22 1 = (,) 2 = ( +, +1) 将 (+1)在 处泰勒展开有 (+1) = ( +) = ()+() + 12!()2 + 13!()3 +(4) 由 () = (,()知 () = + = + () = +2 +2 + + 2 所以 (+1) =
13、()+ +12( + )2 +16( +2 +2 + + )3 +(4) 其中 为 (,()的简写，符号 , ,的含义亦然。 又由二元函数的泰勒展开式有 2 = ( +, +1) = ( + + +1222 +2 +12222 +) 于是 +1 = +11 +22 23 = +1 +2( + + +1222 +2 +12222 +) = +(1 +2) +2( + )2 +2 (122 + +1222)3 +(4) 可得局部截断误差 = (+1)+1 = (11 2) +(122) +(122) 2 +(161222) +(132) +(161222)2 +16( + 2)3 +(4) 为使
14、的阶尽可能高，由 的任意性，应选取 ,1,2，使得 ,2的系数为零，即得方程组 11 2 = 0 122 = 0 122 = 0 由此解出 ,1,2，代入相应表达式即可。注意到 ,1,2无论怎样选取， 3的系数中 + 2 0，故 = (3)，所以通过解上面的方程组得到的方法均是二阶方法。其含有 3 个方程， 4 个未知量，解有无穷个。 (1)取 1 = 2 = 12, = = 1，则得 +1 = +121 +122 1 = (,) 2 = ( +, +1) 即改进欧拉法。 (2)取 1 = 0,2 = 1, = = 1/2，则可得 +1 = +2 1 = (,) 2 = ( +12, +121
15、) 通常也写为 +1 = +( +12, +12(,) 即变形欧拉法。 24 类似于上述二阶方法的推导，可得多种 4 级 4 阶 R-K 法。应用最广泛的是如下标准（经典） 4 级 4 阶 R-K 法： +1 = +16(1 +22 +23 +4) 1 = (,) 2 = ( +12, +121) 3 = ( +12, +122) 4 = ( +, +3) 另外两个常用的 4 级 4 阶 R-K 法的计算公式不再列出。 4 阶以上的 R-K 法计算函数值的工作量较大，对于大量的实际问题 4 阶 R-K 法已可满足精度要求。若有更高的精度要求，可使用 6 级 5 阶英格兰公式（此处不再给出）。
16、R-K 方法的优缺点如下： 优点：精度高，它是显式单步法，计算 +1只需要 一个值，每一步的步长可根据精度的要求单独考虑其大小。 缺点：需要计算多个函数值，当 (,)比较复杂时，计算量较大。 八、 程序框图 原题目是高阶微分方程初值问题，需将其转化为一阶微分方程组来求解。令1 = ,2 = ,3 = ,转化后的方程如下： 1 = 2 2 = 3 3 = 3 +2 1 +2 3 1(0) = 1,2(0) = 3,3(0) = 3 1,2,3,4计算如下： 1 = (,) = (111213) 2 = ( +12, +121) = (212223) 3 = ( +12, +122) = (313
17、233) 4 = ( +, +3) = (414243) 25 程序框图如下： 开始 输入区间 a,b，步长 h 初值 1(),2(),3() 1 = (,) = (111213) 2 = ( +12, +121) = (212223) 3 = ( +12, +122) = (313233) 4 = ( +, +3) = (414243) +1 = +16(1 +22 +23 +4) 输出 () 结束 n=(b-a)/h， i=0 in? i=i+1 26 三、程序使用说明及案例计算结果 该程序只针对案例 9.2 的微分方程，对其他微分方程需修改程序。该程序需要用户输入区间 a,b的左右端点，
18、步长 ,1,2,3的在 a 处的初值，然后可得案例 9.2的 R-K 解和解析解，并计算了绝对误差和相对误差。应该注意输入的 a,b 之差应为为步长 h 的整数倍。案例计算结果如下： 四、 Matlab 程序（ M-文件） clear;clc; a=input(请输入区间左端点 :a=); b=input(请输入区间左端点 :b=); h=input(请输入步长 :h=); y1=input(请输入 y1 的初值 :y1(a)=); y2=input(请输入 y2 的初值 :y2(a)=); y3=input(请输入 y3 的初值 :y3(a)=); y0=f(b); y=RK4(a,b,y1
19、,y2,y3,h); d=y0-y; m=d/y0; fprintf(R-K 法解为 y(b)=%.8fn,y); fprintf(解析解为 y(b)=%.8fn,b,y0); fprintf(绝对误差为 %.8fn,d); fprintf(相对误差为 %.8fn,m); functiony0=f(b) y0=b*exp(b)+2*b-1; 27 end functiony=RK4(a,b,y1,y2,y3,h) x=a; n=(b-a)/h; for i=0:n-1 k11=h*y2; k12=h*y3; k13=h*(y3+y2-y1+2*(x+i*h)-3); k21=h*(y2+k12
20、/2); k22=h*(y3+k13/2); k23=h*(y3+k13/2+y2+k12/2-(y1+k11/2)+2*(x+i*h)+h/2)-3); k31=h*(y2+k22/2); k32=h*(y3+k23/2); k33=h*(y3+k23/2+y2+k22/2-(y1+k21/2)+2*(x+i*h)+h/2)-3); k41=h*(y2+k32); k42=h*(y3+k33); k43=h*(y3+k33+y2+k32-(y1+k31)+2*(x+i*h)+h)-3); y1=y1+(k11+2*k21+2*k31+k41)/6; y2=y2+(k12+2*k22+2*k32+k42)/6; y3=y3+(k13+2*k23+2*k33+k43)/6; end y=y1; end
