西安交通大学计算方法A上机作业

计算方法（A）大作业姓名：班级：专业：学号：共轭梯度法一、算法原理共轭梯度法是把求解线性方程组的问题转化为求解一个与之等价的二次函数极小化的问题，因此从任意给定的初始点出发，沿一组关于矩阵A的共轭方向进行线性搜索，在无舍入无差的假定下，最多迭代n(其中n为矩阵A的阶数）次就可求得二次函数的极小点，也就求得了线性方程组Ax=B的解。

下述定理给出了求系数矩阵A是对称正定矩阵的线性方程组Ax=b的解与求二次函数f(x)=12x T Ax−b T x极小点的等价性。

定理3.4.1设A是n阶对称正定矩阵，则x∗是方程组Ax=b的解的充分必要条件是x∗是二次函数f(x)=12x T Ax−b T x的极小点，即Ax∗=b⟺f(x∗)=minx∈R nf(x)证明：充分性.设x∗是f(x)的极小点，则由无约束最优化问题最优解的必要条件知∇f(x∗)=Ax∗−b即x∗是方程组Ax=b的解。

必要性. 若x∗是方程组Ax=b的解，即A x∗=b，注意到A是对称正定矩阵，故∀x∈R n有f(x)−f(x∗)=12x T Ax−b T x−12x T Ax∗+b T x∗=12(x T Ax−2b T x+x∗T Ax∗)−x∗T Ax∗+b T x∗=12(x T Ax−2(Ax∗)T x+x∗T Ax∗)−(Ax∗−b)T x∗=12(x−x∗)T A(x−x∗)≥0即x∗是f(x)的极小点，进而由A是正定矩阵知，x∗是f(x)的最小点。

证毕。

共轭梯度法在形式上具有迭代法的特征，在给定初始值情况下，根据迭代公式：x(k+1)=x(k)+αk d(k)产生的迭代序列x(1),x(2),x(3)…在无舍入误差假定下，最多经过n次迭代，就可求得f(x) 的最小值，也就是方程Ax=b的解。

共轭梯度法中关键的两点是，确定迭代格式中的搜索向量d(k)和最佳步长αk （αk≥0）。

实际上，搜索方向d(k)是关于矩阵A的共轭向量，在迭代中逐步构造之。

步长αk的确定原则是给定迭代点x(k)和搜索方向d(k)后，要求选取非负实数αk，使得f(x(k)+αk d(k))达到最小，即选择αk，满足f(x(k)+αk d(k))=min0≤αf(x(k)+αd(k))下面首先导出最佳步长k的计算公式。

假设迭代点x(k)和搜索方向d(k)已经给定，αk可以通过一元函数∅(α)=f(x(k)+αd(k))的极小化min∅(α)=0≤αf(x(k)+αd(k))来求得，根据多元复合函数的求导法则得:∅́(α)=∇f(x(k)+αd(k))T d(k)令0=∅́(α)=∇f(x(k)+αd(k))T d(k)=[A(x(k)+αd(k))−b]T d(k)=(Ax(k)−b+αAd(k))T d(k)=(−r(k)+αAd(k))T d(k)其中r(k)=b−Ax(k)。

由此解得最佳步长αk=r(k)T d(k)d(k)T Ad(k)(3.4.4)下面确定搜索方向d(k)。

给定初始向量x(0)后，由于负梯度方向是函数下降最快的方向，故第1次迭代取搜索方向d(0)=r(0)=−∇f(x(0))=b−Ax(0)。

令x(1)=x(0)+α0d(0)其中α0=r(0)T d(0)d(0)T Ad(0)。

第2次迭代时，从x(1)出发的搜索方向不再取r(1)，而是选取d(1)=r(1)+β0d(0)，使得d(1)与d(0)是关于矩阵A的共轭向量，即要求d(1)满足(d(1),Ad(0))=0，由此可求得参数β0:β0=−r(1)T Ad(0) d(0)T Ad(0)然后从x(1)出发，沿d(1)进行搜索得到x(2)=x(1)+α1d(1)其中α1由式（3.4.4）确定。

一般地，设已经求出x(k+1)=x(k)+αk d(k)，计算r(k+1)=b−Ax(k+1)。

令d(k+1)=r(k+1)+βk d(k)，选取βk，使得d(k+1)和d(k)是关于A的共轭向量，即要求d(k+1)满足(d(k+1),Ad(k))=0。

由此可得:βk=−r(k+1)T Ad(k) d Ad从而确定了搜索方向d(k+1)。

综上所述，共轭梯度法的计算公式为d(0)=r(0)=b−Ax(0),αk=r(k)T d(k)d(k)T Ad(k),x(k+1)=x(k)+αk d(k),r(k+1)=b−Ax(k+1),βk=−r(k+1)T Ad(k)d(k)T Ad(k),d(k+1)=r(k+1)+βk d(k).利用定理3.4.3，可以将αk和βk的表达式简化为αk=r(k)T r(k) d Adβk=r(k+1)T r(k+1) r(k)T r(k)定理3.4.3设分别是共轭梯度法中产生的非零残向量和搜索方向，则（1）(r(k),d(k−1))=0（2）(r(k),d(k))=(r(k),r(k))（3）r(0), r(1),…, r(k)(k≤n−1)是正交向量组，d(0), d(1),…, d(k)(k≤n−1)是关于A的共轭向量组。

关于共轭向量法的误差有如下定理：定理3.4.5设是用共轭梯度法求得的迭代序列，则有误差估计式‖x(k)−x∗‖A ≤2(√K)k‖x(0)−x∗‖A其中范数‖x‖A=T,K=Cond2(A).若r(k)=0，则x(k)就是线性方程组的解，故在计算过程中将‖r(k+1)‖≤ε作为迭代终止的条件。

共轭梯度法的计算过程如下：（1）给定初始近似向量x(0)以及精度ε>0；（2）计算r(0)=b−Ax(0)，取d(0)=r(0)；（3）for k=0 to n-1 do（i）αk =r(k)T d(k)d(k)T Ad(k)；（ii）x(k+1)=x(k)+αk d(k)；（iii）r(k+1)=b−Ax(k+1)；（iv）若‖r(k+1)‖≤ε或k+1=n，则输出近似解x(k+1)，停止；否则，转（v）；（v ）βk =‖r (k+1)‖22‖r (k)‖22；（vi ）d (k+1)=r (k+1)+βk d (k)； end do计算经验表明，对于不是十分病态的问题，共轭梯度法收敛较快，迭代次数远小于矩阵的阶数n 。

对于病态问题，只要进行足够多次的迭代（迭代次数大约为矩阵阶数n 的3-5倍）后，一般也能得到满意的结果，因而共轭梯度法是求解高阶稀疏线性方程组的一个有效、常用的方法。

二、 程序框图三、程序使用说明及案例计算结果程序使用说明本共轭梯度法求解线性方程的程序需要用户给定对称正定矩阵A的阶数，误差以及初始向量，然后程序自动调用定义的函数Gongetidu(A,b,x0,epsa)实现求解线性方程组Ax=b。

其中A，b的阶数，初始向量x0和epsa都是可变的。

案例计算结果可以发现，当n=100，200，300时，方程组Ax=b的解x为元素均为1的n阶向量。

四、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的矩阵AA(k,k)=-2;A(k,k+1)=1;A(k+1,k)=1;endA(n,n)=-2;A;b(1,1)=-1; %读取题目3.2的矩阵bb(n,1)=-1;b;x0(1,1)=input('假定初始向量各元素相同，均为：'); %给题目3.2的迭代过程赋初值for kk=2:nx0(kk,1)=x0(1,1);endx=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:(n-1)alpha=(r'*r)/(d'*A*d);x=x+alpha*d;r2=b-A*x;if ((norm(r2)<=epsa)||(k==n-1))disp('方程组的近似解为x=');disp(x);break;endbeta=norm(r2)^2/norm(r)^2;d=r2+beta*d;r=r2;endend三次样条插值三、算法原理分段低次插值多项式的一阶或二阶导数不连续，不能满足许多实际问题的要求。

例如，在飞机机翼外形设计，船体放样等许多实际问题中，要求曲线二阶光滑，即要求函数的二阶导数连续。

为了将一些已知点连成一条充分光滑的曲线，早期的做法是绘图人员用压铁把一条称为样条的富有弹性的细木条或金属条固定在相邻的若干点上，然后沿样条画下一条所需的光滑曲线，这条曲线称为样条曲线。

实际上，它是由分段三次曲线拼接而成的，在连接点处二阶导数连续。

将样条曲线抽象为数学问题称为样条函数。

1、三次样条插值函数的定义设在区间[a,b]上给定n+1个节点x i(a≤x0<x1<⋯<x n≤b)，在节点x i处的函数值y i=f(x i)(i=0,1,⋯,n)。

若函数S(x)满足以下三条：（1）在每个子区间[x i−1,x i]上，S(x)是三次多项式；（2）是个S(x i)=y i(i=0,1,⋯,n)；（3）在区间[a,b]上，S(x)的二阶导数S(x)连续，则称S(x)为函数y=f(x)在区间[a,b]上的三次样条插值函数。

由定义可知，S(x)是分段三次多项式，故在每个子区间[x i−1,x i] (i=1,⋯,n)上，S(x)有4个待定参数，共有n个子区间，所以S(x)共有4n个待定参数。

根据定义中的条件（3）知S−(x i)=S+(x i)，S−(x i)=S+(x i)，i=1,2,⋯,n−1.S−(x i)=S+(x i)，上式共有3n-3个条件，再加上定义条件（2）中的n+1个条件，共有4n-2个条件。

还需要增加两个条件才能确定4n个待定参数，即才能确定S(x)。

所增加的条件称为边界条件或端点条件。

三种常用的边界条件如下：（1）已知f(x)在区间[a,b]的两端点a，b处的二阶导数值f(a),f(b)；（2）已知f(x)在区间[a,b]的两端点a，b处的一阶导数值f(a),f(b)；（3）已知函数f(x)式以T=b-a为周期的周期函数。

2、三次样条插值函数的导出1)导出在子区间[x i−1,x i](i=1,2,⋯,n)上的S(x)的表达式由于S(x)的二阶导数连续，设S(x)在节点x i处的二阶导数值S(x i)=M i(i= 1,2,⋯,n)，其中M i为未知的待定参数。