四阶微分方程Hermite有限元的预处理共轭梯度法
第26卷第2期大 学 数 学V ol.26, .2 2010年4月COLLEGE M AT H EM AT ICS Apr.2010四阶微分方程Hermite有限元的预处理共轭梯度法
赵国忠
(包头师范学院数学科学学院,包头014030)
[摘 要]对一类四阶微分方程两点边值问题的H ermit e有限元方法进行了研究.首先讨论了该方程通常意义下的Galerkin有限元离散,考虑到有限元离散得到的线性方程组的对称正定性,文中采用了预处理最速下降法和共轭梯度方法求解线性方程组,通过选择不同的预处理器,使得求解该方程组的迭代次数有了很大的改观.
[关键词]H ermit e有限元;四阶微分方程;最速下降法;共轭梯度法;预处理器
[中图分类号]O224 [文献标识码]A [文章编号]1672 1454(2010)02 0047 05
1 模型方程及Hermite有限元离散
本文考虑如下的四阶微分方程
( (x)u(x))=f(x), 0 (1) u(0)=u(L)=0, u!(0)=u!(L)=0, 其中 (x)在(0,L)上有界,且 (x)>0( x?(0,L)),f(x)?L2(0,L). 首先建立该方程的Galerkin有限元格式[1,2],为此用具有一定光滑性的函数v乘方程的两端,并在(0,L)上积分,可得 #L0 u v d x- u v L0+ u v!L0=#L0f v d x. 由此可得问题(1)的弱形式为 求u?V,使得 v?V,都有 #L0 u v d x=#L0f v d x,(2)其中 V=H20(0,L)={v?v(r)?L2(0,L),r=0,1,2,v(r)(0)=v(r)(L)=0,r=0,1}. 为得到离散的有限元格式,用剖分 h将区间[0,L]等距分成K个子区间T k=[x k-1,x k],(k=1,%,K), x k=kh,h=L/K.定义有限元空间 V h={v h?C1(0,L),v h?P3(T), T? h,v(r)h(0)=v(r)h(L)=0,r=0,1},(3)从而原问题的Galer kin有限元格式可定义为 求u h?V h,使得 #L0 u h v h d x=#L0f v h d x, v h?V h.(4)以下定义有限元空间的H ermite基函数:对于任给的内部节点x i(i=1,%,K-1),定义其最小支集为i=T i&T i+1,每个节点对应的基函数!i和?i满足 k,!i|T k?P3(T k),?i|T k?P3(T k),并且 [收稿日期]2007 07 02 j =0,%,K ,都有 !i (x j )=#ij ,!i !(x j )=0,?i (x j )=0, ?i !(x j )=#ij . (5) 定义H ermite 基函数的集合 B h ={!i ,?i ,i =1,%,K -1}.(6) 以下采用单元分析法.设标准单元T ^=[0,1]. x ^?[0,1],借助于映射x =hx ^+x k -1,可将标准单元T ^ 变换到单元T k (k =1,%,K -1),由条件(5)可以得到与节点x ^=0对应的基函数^!(0)0,^!(1) 0及与节点x ^=1对应的基函数^!(0)1,^!(1)1为 ^!(0) 0(x ^)=1-3x ^2 +2x ^3 ,^!(1) 0=x ^-2x ^2 +x ^3 ,^!(0)1(x ^)=3x ^2-2x ^3, ^!(1) 1=-x ^2+x ^3.(7) 函数u h ?V h 可以记为 u h (x )= ? K-1 i=1 u i !i (x )+ ? K-1 i=1 u (1)i ?i (x ), (8) 其中u i =u h (x i ),u (1) i =u h !(x i ),i =1,%,K -1.将(8)代入(4)并令v h 取遍(6)式中的所有函数,可得如下的方程组 ?K-1 j=1u j #L 0 !j !i d x +u (1) j #L 0 ?j !i d x =#L f !i d x , ? K-1j=1 u j # L !j ?i d x +u (1) j #L 0 ?j ?i d x =#L f ?i d x , (i =1,%,K -1).(9) 为计算简单,以下假定L =1, 和f 为常数,计算(9)式中的积分并整理可得方程组 Au +Bp =b 1, B T u +Cp =0, (10) 其中u,p ? k -1 ,u =[u 1,%,u K -1]T ,p =[u (1)1,%,u (1)K -1]T ,b 1? K -1 ,b 1的元素为h 4f / . A =diag K -1(-12,24,-12), B =diag K -1(-6,0,6), C =diag K -1(2,8,2)为三对角矩阵,在(10)中消去p ,可得K -1阶线性方程组 (A -BC -1 B T )u =b 1. (11) 令M =A -BC -1 B T ,则M 是对称正定的,方程(11)可以用共轭梯度法求解. 2 预处理共轭梯度法的应用 最速下降法和共轭梯度法(简称CG 法)是求解大型线性代数方程组的有效方法.共轭梯度法是20世纪50年代初期由H estenes 和Stiefel 首先提出的.由于它是理论上有限步终止的,故仅仅被视为 是直接法,而计算上又没有表现出任何优势,所以在其后的20年没有受到重视[3].20世纪70年代初共轭梯度法被首次作为迭代法研究,共轭梯度法的收敛速度与?-1?+1成正比,其中的?=cond (A) =(A ()(A -1(= %max (A ) %min (A ) .由此可知只要求解的线性方程组的系数矩阵仅有少数几个互不相同的 特征值或者是十分良态(即??1)时,共轭梯度法收敛的就会很快.因此在应用共轭梯度法时,首先应该 将方程组转化为系数矩阵仅有少数几个互不相同的特征值或者是十分良态的等价方程组,然后再应用共轭梯度法于转化后的方程组.预处理共轭梯度法[4,5](简称PCG 法)正是基于这一基本思想[6]而产生的.他是先将方程组 Ax =b (12)转化为 A ~ x = b , (13) 其中A ~ =C -1AC -1, x =Cx, b =C -1b,这里要求C 是对称正定的,其目的是通过C 的选择,使A ~ 具有我们 48 大 学 数 学 第26卷 所希望的良好性质.然后应用共轭梯度法于方程(13). 矩阵C 称作预优矩阵,常用的预优矩阵的选取技巧有对角预优矩阵、不完全Cholesky 因子预优矩阵、多项式预优矩阵等.为计算简单起见,以下假定 (x )=f =1.采用最速下降法和共轭梯度法求解,求解过程中首先采用通常的最速下降法和共轭梯度法,然后通过选取对角预优矩阵和一类比较特殊的预优矩阵用预处理的最速下降法和共轭梯度法重新求解线性代数方程组,将相同精度要求下的求解迭代次数进行了比较,并将我们得到的数值解与原方程的精确解进行了比较,具体的计算结果及相关程序如下: functio n zgzpcg 4(K ,K 1) %函数中的K 和K 1分别表示有限元单元剖分的数量和数值求解线性代数方程组的最大迭代次数.h=1/K ;%步长 A=diag (24.*ones(1,K-1),0)+diag(-12.*o nes(1,K -2),1)+diag (-12.*ones(1,K-2),-1);B=diag(0.*ones(1,K -1),0)+diag(6.*ones(1,K -2),1)+diag(-6.*ones(1,K -2),-1);C1=diag(8.*ones(1,K-1),0)+diag(2.*ones(1,K-2),1)+diag(2.*ones(1,K-2),-1);A2=A -B*inv(C1)*B !; %以上部分是为了生成方程组的系数矩阵V=12.*ones(1,K -1);V(1)=10;V(K -1)=10;C=diag (V ,0) A1=A -B*inv(C)*B !;%预优矩阵 b=h ^4.*ones(K -1,1);方程组的右端项x=zero s(K -1,1);%迭代的初值 tol=1.0000e-010;%求解线性代数方程组Ax =b 的精度要求M =A1;t=0:h:1; z1=1/24.*t.^4-1/12.*t.^3+1/24.*t.^2;%微分方程的精确解maxit=K1; [y ,err or ,niter ,flag]=co njg rad(A 2,x,b,M ,max it ,to l) %调用求解方程组A x=b 的函数conjgr ad,其中y 表示方程组的数值解,er ror 表示经过max it 次迭代后方程组数值解的误差,niter 表示达到精度要求to l 所需的迭代次数.z2=[0,y !,0]; plot(t,z1,t,z2,!r--!)%通过图像比较数值解与精确解的差别. 表1 最速下降法迭代次数的比较 K Witho ut P recond diag M 254267024246203750--38100--36200 --34 表2 共轭梯度法迭代次数的比较 49 第2期 赵国忠:四阶微分方程H ermite 有限元的预处理共轭梯度法 在表1和表2中,K 表示有限元剖分的单元数,Without Preco nd 表示没有经过预处理时的情形,Diag 表示取对角预优矩阵的情形,M 表示预优矩阵取为A -B !C -1B T 的情形,其中!C 是一个对角矩阵,其主对角线上的元素满足 ?c ii = ? K-1 j=1 |c ij |, 精度要求为To l=10-10.由表1可以看到,没有预处理和采用对角预处理矩阵时的最速下降法收敛速度非常慢,达到所需精度要求所需的迭代次数较多,随着剖分单元的增加,在普通的微机上已经算不到所需的精度要求.而采用预优矩阵M 时,最速下降法达到所需的精度要求所需的迭代次数锐减.但从表2可以看到,达到相同的精度要求,(预处理)共轭梯度法所需的迭代次数均较少,尤其是剖分单元较多时采用预处理矩阵M ,共轭梯度法的迭代次数有了数量级上的变化.这充分显示了预处理技巧与共轭梯度法结合所产生的强大威力.同时也要看到,采用对角预优矩阵对迭代次数的改善并不明显,因此预优矩阵的选择对预处理共轭梯度法的计算效果有很大的影响.为了进一步比较共轭梯度法和最速下降法的计算效果,我们限定迭代次数取为2000时,得到的计算误差如下,其中n 次迭代后的误差计算公式为 er ror= (r niter (2 (b (2 .(14) 表3 最速下降法迭代次数为2000时的误差 K Without Pr eco nd diag M 25 1.0402 1.00007.0131e-01150 1.2223 1.1823 6.2775e-011100 1.3034 1.2632 3.0623e-011200 1.3416 1.3010 3.8675e-011 表4 共轭梯度法迭代次数为2000时的误差 K Without Pr eco nd diag M 257.1009e-0117.5040e-011 4.3226e-01150 3.3793e-0117.0990e-011 5.4211e-011100 6.9309e-0117.7616e-011 3.6950e-011200 1.7896e-008 4.8960e-011 6.8386e-017 本文进一步讨论了微分方程有限元解与精确解之间的误差,为直观起见,计算结果如下图. 图1 共轭梯度法下微分方程精确解与 图2 最速下降法下微分方程精确解与 数值解的比较(K =25,M axit=100) 数值解的比较(K =25,M ax it=100) 50大 学 数 学 第26卷 [参 考 文 献] [1] 王奇生,邓康.四阶方程两点边值问题H ermite 有限元解的渐进展式与外推[J].高等学校计算数学学报,2006, 28(4):299-306. [2] St ig L arsson,Vidar T ho mee.P artial differential equatio ns wit h numer ical met ho d[M ].北京:科学出版社,2006: 15-73. [3] 白峰杉.数值计算引论[M ].北京:高等教育出版社,2004:180-198. [4] Reid J K.T he use o f conjugat e g radient for sy stems of equat ions possessing pro per ty A [J].SIA M J.N umerical A na lysis,1972,9:325-332. [5] Ker shaw D S.T he incomplete cho leski conjugat e g radient metho d fo r the iterat ive solut ion of systems of linear equations[J].J Co mputatio nal Phy sics,1978,26:43-65. [6] 徐树方,高立,张平文.数值线性代数[M ].北京:北京大学出版社,2000:139-160. Preconditioned Conjugate Gradient Method of Hermite Finite Element Approximations for Two Point Boundary Value Problem of Fourth Order Differential Equations ZH A O Guo z hong (Facult y of M athematics,Baoto u T eacher +s College,Bao tou 014030,China) Abstract:H ermite finite element appro x imatio n fo r tw o po int bo undary value pro blem of fourth or der differential equat ions is discussed.Co nsidering the sy mmet ric and positive of the discrete matr ix ,steepest decent method and pr eco ndit ioned conjug ate gr adient method ar e used to solve the equatio n.Differ ent preconditio ned mat rix es a re used to reduce the iterat ive number. Key words:H ermite finite element;fo ur th o rder differentia l equatio n;steepest decent method;preconditio ned co njug ate g radient met ho d;pr econdit ioner 51 第2期 赵国忠:四阶微分方程H ermite 有限元的预处理共轭梯度法 最优化课程设计--共轭梯度法算法分析与实现(设计程序) 题目共轭梯度法算法分析与实现 班级 / 学号 14140101/2011041401011 学生姓名黄中武指导教师王吉波王微微 课程设计任务书 课程名称最优化方法课程设计院(系) 理学院专业信息与计算科学 课程设计题目共轭梯度法算法分析与实现课程设计时间: 2014 年 6月 16日至 2014 年 6月 27日 课程设计的要求及内容: [要求] 1. 学习态度要认真，要积极参与课程设计，锻炼独立思考能力; 2. 严格遵守上机时间安排; 3. 按照MATLAB编程训练的任务要求来编写程序; 4. 根据任务书来完成课程设计论文; 5. 报告书写格式要求按照沈阳航空航天大学“课程设计报告撰写规范”; 6. 报告上交时间:课程设计结束时上交报告; 7. 严禁抄袭行为，一旦发现，课程设计成绩为不及格。 一、运用共轭梯度法求解无约束最优化问题 要求:1)了解求解无约束最优化问题的共轭梯度法; 2)绘出程序流程图; 3)编写求解无约束最优化问题的共轭梯度法MATLAB程序; 4)利用编写文件求解某无约束最优化问题; 5)给出程序注释。 指导教师年月日 负责教师年月日 学生签字年月日 沈阳航空航天大学 课程设计成绩评定单 课程名称最优化理论与算法课程设计院(系) 理学院专业信息与计算科学课程设计题目共轭梯度法算法分析与实现学号 2011041401011 姓名黄中武指导教师评语: 课程设计成绩 指导教师签字 年月日 最优化方法课程设计沈阳航空航天大学课程设计用纸目录 目录 一、正 文 (1) 二、总结 ............................................................... 8 参考文 献 ............................................................... 9 附录 .. (10) 第 I 页 最优化方法课程设计沈阳航空航天大学课程设计用纸正文 一、正文 一无约束最优化问题的共轭梯度法 在科学技术各领域中，有很多问题都可以归结为偏微分方程问题。在物理专业得力学、热学、电学、光学、近代物理课程中都可遇见偏微分方程。 偏微分方程,再加上边界条件、初始条件构成得数学模型，只有在很特殊情况下才可求得解析解。随着计算机技术得发展,采用数值计算方法,可以得到其数值解。 偏微分方程基本形式 而以上得偏微分方程都能利用PDＥ工具箱求解。 ＰDE工具箱 ＰＤE工具箱得使用步骤体现了有限元法求解问题得基本思路,包括如下基本步骤： 1) 建立几何模型 2）定义边界条件 3) 定义ＰＤE类型与ＰＤE系数 ４）三角形网格划分 5) 有限元求解 ６）解得图形表达 以上步骤充分体现在ＰDE工具箱得菜单栏与工具栏顺序上，如下 具体实现如下。 打开工具箱 输入pdetool可以打开偏微分方程求解工具箱，如下 首先需要选择应用模式,工具箱根据实际问题得不同提供了很多应用模式，用户可以基于适当得模式进行建模与分析。 在Ｏpｔｉons菜单得Ａｐplicａｔｉoｎ菜单项下可以做选择,如下 或者直接在工具栏上选择，如下 列表框中各应用模式得意义为： ①Gｅｎｅric Scalaｒ：一般标量模式(为默认选项）。 ② GeneｒiｃSｙstｅm：一般系统模式. ③ Stｒucｔural Ｍecｈ、,Plane Streｓs：结构力学平面应力。 ④ Ｓtructurａl Meｃh、,Plane Sｔrain：结构力学平面应变。 ⑤Elｅctroｓtａtｉcｓ:静电学。 ⑥ Mａgnｅtｏｓtaｔiｃs：电磁学。 ⑦Ａc Power Ｅｌecｔｒomagneｔics:交流电电磁学。 ⑧CｏnducｔｉvｅＭedｉa DＣ:直流导电介质。 ⑨ Heａｔ Traｎｆer：热传导。 ⑩ Ｄiｆｆｕsioｎ：扩散。 可以根据自己得具体问题做相应得选择，这里要求解偏微分方程，故使用默认值。此外，对于其她具体得工程应用模式，此工具箱已经发展到了ｓoｌMulｔiphｙsｉcｓ软件,它提供了更强大得建模、求解功能。 另外，可以在菜单Opｔioｎｓ下做一些全局得设置，如下 l Grｉd:显示网格 l Ｇrid Spacing…:控制网格得显示位置 l Sｎap：建模时捕捉网格节点，建模时可以打开 l Axeｓ Limｉts…：设置坐标系范围 ｌ Axes Equal:同Mａtlａｂ得命令ａxes eｑｕal命令 建立几何模型 使用菜单Ｄｒaw得命令或使用工具箱命令可以实现简单几何模型得建立，如下 各项代表得意义分别为 ｌ绘制矩形或方形； l 绘制同心矩形或方形； 16梯度法和共轭梯度法基本原理和特点？ 梯度法又称最速下降法，基本原理是在迭代点附近采用使目标函数值下降最快的负梯度方向作为搜索方向， 求目标函数的极小值，特 点；迭代计算简单，只需求一阶偏导数，所占的存储单 元少，对初始点的要求不 高，在接近极小点位置时收敛速度很慢，共轭的特点为 在梯度法靠近极值点收敛速度放慢时，它可以构造共轭方向使其收敛速度加快， 迭代计算比较简单，效果 好，在每一步迭代过程中都要构造共轭的、方向，比较繁琐。 17迭代终止准则有哪三种？ 1）当设计变量在相邻两点之间的移动距离充分小时，可用相邻两点的矢量差的模作为终止的判据， 2）当相邻两点目标函数值之差达到充分小时，可 用两次迭代的目标函数之 差作为终止判据。 3）当迭代点逼近极值点时，目标函数在该点的梯度已达到充分小时，可用梯度的模作为 终止判据。 18 .无约束设计法，1）powell法，它是在下降迭代过运算中只需计算和比较目标函数值的大小，不需计算偏导数的方法，是较好的一种直接搜索 算法。 2）梯度法，又称最速下降法，它是采用使目标 函数值下降最快的负梯度方向作为搜索方向来求目标函数的极小值。 3）共轭梯度法，又称FR法，是利用目标函数的梯度确定共轭方向，使得计算简便而效 果好，只需利用相邻两点的梯度就可以构造一个共轭方向，这种方式产生共轭方 向并进行迭代的算法称为 共轭梯度法。 4）变尺度法，又称DFP法，为了得到既有快速收敛的性质，又能避免计算二阶导数矩阵及逆矩阵，减少计算工作量。迭代公式X=X+aS, 19有约束设计法？ 1）复合形法，在可行域中选取k个设计点作为 初始复合形的顶点，然后比较复合形个各项 目标函数值的大小，其中目标函数值最大的点为坏点，以坏点之外其余各点的中心为映 射中心，寻坏点的映射点，以映射点替换坏点，并与原复合 型除坏点之外其余各点构成就k 顶点的新的复合型，这样反复迭代直到达到精度找到最优点， 2）简约梯度法，用来解决线性约束非线性规划问题。3）罚函数法，是把一个有约束的问 题转化为一系列无约束的问 题求解，逐渐逼近最优值。 . 可靠性工程包括的三个方面？ 1可靠性设计，包括设计方面的分析，对比评价，必要时也包 括可靠性实验，生产制造中的质量控制设计 及使用维修规程的设计。 2可靠性分析，主要是失效分析，也包括故障分析 3可靠性数学， 这是数理统计方法在开展 可靠性工作中发展起来的 数学分支。 常用的可靠 偏微分方程数值解 （Numerical Methods for Partial Differential Equations） 课程代码：10210801 学位课程/非学位课程：非学位课程 学时/学分：46/3 课程简介： 《偏微分方程数值解》是数学类专业必修的一门专业课。主要内容包括：变分形式和Galerkin有限元法、椭圆型方程的差分方法、抛物型方程的差分方法、双曲型方程的差分方法、离散方程的解法。通过本课程的学习，使学生掌握求解偏微分方程数值解的基本方法，能够根据具体的微分方程使用合适的计算方法。 一、教学目标 1、知识水平教学目标 偏微分方程数值解课程的教学，要使学生掌握椭圆型微分方程、抛物型微分方程、双曲型微分方程等典型方程的差分方法，了解与之相关的理论问题，理解变分原理、有限元方法以及离散方程的解法，理解各种计算方法的收敛条件和收敛速度。 2、能力培养目标 通过偏微分方程数值解课程教学,应注意培养学生以下能力： （1）连续问题离散化能力——掌握科学的思维方法，能够使用差分方法和有限元方法的各种格式对三类典型方程进行离散化处理。 （2）算法分析与设计能力——结合各类偏微分方程的特点，设计各种计算方法，对计算方法的收敛条件和收敛速度等进行分析，具体设计易于上机实现的算法。（3）离散方程组的快速求解能力——理解离散方程组的特点，使用数学软件编程，具体上机实现，进行数值模拟的动手能力。 3、素质培养目标 通过数学物理方程课程教学,应注重培养学生以下素质： （1）具体问题有限化——善于对现实世界中得到的偏微分方程进行有限差分、有限元分析的有限化思想素养。 （2）数值解法定性化——通过学习，引导学生树立偏微分方程数值求解的基本原则，培养学生对数值方法中的稳定性、收敛性和误差等进行定性分析的素质。（3）算法实现程序化——培养学生的创造性和具体实现程序化的思维，使学生学会用数学中算法的观点思考实际问题，用程序和计算机解决数学问题。 二、教学重点与难点 1、教学重点：椭圆型、抛物型、双曲型等微分方程的差分方法，有限元方法。 2、教学难点：各种计算方法的稳定性、收敛性和误差分析，变分形式。 三、教学方法与手段 以教师讲授为主，安排上机实验，辅以习题课、课堂讨论、小论文，注重理论联系实际。 四、教学内容与目标 教学内容教学目标课时分配 （46学时） 1. 边值问题的变分形式 6 二次函数的极值掌握 两点边值问题掌握 共轭梯度法 1.算法思想： 共轭梯度法是利用目标函数梯度逐步产生共轭方向作为线搜索方向的方法，每次搜索方向都是在目标函数梯度的共轭方向，搜索步长通过一维极值算法确定。 2.算法步骤： 用共轭梯度法求无约束多维极值问题min (),n f x x R ∈的算法步骤如下： （1） 给定初始点(0)x ，及精度0ε>； （2） 若 (0)()f x ε ?≤，停止，极小值点为(0)x ，否则转步骤（3）； （3） 取(0)(0)()p f x =-?，且置0k =； （4） 用一维搜索法求k t ，使得()()()()()0 ()min k k k k k t f x t p f x tp ≥+=+，令，(1)()()k k k k x x t p +=+，转步骤5； （5） 若 (1)()k f x ε +?≤，停止，极小值点为(1)k x +，否则转步骤（6）； （6） 若1k n +=，令(0)()n x x =，转步骤（3），否则转步骤（7）； （7） 令(1)(1)()()k k k k p f x p λ++=-?+，2(1)2 () () () k k k f x f x λ+?=?，置1k k =+，转 步骤（4）。 3.算法源程序： #include #define N 10 #define eps pow(10,-6) double f(double x[],double p[],double t) { double s; s=pow(x[0]+t*p[0],2)+25*pow(x[1]+t*p[1],2); return s; } /*以下是进退法搜索区间源程序*/ void sb(double *a,double *b,double x[],double p[]) { double t0,t1,t,h,alpha,f0,f1; int k=0; t0=2.5; /*初始值*/ h=1; /*初始步长*/ alpha=2; /*加步系数*/ f0=f(x,p,t0); t1=t0+h; f1=f(x,p,t1); while(1) { #include printf("\n"); printf("input the Maxtr X\n"); /*给X0输入一个初始向量*/ for(i=0;i 一、共轭梯度法 共轭梯度法（Conjugate Gradient）是共轭方向法的一种，因为在该方向法中每一个共轭向量都是依靠赖于迭代点处的负梯度而构造出来的，所以称为共轭梯度法。由于此法最先由Fletcher和Reeves （1964）提出了解非线性最优化问题的，因而又称为FR 共轭梯度法。由于共轭梯度法不需要矩阵存储，且有较快的收敛速度和二次终止性等优点，现在共轭梯度法已经广泛地应用于实际问题中。共轭梯度法是一个典型的共轭方向法，它的每一个搜索方向是互相共轭的，而这些搜索方向d仅仅是负梯度方向与上一次迭代的搜索方向的组合，因此，存储量少，计算方便，效果好。 二、共轭梯度法的原理 设有目标函数 f（X）=1/2X T HX+b T X+c 式1 式中，H作为f（X）的二阶导数矩阵，b为常数矢量，b=[b1,b2,b3,...b n]T 在第k次迭代计算中，从点X（k）出发，沿负梯度方向作一维搜索，得 S(K)=-?f（X（k））式2 X(k+1)=X(k)+ɑ(k)S(k) 式3 在式中，ɑ（k）为最优步长。 设与S(k）共轭的下一个方向S(k+1)由点S(k）和点X(k+1)负梯度的线性组 合构，即 S (k+1)=-?f （X （k+1））+β(k)S (k) 式4 根据共轭的条件有 [S (k)]T ?2f （X (k)）S (k+1)=0 式5 把式2和式4带入式5，得 -[?f(X (k))]T ?2f （X (k)）[-?f （X （k+1））+β(k)S (k) ]=0 式6 对于式1，则在点X (k)和点X (k+1)的梯度可写为 ?f(X (k))=HX (k)+b 式7 ?f （X （k+1）)=HX (k+1)+b 式8 把上面两式相减并将式3代入得 ɑ(k)H S (k)=?f （X （k+1）)-?f(X (k)) 式9 将式4和式9两边分别相乘，并代入式5得 -[?f （X （k+1）)+β(k)?f(X (k))]T [?f （X （k+1）)-?f(X (k)]=0 式10 将式10展开，并注意到相邻两点梯度间的正交关系，整理后得 β （k ） =2 2 ||))((||||))1((||k X f k X f ?+? 式11 把式11代入式4和式3，得 S （k+1）=-?f （X (k)）+β （k ） S （k ） X (k+1)=X （k ）+ɑ（k ）S （k ） 由上可见，只要利用相邻两点的梯度就可以构造一个共轭方向。以这种方式产生共轭方向并进行迭代运算的方法，即共轭梯度法。 最优化方法第四次作业 题目：利用FR-共轭梯度法求解无约束优化问题222 12122min ()44412x R f x x x x x x ∈=+--。初始点(0)(0.5,1).T x =- () ()T k k T k k k k k k k g g g g k d g k g d 111110.0,;0,-----=???≥+-=-=ββ 一、程序 function [x,val,k]=frcg(fun,gfun,x0) %功能：用FR 共轭梯度法求解无约束问题min f （x ） %输入：x0是初始点，fun,gfun 分别是求目标函数和梯度 %输出：x,val 分别是近似最优点和最优值，k 是迭代次数 maxk=5000; rho=0.6; sigma=0.4; k=0; epsilon=1e-4; n=length(x0); while (k while (m<20) %用Armijo 搜索求步长 if (feval(fun,x0+rho^m*d) 车辆薄板有限元分析中的多因子不完全分解 预处理解法 姚 松 田红旗 1中南大学轨道交通安全教育部重点实验室，湖南长沙，410075 dynacn@https://www.360docs.net/doc/d813008790.html, 摘要：薄板是轨道车辆结构的主要形式，本文基于离散Kirchhoff 假设的DKT 弯曲板单元推导了四边形弯曲板单元DKQ 的构造过程，并进一步阐述了用于一般薄板问题分析的平板单元的构造。提出了一种“多因子不完全分解” 的预处理方法，与共轭梯度迭代法结合能够大大加快薄板问题大型稀疏方程组的收敛速度，经过数值试验，说明该方法是稳定可靠的。该方法避免了常规不完全分解不适用于薄板这样的 “病态”结构的情况。在此基础上，编写了一般薄板问题分析的有限元程序，程序对结构刚度矩阵采用压缩存贮的方法，节约了大量内存空间。本文还对分解算法中的可选参数进行了优化研究。通过一个数值试验，本程序计算结果与商业有限元软件ANSYS5.7的结果完全一致。 关键词：薄板结构，DKQ 单元，预处理，不完全分解，共轭梯度法 1 概 述 有限元单元法已经成为结构分析的重要方法，薄板结构是轨道车辆的主要结构形式，因此薄板结构有限元分析已成为车辆结构分析中的重大课题。早期的弯曲板单元大多基于经典的薄板理论，在以该理论为基础的板单元的能量泛函中，包含位移的二阶偏导数，要求位移为类连续。这给构造板单元带来了困难，由此研究人员将注意力转向了中厚板单元，大多采用中厚板理论，其能量泛函仅包含位移的一阶导数，只要求位移是类连续，但是用厚板理论建立的单元仅对中厚板有效，当板逐渐变薄时，单元刚度矩阵中的剪切项占主导地位，计算出的弯曲变形远小于实际变形；当板非常薄时，求得的位移趋向于零，从而产生了“剪切闭锁”现象。 1C issner Mindlin Re ?0C 基于离散的假设， 通过挠度和转角分别独立插值，然后在若干个离散点上强迫挠度与转角满足薄板经典理论中的约束，构造出三角形（DKT ）和四边形（DKQ ）薄板弯曲单元，其泛函的表达式又回复为经典薄板理论的泛函表达式，又自然解决了“剪切闭锁现象”问题。多个文献表明DKT 元与DKQ 元在求解薄板弯曲问题时都显示出良好的性能，具有较高的精度。在对实际车辆结构进行有限元分析时，由于结构受力复杂，在承受板平面内的载荷的同时，也有可能板平面外的载荷，因此在进行分析时所采用的平板单元是平面应力单元与DKT 弯曲单元的组合而成。由于三角形平面应力单元为常应变单元，为了提高分析的精度，在本文中我们讨论由四边形膜单元和DKQ 单元组合而成的平板单元。 Kirchhoff Kirchhoff 1 教育部博士点基金（20020533007）项目资助 1https://www.360docs.net/doc/d813008790.html, 编号:_ 09 《最优化方法》 课程设计 题目：共轭梯度算法分析与实现 院系：数学与计算科学学院 专业：数学与应用数学 姓名学号： 指导教师： 日期：2013 年12 月23 日 摘要 在最优化计算中，共轭梯度法是非常重要的一种方法。共轭梯度法是一种改进的最速下降法，介于最速下降法与牛顿法之间的一种无约束优化算法,是为求解目标函数为二次函数的问题而设计的一类算法。它利用目标函数的梯度逐步产生共轭方向并将其作为搜索方向的方法，收敛速度快。共轭梯度法仅需利用一阶导数信息，避免了牛顿法需要存储和计算Hesse矩阵并求逆的缺点，具有二次终止性。 关键词：共轭梯度法；牛顿法；二次函数；无约束优化 Abstract In the calculation of optimization method, conjugate gradient method is a very important one. The conjugate gradient method is a unconstrained optimization method between the steepest descent method and Newton method, and sove the objective function for the original quadratic function problems and design for a class of algorithm. Conjugate gradient method using only first derivative information, to avoid the Newton method requires storage and computing the inverse Hesse matrix and shortcomings, this method has the quadratic termination. Keywords: Conjugate gradient method; Newton method;Unconstrained optimization 万方数据 万方数据 大型复线性方程组预处理双共轭梯度法 作者：张永杰， 孙秦， ZHANG Yong-jie， SUN Qin 作者单位：西北工业大学航空学院,西安,710072 刊名： 计算机工程与应用 英文刊名：COMPUTER ENGINEERING AND APPLICATIONS 年，卷(期)：2007,43(36) 被引用次数：1次 参考文献(6条) 1.Wu J P;Wang z H;Li X M High-performanco solution and paral lel computation of sparse linear equations 2004 2.Xu S F Theories and Methods on matrix computations 1995 3.Jin J M The finite element method in electromagneties 2002 4.Zhang Yongjie;Sun Qin A new ICCG method of large scale sparse linear equations[期刊论文]-Journal on Numerical Methods and Computer Applications 2007(02) 5.Betmmens Robert Itemtive solution methods 2004 6.Wu J P;Wang Z H Problems and improvements to the incomplete Cholesky decomposition with thresholds [期刊论文]-Journal on Numerical Methods and Computer Applications 2003(03) 引证文献(1条) 1.明星.苑秉成.刘建国基于共轭梯度的宽带相关处理快速算法[期刊论文]-系统工程与电子技术 2010(12) 本文链接：https://www.360docs.net/doc/d813008790.html,/Periodical_jsjgcyyy200736007.aspx M a t l a b PDE工具箱有限元法求解偏微分 方程 在科学技术各领域中，有很多问题都可以归结为偏微分方程问题。在物理专业的力学、热学、电学、光学、近代物理课程中都可遇见偏微分方程。 偏微分方程，再加上边界条件、初始条件构成的数学模型，只有在很特殊情况下才可求得解析解。随着计算机技术的发展，采用数值计算方法，可以得到其数值解。 偏微分方程基本形式 而以上的偏微分方程都能利用PDE工具箱求解。 PDE工具箱 PDE工具箱的使用步骤体现了有限元法求解问题的基本思路，包括如下基本步骤： 1) 建立几何模型 2) 定义边界条件 3) 定义PDE类型和PDE系数 4) 三角形网格划分 5) 有限元求解 6) 解的图形表达 以上步骤充分体现在PDE工具箱的菜单栏和工具栏顺序上，如下 具体实现如下。 打开工具箱 输入pdetool可以打开偏微分方程求解工具箱，如下 首先需要选择应用模式，工具箱根据实际问题的不同提供了很多应用模式，用户可以基于适当的模式进行建模和分析。 在Options菜单的Application菜单项下可以做选择，如下 或者直接在工具栏上选择，如下 列表框中各应用模式的意义为： ① Generic Scalar：一般标量模式（为默认选项）。 ② Generic System：一般系统模式。 ③ Structural Mech.，Plane Stress：结构力学平面应力。 ④ Structural Mech.，Plane Strain：结构力学平面应变。 ⑤ Electrostatics：静电学。 ⑥ Magnetostatics：电磁学。 ⑦ Ac Power Electromagnetics：交流电电磁学。 ⑧ Conductive Media DC：直流导电介质。 ⑨ Heat Tranfer：热传导。 ⑩ Diffusion：扩散。 可以根据自己的具体问题做相应的选择，这里要求解偏微分方程，故使用默认值。此外， 对于其他具体的工程应用模式，此工具箱已经发展到了Comsol Multiphysics软件，它提供了更强大的建模、求解功能。 另外，可以在菜单Options下做一些全局的设置，如下 l Grid：显示网格 l Grid Spacing…：控制网格的显示位置 l Snap：建模时捕捉网格节点，建模时可以打开 l Axes Limits…：设置坐标系范围 l Axes Equal：同Matlab的命令axes equal命令 建立几何模型 使用菜单Draw的命令或使用工具箱命令可以实现简单几何模型的建立，如下 油藏数值模拟技术现状与发展趋势 摘要：介绍了当前国内外油藏数值模拟的现状，简述了并行算法、网格技术、粗化技术、数值解法、动态油藏模型建立、动态跟踪模拟及三维显示等技术，指出了数值模拟的发展趋势。 关键词：并行算法；网格技术；网格粗化；分阶段模拟；动态跟踪模拟；数值解法 引言 近年来，随着计算机、应用数学和油藏工程学科的不断发展，油藏数值模拟方法得到不断的改进和广泛应用。通过数值模拟可以搞清油藏中流体的流动规律、驱油机理及剩余油的空间分布；研究合理的开发方案，选择最佳的开采参数,以最少的投资,最科学的开采方式而获得最高采收率及最大经济效益[1]。经过几十年的发展,该技术不断成熟和完善并呈现出一些新的特点。 1 国内外现状 1.1 并行算法 并行算法是一些可同时执行的诸进程的集合这些进程互相作用和协调动作从而达到给定问题的求解[2]。并行算法首先需合理地划分模块,其次要保证对各模块的正确计算,再次为各模块间通讯安排合理的结构,最后保证各模块计算的综合效果并行机及并行软件的开发和应用将极大地提高运算速度,以满足网格节点不断增多的油藏数值模型。在并行计算机上使用并行数值解法是提高求解偏微分方程的计算速度,缩短计算时间的一个重要途径[3,4]。在共享内存的并行机上把一个按向量处理的通用油藏模拟器改写成并行处理是容易的,但硬件扩充难；分布内存并行机编程较共享式并行机困难,但硬件扩充容易,关键是搞好超大型线形代数方程组求解的并行化。并行部分包括输入输出、节点物性、构造矩阵、节点流动及井筒等。 1.2 网格技术 为了模拟各种复杂的油藏、砂体边界或断层渗透率在垂向或水平方向的各向异性,以及近井地区的高速、高压力梯度的渗流状态,近年来在国外普遍发展了各种类型的局部网格加密及灵巧的网格技术。这种系统大体可以分为二类：一类称控制体积有限元网格(CVFE),这是将油藏按一定规则剖分为若干个三角形以后,把三角形的中心和各边的中点连接起来所形成的网格。另一类则称垂直等分线排比网格(PEBI),其剖分方法是将油藏分成若干三角形后,使三角形各边的垂直等分线相交而形成网格。这些方法在处理复杂几何形状油藏及进行局部网格加密时简单而一致。在多相流情况下,参照某一给定的几何准则时该方法是单调的,这保证了其稳定性和收敛性。这两种方法都能以直观的控制体积的概念出发并且采用一致的上游权而推导得出这些方法对网格的方向不敏感,在某些情况下比九点差分格式的效果好。 1.3 计算机辅助历史拟合技术 最速下降法and 共轭梯度法 分类：matlab 2011-04-17 17:02 961人阅读评论(2) 收藏举报算法出版优化 注明：程序中调用的函数jintuifa.m golddiv.m我在之前的笔记中已贴出 最速下降法 %最速下降法求解f = 1/2*x1*x1+9/2*x2*x2的最小值，起始点为x0=[9 1] %算法根据最优化方法（天津大学出版社）97页算法3.2.1编写 %v1.0 author: liuxi BIT %format long syms x1 x2 alpha; f = 1/2*x1*x1+9/2*x2*x2;%要最小化的函数 df=jacobian(f,[x1 x2]);%函数f的偏导 epsilon=1e-6;%0.000001k=0; x0=[9 1];%起始点 xk=x0; gk=subs(df,[x1 x2],xk);%起始点的梯度 gk=double(gk); while(norm(gk)>epsilon)%迭代终止条件||gk||<=epsilon pk=-gk;%负梯度方向 f_alpha=subs(f,[x1 x2],xk+alpha*pk);%关于alpha的函数 [left right] = jintuifa(f_alpha,alpha);%进退法求下单峰区间 [best_alpha best_f_alpha]=golddiv(f_alpha,alpha,left,right);%黄金分割法求最优步长xk=xk+best_alpha*pk; k=k+1; gk=subs(df,[x1 x2],xk); gk=double(gk); end best_x=xk;%最优点 best_fx=subs(f,[x1 x2],best_x);%最优值 共轭梯度法 最速下降法 最速下降法又称为梯度法，是1847 年由著名数学家Cauchy 给出的，它是解析法中最古老的一种，其他解析方法或是它的变形，或是受它的启发而得到的，因此它是最优化方法的基础。作为一种基本的算法，他在最优化方法中占有重要地位。其优点是工作量少，存储变量较少，初始点要求不高；缺点是收敛慢，效率不高，有时达不到最优解。非线性规划研究的对象是非线性函数的数值最优化问题。它的理论和方法渗透到许多方面，特别是在军事、经济、管理、生产过程自动化、工程设计和产品优化设计等方面都有着重要的应用。而最速下降法正是n元函数的无约束非线性规划问题min f (x)的一种重要解析法，研究最速下降法原理及其算法实现对我们有着极其重要的意义。 最速下降法 1.最速下降方向 函数f(x)在点x处沿方向d的变化率可用方向导数来表示。对于可微函数，方向导数等于梯度与方向的内积，即： Df(x;d) = ▽f(x)Td, 因此，求函数f(x)在点x处的下降最快的方向，可归结为求解下列非线性规划：min ▽f(x)Td s.t. ||d|| ≤ 1 当 d = -▽f(x) / ||▽f(x)|| 时等号成立。因此，在点x处沿上式所定义的方向变化率最小，即负梯度方向为最速下降方向。 2.最速下降算法 最速下降法的迭代公式是 x(k+1) = x(k) + λkd(k) , 其中d(k)是从x(k)出发的搜索方向，这里取在x(k)处的最速下降方向，即 d = -▽f(x(k)). λk是从x(k)出发沿方向d(k)进行一维搜索的步长，即λk满足 f(x(k) + λkd(k)) = min f(x(k)+λd(k)) (λ≥0). 计算步骤如下： (1)给定初点x(1) ∈ Rn，允许误差ε> 0，置k = 1。 (2)计算搜索方向d = -▽f(x(k))。 (3)若||d(k)|| ≤ε，则停止计算；否则，从x(k)出发，沿d(k)进行一维搜索，求λk，使 f(x(k) + λkd(k)) = min f(x(k)+λd(k)) (λ≥0). (4)令x(k+1) = x(k) + λkd(k) ，置k = k + 1，转步骤(2)。 https://www.360docs.net/doc/d813008790.html,/yangxi/archive/2011/10/20/2219408.html 梯度下降法 牛顿法 迭代公式：(1)2()1()[()]()k k k k x x f x f x +-=-?? Matlab 代码： function [x1,k] =newton(x1,eps) hs=inline('(x-1)^4+y^2'); 写入函数 ezcontour(hs,[-10 10 -10 10]); 建立坐标系 hold on; 显示图像 syms x y 定义变量 f=(x-1)^4+y^2; 定义函数 grad1=jacobian(f,[x,y]); 求f 的一阶梯度 grad2=jacobian(grad1,[x,y]); 求f 的二阶梯度 k=0; 迭代初始值 while 1 循环 grad1z=subs(subs(grad1,x,x1(1)),y,x1(2)); 给f 一阶梯度赋初值 grad2z=subs(subs(grad2,x,x1(1)),y,x1(2)); 给f 二阶梯度赋初值 x2=x1-inv(grad2z)*(grad1z)'; 核心迭代公式 if norm(x1-x2) end end end 优点：在极小点附近收敛快 缺点：但是要计算目标函数的hesse 矩阵 最速下降法 1. ：选取初始点xo ，给定误差 2. 计算一阶梯度。若一阶梯度小于误差，停止迭代，输出 3. 取()()()k k p f x =? 4. 10 t ()(), 1.min k k k k k k k k k k t f x t p f x tp x x t p k k +≥+=+=+=+进行一维搜索，求，使得令转第二步 例题： 求min (x-2)^4+(x-2*y)^2.初始值（0,3）误差为0.1 (1)编写一个目标函数，存为f.m function z = f( x,y ) z=(x-2.0)^4+(x-2.0*y)^2; end (2)分别关于x 和y 求出一阶梯度，分别存为fx.m 和fy.m function z = fx( x,y ) z=2.0*x-4.0*y+4.0*(x-2.0)^3; end 和 function z = fy( x,y ) 共轭梯度法及其基本性质 预备知识 定义1 设是对称正定矩阵。称是A-共轭的，是指 性质1 设有是彼此共轭的维向量，即 则一定是线性无关的。 ［证明］若有一组数满足 则对一切一定有 注意到，由此得出：即所有的＝０．因此， 是线性无关的． 性质２设向量是线性无关的向量组，则可通过它们的线性组合得出一组向量，而是两两共轭的． ［证明］我们用构造法来证实上面的结论． Ｓ０：取； Ｓ１：令，取． …… Ｓm：令 取 容易验证：符合性质２的要求． 性质３设是两两Ａ－共轭的，是任意指定的向量，那么从出发，逐次沿方向搜索求的极小值，所得序列，满足： ． ［证明］由下山算法可知，从出发，沿方向搜索，获得 从而 性质４设是两两Ａ共轭的，则从任意指定的出发，依次沿搜索，所得序列满足： （１） （２），其中是方程组(5.1.1)的解． ［证明］（１）是性质３的直接推论，显然成立． （２）由于是两两Ａ共轭的，故是线性无关的．所以对于向量可用线性表出，即存在一组数使 由于及，得出 ， 于是，再由得出 于是，与得出一样地，我们可以陆续得出： 对比和的表达式可知， 证明完毕 性质４是性质３的直接推论．但它给出了一种求（５.１.１）的算法，这种算法称之为共轭方向法．结合性质２，我们可以得到如下的性质５． 性质５设是上的一组线性无关的向量，则从任意指定的出发，按以下迭代产生的序列： Ｓ１：取，，； Ｓ２：计算，取； 计算，得出； 如此进行下去，直到第n步： Ｓn：计算取 计算，得出． 显然： 根据性质４可知，不论采用什么方法，只要能够构造个两两Ａ共轭的向量作为搜索方向，从任一初始向量出发，依次沿两两Ａ共轭的方向进行搜索，经 步迭代后，便可得到正定方程组的解． 最优化课程设计 题目：共轭梯度法 姓名：田鑫 指导老师：智红英 学号： 201118030216 班级：信息与计算科学111802班 共轭梯度法（Conjugate Gradient） 是介于最速下降法与牛顿法之间的一个方法，它仅需利用一阶导数信息，但克服了最速下降法收敛慢的缺点，又避免了牛顿法需要存储和计算Hesse矩阵并求逆的缺点，共轭梯度法不仅是解决大型线性方程组最有用的方法之一，也是解大型非线性最优化最有效的算法之一。在各种优化算法中，共轭梯度法是非常重要的一种。其优点是所需存储量小，具有步收敛性，稳定性高，而且不需要任何外来参数。 设我们要求解下列线性系统 其中n-×-n矩阵A是对称的（也即，A T = A），正定的（也即，x T Ax > 0对于所有非0向量x属于R n），并且是实系数的。 将系统的唯一解记作x*。 最后算法 经过一些简化，可以得到下列求解Ax = b的算法，其中A是实对称正定矩阵。x := 0 k := 0 r := b repeat until r k is "sufficiently small": k := k + 1 if k = 1 p := r0 1 else end if x := x k-1 + αk p k k r := r k-1 - αk A p k k end repeat 结果为x k 共轭梯度法程序源代码 #include 最优化方法实验报告optimization method Experiment Report 学生所在学院：理学院 学生所在班级：信息1 学生姓名： 教务处 20014年5 月 最优化方法实验报告书 说明:1.下面程序在MATLAB R2012a 中均能正常运行。 2.程序之间有关联。 实验一熟悉MATLAB基本功能（２学时）实验的目的和要求：在本次实验中，通过亲临使用MATLAB，对该软件做一全面了解并掌握重点内容。 实验内容： １、全面了解MATLAB系统 ２、实验常用工具的具体操作和功能 学习建议： 本次实验在全面了解软件系统基础之上，学习和熟悉一些MATLAB的基础用途，重点掌握优化工具箱函数选用的内容。 重点和难点： 优化工具箱函数选用。 数学模型： 其中f,x,b,beq,lb和ub为向量，A和Aeq为矩阵。 语法： x = linprog(f,A,b,Aeq,beq) x = linprog(f,A,b,Aeq,beq,lb,ub) x = linprog(f,A,b,Aeq,beq,lb,ub,x0) x = linprog(f,A,b,Aeq,beq,lb,ub,x0,options) [x,fval] = linprog(...) [x,fval,exitflag] = linprog(...) [x,fval,exitflag,output] = linprog(...) [x,fval,exitflag,output,lambda] = linprog(...) 描述： x = linprog(f,A,b)求解问题 min f'*x，约束条件为A*x <= b。 x = linprog(f,A,b,Aeq,beq)求解上面的问题，但增加等式约 束，即Aeq*x = beq。若没有不等式存在，则令A=[]、b=[]。最优化课程设计--共轭梯度法算法分析与实现
Matlab PDE工具箱有限元法求解偏微分方程
16梯度法和共轭梯度法基本原理和特点
微分方程数值解--大纲
共轭梯度法
共轭梯度法C语言(西安交大)
共轭梯度法程序
作业4-FR共轭梯度法
车辆薄板有限元分析中的多因子不完全分解预处理解法
共轭梯度算法分析与实现
大型复线性方程组预处理双共轭梯度法
MatlabPDE工具箱有限元法求解偏微分方程教学提纲
油藏数值模拟技术现状与发展趋势
共轭梯度法
消去法实验报告19
最优化牛顿法最速下降法共轭梯度法matlab代码
共轭梯度法及其基本性质
共轭梯度法c++程序
最优化方法实验报告