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

计算方法（A）大作业姓名：班级：专业：学号：1共轭梯度法一、算法原理共轭梯度法是把求解线性方程组的问题转化为求解一个与之等价的二次函数极小化的问题，因此从任意给定的初始点出发，沿一组关于矩阵A的共轭方向进行线性搜索，在无舍入无差的假定下，最多迭代n(其中n为矩阵A的阶数）次就可求得二次函数的极小点，也就求得了线性方程组Ax=B的解。下述定理给出了求系数矩阵A是对称正定矩阵的线性方程组Ax=b的解与求二次函数𝑓(𝑥)=12𝑥𝑇𝐴𝑥−𝑏𝑇𝑥极小点的等价性。定理3.4.1设A是n阶对称正定矩阵，则𝑥∗是方程组Ax=b的解的充分必要条件是𝑥∗是二次函数𝑓(𝑥)=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)…在无舍入误差假定下，最多经过n次迭代，就可求得f(x)的最小值，也就是方程Ax=b的解。共轭梯度法中关键的两点是，确定迭代格式中的搜索向量𝒅(𝒌)和最佳步长𝜶𝒌（𝜶𝒌≥0）。实际上，搜索方向d(𝑘)是关于矩阵A的共轭向量，在迭代中逐步构造之。步长𝛼𝑘的确定原则是给定迭代点𝑥(𝑘)和搜索方向d(𝑘)后，要求选取非负实数𝛼𝑘，使得f(𝑥(𝑘)+𝛼𝑘d(𝑘))达到最小，即选择𝛼𝑘，满足𝑓(𝑥(𝑘)+𝛼𝑘𝑑(𝑘))=𝑚𝑖𝑛0≤𝛼𝑓(𝑥(𝑘)+𝛼𝑑(𝑘))2下面首先导出最佳步长k的计算公式。假设迭代点𝑥(𝑘)和搜索方向d(𝑘)已经给定，𝛼𝑘可以通过一元函数∅(𝛼)=𝑓(𝑥(𝑘)+𝛼𝑑(𝑘))的极小化𝑚𝑖𝑛∅(𝛼)=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)=𝑟(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综上所述，共轭梯度法的计算公式为𝑑(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≤n−1)是正交向量组，𝑑(0),𝑑(1),…,𝑑(𝑘)(k≤n−1)是关于A的共轭向量组。关于共轭向量法的误差有如下定理：定理3.4.5设是用共轭梯度法求得的迭代序列，则有误差估计式‖𝑥(𝑘)−𝑥∗‖𝐴≤2(√𝐾−1√𝐾+1)𝑘‖𝑥(0)−𝑥∗‖𝐴其中范数‖𝑥‖𝐴=√𝑥𝑇𝐴𝑥,𝐾=𝐶𝑜𝑛𝑑2(𝐴).若𝑟(𝑘)=0，则𝑥(𝑘)就是线性方程组的解，故在计算过程中将‖𝑟(𝑘+1)‖≤𝜀作为迭代终止的条件。共轭梯度法的计算过程如下：（1）给定初始近似向量𝑥(0)以及精度ε0；（2）计算𝑟(0)=𝑏−𝐴𝑥(0)，取𝑑(0)=𝑟(0)；（3）fork=0ton-1do（i）𝛼𝑘=𝑟(𝑘)𝑇𝑑(𝑘)𝑑(𝑘)𝑇𝐴𝑑(𝑘)；（ii）𝑥(𝑘+1)=𝑥(𝑘)+𝛼𝑘𝑑(𝑘)；（iii）𝑟(𝑘+1)=𝑏−𝐴𝑥(𝑘+1)；（iv）若‖𝑟(𝑘+1)‖≤𝜀或k+1=n，则输出近似解𝑥(𝑘+1)，停止；否则，转（v）；4（v）𝛽𝑘=‖𝑟(𝑘+1)‖22‖𝑟(𝑘)‖22；（vi）𝑑(𝑘+1)=𝑟(𝑘+1)+𝛽𝑘𝑑(𝑘)；enddo计算经验表明，对于不是十分病态的问题，共轭梯度法收敛较快，迭代次数远小于矩阵的阶数n。对于病态问题，只要进行足够多次的迭代（迭代次数大约为矩阵阶数n的3-5倍）后，一般也能得到满意的结果，因而共轭梯度法是求解高阶稀疏线性方程组的一个有效、常用的方法。二、程序框图5三、程序使用说明及案例计算结果程序使用说明本共轭梯度法求解线性方程的程序需要用户给定对称正定矩阵A的阶数，误差以及初始向量，然后程序自动调用定义的函数Gongetidu(A,b,x0,epsa)实现求解线性方程组Ax=b。其中A，b的阶数，初始向量x0和epsa都是可变的。案例计算结果可以发现，当n=100，200，300时，方程组Ax=b的解x为元素均为1的n阶向量。6四、Matlab程序（M-文件）clearAbx0;clc;%程序运行前首先清除A，b和X0的数值，以免影响计算n=input('请输入对称正定矩阵A的阶数n=');epsa=input('请输入误差=');A=zeros(0,0);%给矩阵A预先配置内存空间，减少耗时fork=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的迭代过程赋初值forkk=2:nx0(kk,1)=x0(1,1);endx=Gongetidu(A,b,x0,epsa);%调用共轭梯度法求线性方程组的函数functionx=Gongetidu(A,b,x0,epsa)n=size(A,1);x=x0;r=b-A*x;d=r;fork=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;endend7三次样条插值三、算法原理分段低次插值多项式的一阶或二阶导数不连续，不能满足许多实际问题的要求。例如，在飞机机翼外形设计，船体放样等许多实际问题中，要求曲线二阶光滑，即要求函数的二阶导数连续。为了将一些已知点连成一条充分光滑的曲线，早期的做法是绘图人员用压铁把一条称为样条的富有弹性的细木条或金属条固定在相邻的若干点上，然后沿样条画下一条所需的光滑曲线，这条曲线称为样条曲线。实际上，它是由分段三次曲线拼接而成的，在连接点处二阶导数连续。将样条曲线抽象为数学问题称为样条函数。1、三次样条插值函数的定义设在区间[a,b]上给定n+1个节点𝑥𝑖(𝑎≤𝑥0𝑥1⋯𝑥𝑛≤𝑏)，在节点𝑥𝑖处的函数值𝑦𝑖=𝑓(𝑥𝑖)(𝑖=0,1,⋯,𝑛)。若函数S(x)满足以下三条：（1）在每个子区间[𝑥𝑖−1,𝑥𝑖]上，S(x)是三次多项式；（2）是个S(𝑥𝑖)=𝑦𝑖(𝑖=0,1,⋯,𝑛)；（3）在区间[a,b]上，S(x)的二阶导数𝑆(𝑥)̈连续，则称S(x)为函数y=f(x)在区间[a,b]上的三次样条插值函数。由定义可知，S(x)是分段三次多项式，故在每个子区间[𝑥𝑖−1,𝑥𝑖](𝑖=1,⋯,𝑛)上，S(x)有4个待定参数，共有n个子区间，所以S(x)共有4n个待定参数。根据定义中的条件（3）知𝑆−(𝑥𝑖)=𝑆+(𝑥𝑖)，𝑆−̇(𝑥𝑖)=𝑆+̇(𝑥𝑖)，𝑖=1,2,⋯,𝑛−1.𝑆−̈(𝑥𝑖)=𝑆+̈(𝑥𝑖)，上式共有3n-3个条件，再加上定义条件（2）中的n+1个条件，共有4n-2个条件。还需要增加两个条件才能确定4n个待定参数，即才能确定S(x)。所增加的条件称为边界条件或端点条件。三种常用的边界条件如下：（1）已知f(x)在区间[a,b]的两端点a，b处的二阶导数值𝑓(𝑎)̈,𝑓(𝑏)̈；（2）已知f(x)在区间[a,b]的两端点a，b处的一阶导数值𝑓(𝑎)̇,𝑓(𝑏)̇；（3）已知函数f(x)式以T=b-a为周期的周期函数。2、三次样条插值函数的导出1)导出在子区间[𝒙𝒊−𝟏,𝒙𝒊](𝒊=𝟏,𝟐,⋯,𝒏)上的S(x)的表达式由于S(x)的二阶导数连续，设S(x)在节点𝑥𝑖处的二阶导数值𝑆(𝑥𝑖)̈=𝑀𝑖(𝑖=1,2,⋯,𝑛)，其中𝑀𝑖为未知的待定参数。由S(x)是分段的三次多项式知，𝑆(𝑥)̈是分段线性函数，𝑆(𝑥)̈在子区间[𝑥𝑖−1,𝑥𝑖]上可表示为8𝑆(𝑥)̈=𝑥−𝑥𝑖𝑥𝑖−1−𝑥𝑖𝑀𝑖−1+𝑥−𝑥𝑖−1𝑥𝑖−𝑥𝑖−1𝑀