共轭梯度法matlab实例
运筹学实验报告(F-R共轭梯度法、Wolfe简约梯度法)
一、实验目的:1、掌握求解无约束最优化问题的 F-R 共轭梯度法，以及约束最优化问题 Wolfe 简约梯度法。
2、学会用MATLAB 编程求解问题，并对以上方法的计算过程和结果进行分析。
二、实验原理与步骤: 1、F-R 共轭梯度法基本步骤是在点)(k X 处选取搜索方向)(k d , 使其与前一次的搜索方向)1(-k d关于A 共轭，即(1)()(1),0k k k d d Ad --<>=然后从点)(k X 出发，沿方向)(k d 求得)(X f 的极小值点)1(+k X ,即)(m in )()()(0)1(k d X f X f k k λλ+=>+如此下去, 得到序列{)(k X }。
不难求得0,)1()(>=<-k k Ad d的解为)()1()1()()()()1(,,k k k k k k k d Ad d d AX b XX><>-<+=--+注意到)(k d 的选取不唯一，我们可取)1(1)()()(--+-∇=k k k k d X f d β由共轭的定义0,)1()(>=<-k k Add 可得： ><><-=----)1()1()1()(1,,k k k k k Ad d Ad r β共轭梯度法的计算过程如下：第一步：取初始向量)0(X , 计算⎪⎪⎩⎪⎪⎨⎧+=><><-=-=-∇==(0)0(0)(1))0()0()0()0(0(0)(0)(0)(0)d X X ,,X )X (r d λλAd d Ad r A b f第1+k 步：计算⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧+=><><-=+=><><-=-=-∇=+------(k)0(k)1)(k )()()()()1(1(k))()1()1()1()(1(k)(k)(k)d X X ,,r ,,X )X (r λλββk k k k k k k k k k k k k Ad d Ad r d d Ad d Adr A b f2、Wolfe 简约梯度法Wolfe 基本计算步骤：第一步：取初始可行点 x 0∈X l ,给定终止误差ε>0 ，令k:=0；第二步：设 I B k是x k 的 m 个最大分量的下标集，对矩阵A 进行相应分解 A =(B k ,N k )；第三步：计算 ∇f(x k)=(∇B f(x k )∇Nf(x k )) ,然后计算简约梯度r N k=−(B k −1N k )T ∇B f(x k )+∇N f(x k )；第四步：构造可行下降方向 p k . 若||p k ||≤ε ，停止迭代，输出x k 。
求解无约束优化问题的共轭梯度法
求解无约束优化问题的共轭梯度法李芳梅，姚瑞哲指导教师：李良摘要：本文主要针对无约束优化问题，利用共轭梯度法（CG方法）求解此类问题，并得出其迭代次数及问题的解。
论文对此种方法给出了具体事例，并对例子进行了matlab软件实现。
1.引言共轭梯度法时介于最速下降法与牛顿法之间的一个方法。
它仅需利用一阶导数信息，但克服了最速下降法收敛慢的缺点，又避免了牛顿法需要存储和计算Hesse矩阵并求逆的缺点。
共轭梯度法最早是由Hestenes和Stiefel（1952）提出来的，用于解正定系数矩阵的线性方程组。
由于共轭梯度法不需要矩阵存储，且有较快的收敛速度和二次终止性等优点，现已广泛应用与实际问题中。
2.共轭梯度法共轭梯度法是共轭方向法的一种。
所谓共轭方向法就是其所有的搜索方向都是相互共轭的方法。
定义设G是n×n对称正定矩阵，d1，d2是n维非零向量。
如果d1T G d2=0，则称向量d1和d2是G-共轭的，简称共轭的。
设d1，d2，…，d m是R n中任一组非零向量，如果d i T G d j=0 （i≠j），则称d1，d2，…，d m是G-共轭的，简称共轭的。
显然，如果d1，d2，…，d m是共轭的，则它们是线性无关的。
算法（共轭梯度法）步1.初始步：给出x0，ε>0,计算r0=G x0-b，令d0=-r0，k:=0. 步2.如果‖r k‖≤ε，停止.步3.计算αk=r k T r k/d k T Gd k,x k+1=x k+αk d k，r k+1=r k+αk Gd k，βk=r k+1T r k+1/r k T r k，d k+1=-r k+1+βk d k.步4.令k:=k+1，转步2.1.共轭梯度法的matlab实现（数值例子）首先建立如下的m.文件function[x,iter]=cgopt(G,b,x0,max_iter,tol)x=x0;fprintf('\n x0=');fprintf('%10.6f',x0);r=G*x-b;%残差d=-r;for k=1:max_iterif norm(r)<=toliter=k-1;fprintf('\n Algorithm finds a solution!');return endalpha=(r'*r)/(d'*G*d);%收敛速度 xx=x+alpha*d; rr=r+alpha*G*d; beta=(rr'*rr)/(r'*r); d=-rr+beta*d; x=xx; r=rr;fprintf('\n x%d=',k); fprintf('%10.6f',x); enditer=max_iter; return然后建立CG_main 的m.文件，带入数值例子function CG_main()12345123451234512345123451023412923277351432312174351512x x x x x x x x x x x x x x x x x x x x x x x x x ++++=⎧⎪+-+-=-⎪⎪-++-=⎨⎪+++-=-⎪---+=⎪⎩ G=[10 1 2 3 4;1 9 -1 2 -3;2 -1 7 3 -5;3 2 3 12 -1;4 -3 -5 -1 15]; b=[12 -27 14 -17 12]'; x0=[0 0 0 0 0]'; max_iter=1000; tol=1e-6;fprintf('\n');fprintf('Conjugate Gradiant Method:\n'); fprintf('=================\n');[x,iter]=cgopt(G,b,x0,max_iter,tol);fprintf('Iterative number:\n %d\n',iter); fprintf('Solution:\n'); fprintf('%10.6f',x);fprintf('\n================\n');4.结论实际上，共轭梯度法是最速下降法的一种改进。
matlab梯度算法
matlab梯度算法Matlab梯度算法在数学和计算机科学中，梯度是指一个多元函数在某一点上的变化率或斜率。
梯度算法是一种优化算法，用于找到函数的最小值或最大值。
在Matlab中，有多种方法可以使用梯度算法来优化函数，包括梯度下降和共轭梯度法。
本文将详细介绍Matlab中的梯度算法，并逐步讲解其原理和应用。
I. 梯度下降法梯度下降法是一种基于迭代的优化算法，通过计算函数的梯度来更新参数的值，以逐步接近函数的最小值。
在Matlab中，可以使用"gradientDescent"函数来实现梯度下降法。
1. 实现梯度下降法首先，我们需要定义一个优化目标函数，例如：f(x) = x^2 + 2x + 1。
然后，定义其梯度函数为g(x) = 2x + 2。
接下来，我们可以使用以下代码来计算梯度下降：matlab定义优化目标函数f = (x) x^2 + 2*x + 1;定义梯度函数g = (x) 2*x + 2;初始化参数x0 = 0;设置学习率和迭代次数alpha = 0.01;iterations = 100;梯度下降法for i = 1:iterationsx0 = x0 - alpha * g(x0);end打印最优解disp(['Optimal solution: ', num2str(x0)]);在这个例子中，我们使用了学习率(alpha)为0.01，迭代次数(iterations)为100。
通过不断更新参数x0的值，最终得到了最优解。
2. 梯度下降法的原理梯度下降法的核心思想是利用函数在当前点的梯度信息来更新参数的值，以便能够向着函数的最小值前进。
具体来说，算法的步骤如下：a. 初始化参数的值：选择一个初始参数的值作为起始点。
b. 计算梯度：计算函数在当前点的梯度，即求解函数关于参数的偏导数。
c. 更新参数：根据当前点的梯度和学习率，通过减去梯度的乘积来更新参数的值。
最优化方法实验报告(2)
最优化方法实验报告Numerical Linear Algebra And ItsApplications学生所在学院：理学院学生所在班级：计算数学10-1学生姓名：甘纯指导教师：单锐教务处2013年5月实验三实验名称：无约束最优化方法的MATLAB实现实验时间: 2013年05月10日星期三实验成绩：一、实验目的：通过本次实验的学习，进一步熟悉掌握使用MATLAB软件，并能利用该软件进行无约束最优化方法的计算。
二、实验背景：（一）最速下降法1、算法原理最速下降法的搜索方向是目标函数的负梯度方向，最速下降法从目标函数的负梯度方向一直前进，直到到达目标函数的最低点。
2、算法步骤用最速下降法求无约束问题n R()min的算法步骤如下：xxf，a ）给定初始点)0(x ，精度0>ε，并令k=0；b ）计算搜索方向)()()(k k x f v -∇=，其中)()(k x f ∇表示函数)(x f 在点)(k x 处的梯度；c ）若ε≤)(k v ，则停止计算；否则，从)(k x 出发，沿)(k v 进行一维搜索，即求k λ，使得)(min )()()(0)()(k k k k v x f v x f λλλ+=+≥； d ）令1,)()()1(+=+=+k k v x x k k k k λ，转b ）。
（二）牛顿法1、算法原理牛顿法是基于多元函数的泰勒展开而来的，它将)()]([-)(1)(2k k x f x f ∇∇-作为搜索方向，因此它的迭代公式可直接写出来：)()]([)(1)(2)()(k k k k x f x f x x ∇∇-=-2、算法步骤用牛顿法求无约束问题n R x x f ∈),(min 的算法步骤如下：a ）给定初始点)0(x ,精度0>ε，并令k=0；b ）若ε≤∇)()(k x f ，停止，极小点为)(k x ，否则转c ）；c ）计算)()]([,)]([)(1)(2)(1)(2k k k k x f x f p x f ∇∇-=∇--令；d ）令1,)()()1(+=+=+k k p x x k k k ，转b ）。
共轭梯度法matlab程序
共轭梯度法matlab程序共轭梯度法是一种用于求解无约束优化问题的迭代算法。
以下是一个简单的MATLAB 程序示例，用于实现共轭梯度法：matlab复制代码function[x, fval, iter] = conjugate_gradient(A, b, x0, tol,max_iter)% A: 矩阵% b: 向量% x0: 初始解% tol: 容忍度% max_iter: 最大迭代次数r = b - A*x0;p = r;x = x0;fval = A*x - b;iter = 0;while (norm(r) > tol) && (iter < max_iter)Ap = A*p;alpha = dot(p, r) / dot(p, Ap);x = x + alpha*p;r = r - alpha*Ap;if iter == 0fval_new = dot(r, r);elsebeta = dot(r, r) / dot(p, Ap);p = r + beta*p;fval_new = dot(r, r);endif fval_new < fvalfval = fval_new;enditer = iter + 1;endend该程序接受一个矩阵A、一个向量b、一个初始解x0、一个容忍度tol和一个最大迭代次数max_iter作为输入，并返回最终解x、目标函数值fval和迭代次数iter。
使用该函数时，您需要将要优化的目标函数转换为矩阵形式，并调用该函数来找到最优解。
例如，假设您要最小化一个函数f(x)，可以将该函数转换为矩阵形式A*x - b，其中A是目标函数的雅可比矩阵，b是目标函数的常数项向量。
然后，您可以调用该函数来找到最优解，例如：matlab复制代码A = jacobian(f); % 计算目标函数的雅可比矩阵b = [1, 2, 3]; % 目标函数的常数项向量x0 = [0, 0, 0]; % 初始解向量tol = 1e-6; % 容忍度max_iter = 1000; % 最大迭代次数[x, fval, iter] = conjugate_gradient(A, b, x0, tol, max_iter); % 调用共轭梯度法函数求解最优解。
共轭梯度法实验报告
数值代数实验报告一、实验名称：用共轭梯度法解线性方程组。
二、实验目的：进一步熟悉理解掌握共轭梯度法解法思路，提高matlab 编程能力。
三、实验要求：已知线性方程矩阵，应用共轭梯度法在相关软件编程求解线性方程组的解。
四、实验原理：1．共轭梯度法：考虑线性方程组Ax b =的求解问题，其中A 是给定的n 阶对称正定矩阵，b 是给定的n 维向量，x 是待求解的n 维向量.为此,定义二次泛函()2T T x x Ax b x ϕ=-.定理1 设A 对称正定,求方程组Ax b =的解，等价于求二次泛函()x ϕ的极小值点. 定理1表明，求解线性方程组问题就转化为求二次泛函()x ϕ的极小值点问题.求解二次函数极小值问题，通常好像盲人下山那样，先给定一个初始向量0x ，确定一个下山方向0p ，沿着经过点0x 而方向为0p 的直线00x x p α=+找一个点1000x x p α=+，使得对所有实数α有()()00000x p x p ϕαϕα+≤+，即在这条直线上1x 使()x ϕ达到极小.然后从1x 出发，再确定一个下山的方向1p ，沿着直线11x x p α=+再跨出一步，即找到1α使得()x ϕ在2111x x p α=+达到极小：()()11111x p x p ϕαϕα+≤+.重复此步骤，得到一串012,,,ααα 和 012,,,p p p ，称k p 为搜索方向，k α为步长.一般情况下，先在k x 点找下山方向k p ，再在直线k k x x p α=+上确定步长k α使()(),k k k k k x p x p ϕαϕα+≤+最后求出1k k k k x x p α+=+.然而对不同的搜索方向和步长，得到各种不同的算法.由此，先考虑如何确定k α.设从k x 出发，已经选定下山方向k p .令 ()()k k f x p αϕα=+()()()2TT k k k k k k x p A x p b x p ααα=++-+()22TT k k k k k p Ap r p x ααϕ=-+，其中k k r b Ap =-.由一元函数极值存在的必要条件有()220TT k k k k f p Ap r p αα'=-=所确定的α即为所求步长k α，即 T k kk Tk kr p p Ap α=. 步长确定后，即可算出1k k k k x x p α+=+.此时，只要0T k k r p ≠，就有()()()()1k k k k k k x x x p x ϕϕϕαϕ+-=+-()2220T k k TT k kk k k k Tk kr p p Ap r p p Ap αα=-=-<即()()1k k x x ϕϕ+<.再考虑如何确定下山方向k p .易知负梯度方向是()x ϕ减小最快的方向，但简单分析就会发现负梯度方向只是局部最佳的下山方向，而从整体来看并非最佳.故采用新的方法寻求更好的下山方向——共轭梯度法. 下面给出共轭梯度法的具体计算过程:给定初始向量0x ，第一步仍选用负梯度方向为下山方向，即00p r =，于是有00010001000,,T T r r x x p r b Ax p Ap αα==+=-.对以后各步，例如第k+1步（k ≥1）,下山方向不再取k r ，而是在过点由向量k r 和1k p -所张成的二维平面21{|,,}k k k x x x r p R πξηξη-==++∈内找出使函数ϕ下降最快的方向作为新的下山方向k p .考虑ϕ在2π上的限制：()1,()k k k x r p ψξηϕξη-=++11()()T k k k k k k x r p A x r p ξηξη--=++++12()T k k k b x r p ξη--++. 计算ψ关于,ξη的偏导得： ()()11112,2,T T T k k k k k kT T k k k k r Ar r Ap r r r Ap p Ap ψξηξψξηη----∂=+-∂∂=+∂其中最后一式用到了10T k k r p -=，这可由k r 的定义直接验证.令 0ψψξη∂∂==∂∂， 即知ϕ在2π内有唯一的极小值点001k k k x x r p ξη-=++，其中0ξ和0η满足 00101011,0.T T T k k k k k k T Tk k k k r Ar r Ap r r r Ap p Ap ξηξη----⎧+=⎨+=⎩由于0k r ≠必有00ξ≠，所以可取()01001k k k k p x x r p ηξξ-=-=+作为新的下山方向.显然，这是在平面2π内可得的最佳下山方向.令010k ηβξ-=，则可得1111.T k k k T k k r Ap p Ap β----=-注：这样确定的k p 满足10Tkk p Ap -=，即k p 与1k p -是相互共轭的. 总结上面的讨论，可得如下的计算公式：T k kk Tk kr p p Ap α= ， 1k k k k x x p α+=+， 11k k r b Ax ++=-，1T k kk Tk kr Ap p Ap β+=-， 11k k k k p r p β++=+. 在实际计算中，常将上述公式进一步简化，从而得到一个形式上更为简单而且对称的计算公式.首先来简化1k r +的计算公式：11()k k k k k k k k r b Ax b A x p r Ap αα++=-=-+=-.因为k Ap 在计算k α是已经求出，所以计算1k r +时可以不必将1k x +代入方程计算，而是从递推关系1k k k r b Ap α+=-得到.再来简化k α和k β的计算公式.此处需要用到关系式1110,T T T k k k k k k r r r p r p +-+=== 1,2,k =.从而可导出1111,T T k k k k r r r α+++=-， ()111TT T k k k k k k k k k p Ap p r r p r αα+=-=()1111T Tk k k k k k k kr r p r r βαα--=+=.由此可得,T k k k T k k r r p Ap α=， 11.T k k k T k kr r r r β++=.从而有求解对称正定方程组的共轭梯度法算法如下：0x =初值00r b Ax =-；0k =while 0k r ≠1k k =+ if 1k = 00p r =else21122T T k k k k k r r r r β-----= 1122k k k k p r p β----=+ end11111T Tk k k k k r r p Ap α-----=111k k k k x x p α---=+111k k k k r r Ap α---=-endk x x =注：该算法每迭代一次仅需要使用系数矩阵A 做一次矩阵向量积运算. 定理2 由共轭梯度法得到的向量组{}i r 和{}i p 具有如下基本性质： （1）0T i j p r =， 0;i j k ≤<≤ （2）0T i j r r =， i j ≠，0,;i j k ≤≤ （3）0T i j p Ap =， i j ≠，0,;i j k ≤≤ （4）000{,,}{,,}(,,1)k k span r r span p p A r k κ==+，其中0000(,,1){,,,}k A r k span r Ar A r κ+=，通常称之为Krylov 子空间.下面给出共轭梯度法全局最优性定理：定理3 用共轭梯度法计算得到的近似解k x 满足()(){}00min :(,,)k x x x x A r k ϕϕκ=∈+或{}**00min :(,,)k AAx x x x x x A r k κ-=-∈+，其中Ax=*x 是方程组Ax b =的解，0(,,)A r k κ是由所定义的Krylov 子空间. 定理2表明，向量组0,,k r r 和0,,k p p 分别是Krylov 子空间0(,,1)A r k κ+的正交基和共轭正交基.由此可知，共轭梯度法最多n 步便可得到方程组的解*x .因此，理论上来讲，共轭梯度法是直接法.然而实际使用时，由于误差的出现，使k r 之间的正交性很快损失，以致于其有限步终止性已不再成立.此外，在实际应用共轭梯度法时，由于一般n 很大，以至于迭代()O n 次所耗费的计算时间就已经使用户无法接受了.因此，实际上将共轭梯度法作为一种迭代法使用，而且通常是k r 是否已经很小及迭代次数是否已经达到最大允许的迭代次数max k 来终止迭代.从而得到解对称正定线性方程组的实用共轭梯度法，其算法如下：x =初值0;k =;r b Ax =-T r r ρ=while)()max2band k kε><1k k =+if 1k = p r = else;p r p βρρβ==+ end;;T Ap p x x p ωαρωα===+ ;;T r r r r αωρρρ=-== end算法中，系数矩阵A 的作用仅仅是用来由已知向量p 产生向量Ap ω=，这不仅可以充分利用A 的稀疏性，而且对某些提供矩阵A 较为困难而由已知向量p 产生向量Ap ω=又十分方便的应用问题是十分有益的。
用MATLAB实现最速下降法
实验的题目和要求一、所属课程名称:最优化方法二、实验日期：2０10年5月10日~2010年5月1５日三、实验目的掌握最速下降法，牛顿法和共轭梯度法的算法思想,并能上机编程实现相应的算法。
二、实验要求用ＭATLAB 实现最速下降法，牛顿法和共轭梯度法求解实例。
四、实验原理最速下降法是以负梯度方向最为下降方向的极小化算法，相邻两次的搜索方向是互相直交的。
牛顿法是利用目标函数)(x f 在迭代点k x 处的Tayl ｏr 展开式作为模型函数，并利用这个二次模型函数的极小点序列去逼近目标函数的极小点。
共轭梯度法它的每一个搜索方向是互相共轭的,而这些搜索方向k d 仅仅是负梯度方向k g -与上一次接待的搜索方向1-k d 的组合。
五．运行及结果如下：最速下降法:题目:f ＝（x-２)^２+(y －4）＾2M 文件：ｆuncti ｏｎ [R,n]=s ｔｅｅｌ(x0，y0,eps)ｓyms x ;ｓｙms y ；f=(x －2)^2+（y-4)^２;ｖ＝[x ，ｙ];j=jacobia ｎ(ｆ，v ）;T ＝[ｓu ｂs （j(1),x,x ０),ｓub ｓ（ｊ（2），y,y ０）］;t ｅmp=s ｑrt((T （1））^２+（Ｔ（2))^2);x1＝ｘ0;ｙ１=ｙ0;n=0;syms kk ；wh ｉle (t ｅmp>ep ｓ)d=－T;f1＝x １+kk*d(1);ｆ2＝y1＋kk*d(2);fＴ=[subｓ(j（1）,x,f1),subs(ｊ(2），y,f2)];fun=ｓqrｔ((ｆT(1)）＾2＋(fＴ(２))＾2)；Mini＝Gold（fun,0，1,0.0０001）;x0=x1+Ｍｉnｉ*d（1）;ｙ0=ｙ1＋Mini*d(2);T=[subs（ｊ（1）,x,x０),subs(j(2),y,y0)］；tｅmp=ｓｑrｔ((T(１))^2+(T(２))^２);x1=x0;y1=ｙ0；n=n+1;endR=[ｘ0,ｙ0]调用黄金分割法:Ｍ文件:ｆuｎctiｏnＭiｎi＝Ｇｏld(f,a０，ｂ0,eps）symｓx；format long；sｙms kk；ｕ＝ａ０+0.３82*(b0-a0）;v＝a0＋０.61８＊(b０－a0);k=0;a＝a0;b=b0;arraｙ(k+1,１）=a；arrａy(ｋ+1，2)=ｂ;while(（b-a)/(b0-ａ0)>=eps）Fu=subs（ｆ,ｋk，ｕ);Fv=ｓｕbs(f，kk,v)；if(Fｕ<＝Fv）b=v；v＝u;u=a＋0．382*(ｂ-a)；k=k+１;ｅｌｓeiｆ(Fu>Fｖ)a=u；u=ｖ;ｖ＝a+0.６18＊（b－a）；k＝k+1；endaｒray(k＋1,１)=ａ；arｒay(k＋１,2）=b;eｎdＭini=（ａ+b)/2;输入：[Ｒ,n]=sｔeel(0，1,0.０0０１)Ｒ= １．99９9941366７642 3．99９99１20501４63 Ｒ= 1.9９999４13667６42 3.９９99９1２０50１4６３n = 1牛顿法：题目:f＝(ｘ-2)^2＋(y-4)＾2M文件:syms x1x2;f＝(x１-２)＾2＋（x2-4)^２;v=[x１,ｘ２];df＝ｊacｏｂiａn（f，ｖ);df=df.＇;G＝jacｏｂｉaｎ(df,v);epｓon=1e-12；x0=［0，0]＇;g1=sｕbs（df，｛ｘ1,x2},{x0(1,1）,ｘ0(2，1)});G１=subs(G,{x1,ｘ2｝,{x0(1,1),x０（2,1)}）;k=0;ｍuｌ_coｕnt=0;sum_couｎt=0;mul_count=mｕl_cｏｕnt+12;sum_couｎt=sum_cｏunｔ+６;wｈile(norm（g1）>epson)p＝-G1\g1;x0=x0＋p;g1=subs(df，｛ｘ1，ｘ2},{ｘ0(1，１）,x０（2,1)});G1＝subｓ（G,{x1,x２},{x0（1,1），x0(２,1)}）;ｋ=k+1;mｕｌ＿ｃｏunｔ=mul_ｃounｔ+1６;suｍ＿cｏunt=sｕm_ｃounｔ+１1；ｅｎd;kｘ0mul＿couｎtsum_ｃouｎｔ结果::k ＝ 1ｘ０=24ｍuｌ_coｕnt ＝2８ｓum_ｃｏunt = 17共轭梯度法：题目：f＝(x－2)^２＋(y-4）^2M文件：funｃtion f＝conjuｇaｔｅ_gｒad_２d（x0,t）ｘ=x0;ｓyｍs ｘi yi af＝(xi-２)^2+(ｙi-４)＾２;ｆx=diff(f，xｉ）;fy=ｄifｆ(ｆ,ｙi);fx=ｓubs(fx,｛xi,yi}，x0）;fy=subs(fy，{xi，yi}，x0);fi=[ｆx,fy］;count=0；whiｌe doublｅ(ｓqrｔ(fx＾2+fy^2))>tｓ=-fi;ｉfｃount<=0s=-fｉ;elｓes=s１;eｎdｘ＝x+a*s;f=sｕｂｓ(f,{ｘｉ,yi},x);f1=diff（f)；f１=solvｅ(f1);ifｆ1~=0ａi＝double(f1);eｌsebreakx，f=sｕbs(f,{xi,yi},x),counｔendx=suｂs（x,a,aｉ）;ｆ＝xi－xi^2+2*xi*yi+yｉ^2；fxi=ｄｉff(f,ｘi);fyi＝difｆ(f,ｙi);fxｉ＝ｓｕbs(fｘi,{xi,ｙi},x）;fyｉ=subs(fyｉ，{xi，ｙi}，x);ｆｉｉ＝［fxｉ,fyi];d=(fxi＾2＋ｆyi＾２)/(fx^２+fy^2）；s1=-fｉi＋d＊s;couｎt=coｕnt+1;ｆx＝ｆxｉ；fy=ｆyi;endx,f=ｓubｓ(ｆ,{ｘｉ,yi｝,x),cｏunt输入:cｏnjugatｅ_grad_２d([０,0］,0．0００1)结果：x ＝0.2４998825499785 -０.2４9999987４1273f= 0．124９９９9９98617６cｏuｎｔ=１0aｎs= 0．1２４９9999986176六、结论如下：最速下降法越接近极小值,步长越小，前进越慢。
数值计算实例MATLAB实现(附带详细源码)
数值计算实例MATLAB实现附带详细源码1.在化学反应中，A 的一个分子和 B 的一个分子结合形成物质 C 的分子。
若在时刻t 时，物质 C 的浓度为() y t ，则其是下述初值问题的解()()() ,00y k a y b y y '=--=其中k 为正常数，a 和 b 分别表示 A 和 B 的初始浓度。
假设k = 0.01, a =70毫摩/升， b = 50 毫摩/升. 该方程的真解为0.20.2350(1)()75t te y t e---=- （1）自己编写程序，使用四阶经典Runge-Kutta （龙格-库塔法），以步长为0.5h =，在区间[0, 20]上给出() y t 的近似解； （2）列表给出真解和近似解的比较；（3）讨论当t →∞时，近似解的变化趋势，并分析该数值结果。
解：数学原理：四阶经典Runge-Kutta （龙格-库塔法）112341213243(22)6(,)(,)22(,)22(,)m m m m m m m m m m hu u k k k k k f t u h hk f t u k h hk f t u k k f t h u hk +=++++==++=++=++程序设计见附录 结果如下表：（3）近似解变化趋势当t→∞时，由以下极限方程可知：0.20.2350(1)()75lim()tttey tey t--→∞⎧-=-⎪⎨⎪⎩随着t→∞，近似值越来越接近真实值，极限的真实值为50，lim()50ty t→∞=，变化趋势也可由一下曲线图表示：感想：四阶Runge-Kutta法计算的结果精度非常好，其结果与真实解误差不大。
2.考虑定义在闭区间[−5, 5]上的函数()2112()5f x x -=+ ；（1）利用等距节点构造次数分别为 n = 4,8,16, 32 的插值多项式()n p x ，并分别画()()()()481632,,,p x p x p x p x ；（2）利用chebyshev 零点构造次数分别为 n = 4,8,16, 32 的插值多项式()n pp x()()()()481632,,,pp x pp x pp x pp x ；（3）画出当 n = 32 时，两种插值多项式的比较图，误差图，并给出相应的误差估计;（4）在这个问题中能观察到龙格现象吗？ 解：数学原理：拉格朗日插值多项式：001122()()()()()n n n L x l x y l x y l x y l x y =+++011011()()()()(),0,1,2,()()()()k k n k k k k k k k n x x x x x x x x l x k n x x x x x x x x -+-+----==----0()()()nn n in k k k k k j k jj kx x L x l x y y x x ===≠-==-∑∑∏程序设计见附录(1) 利用等距节点构造次数分别为 n = 4,8,16, 32 的插值多项式如下： ()43240.00160.00.0640.60061400p x x x x x ++=++()876542830.00280.00640.02500.02500.00640.00260.000168.001p x x x x x x x x x ++++++=++()1615141312161110987654320.00210.00280.00410.0064 60.01120.02500.09290.09290.02050 0.01120.00640.00410.002.00160180.021.000p x x x x x x x x x x x x x x x x x ++++++++++++++=++()3231302928272632252423222120191817160001600018000210002400028000340004100050006400083001120016100250004350092902906029p x .x .x .x .x .x .x .x .x .x .x .x .x .x .x .x .x x .=+++++++++++++++++151413121110987654320600929004350025000161001120008300064000500041000340002800024000210001800016x .x .x .x .x .x .x .x .x .x .x .x .x .x .x .+++++++++++++++（2）利用chebyshev 零点构造次数分别为 n = 4,8,16, 32 的插值多项式如下：()43240.00160.00320.00320.0016x x p x x p x =++++()87654328+0.00190.00320.01080.01080.00320.00196=0.0.0106001pp x x x x x x x x x +++++++()161514131211109168765432=0.0016 0.0017 0.0019 0.00230.00320.00520.01080.0403 1.00000.04030.01080.00520.00320.00230.0019 0.0017 0.0016 pp x x x x x x x x x x x x x x x x x ++++++++++++++++()323130292827263225242322212019181700016000160001700017000190002100023000270003200040000520007100108001860040301428pp x .x .x .x .x .x .x .x .x .x .x .x .x .x .x .x .x x =+++++++++++++++++16151413121110987654320142800403001860010800071000520004000320002700023000210001900017000170001600016.x .x .x .x .x .x .x .x .x .x .x .x .x .x .x .+++++++++++++++++(3)两种插值多项式的比较误差图如下(a)等距插值误差 (b) chebyshev零点插值误差(4) 等距插值在高次插值中能观察到龙格现象，而chebyshev零点插值观察不到龙格现象。