数值计算资料

3.7 多项式曲线拟合的递归最小二乘法为了便于计算机计算，下面介绍一种基于递归最小二乘法的多项式拟合方法，其原理如下：对于给定的一组数据，求作一次多项式 (3.7-1)设， (3.7-2)则式(3.7-1)改写为矩阵方程： (3.7-3)很显然，当时，式(3.7-3)是一个超定方程；当时，式(3.7-3)是一个欠定方程。为了使式(3.7-3) 的待定系数有唯一解，对式(3.7-3)进行如下变换： (3.7-4) 设，则式(3.7-4)改写为 (3.7-5)设 (3.7-6)则基于递归最小二乘法的多项式拟合方法如下：Step1 随机给定，设，给定，给定任意小正实数，由式(3.7-2)计算矩阵、，计算向量的长度，令，设置单位矩阵：，；Step2由式(3.7-6)计算多项式拟合值：；计算误差函数；计算性能指标：；Step3 计算向量：；Step4更新向量： ；Step5 更新待定系数向量：； 。如果，返回Step2重复上述计算。如果结束上述计算，否则，令返回Step2重复上述计算。为了验证该算法的有效性，给出次曲线拟合的结果如下：【例3.7】采用递归最小二乘算法进行曲线拟合，样本数据如下表所示：例3.7 样本数据表012345/1.9401.9501.9601.9701.9801.990/ 132.165151.326179.323203.302226.542249.633例3.7的源程序如下（源程序名：ex3p7.m）:clear,format longm=input(m=?); %多项式拟合曲线的次数tol=input(tol=?);J=0;A=zeros(m+1,1);x=1.940,1.950,1.960,1.970,1.980,1.990;yi=132.165,151.326,179.323,203.302,226.542,249.633; N=length(x);p=0:m;for j=1:N X(j,:)=x(j).p;endXb=X*X;Yb=X*yi;M=length(Yb);afa=1e+6;I=eye(M);Pk=afa*I;for n=1:10 % 最大迭代次数 for k=1:M % 对样本数据进行训练 fk=Xb(k,:)*A; erk=Yb(k)-fk; J=J+0.5*erk.2; Qk=Pk*Xb(k,:)./(1+Xb(k,:)*Pk*Xb(k,:); A=A+Qk*erk; % 更新待定系数向量 Pk=(I-Qk*Xb(k,:)*Pk; end if J=tol A=,A J=,J break end J=0;end yk=X*A;x1=1.940,1.950,1.960,1.970,1.980,1.990,2.000,2.010,2.020;%预测样本数据N1=length(x1);for j=1:N1 X1(j,:)=x1(j).p;endyk1=X1*A; %预测结果plot(x*1e+3,yi*1e+6,ok,x1*1e+3,yk1*1e+6,k);grid运行上述程序，当时，拟合结果分别如图3.6-3、图3.6-4和图3.6-5所示。 图3.6-3 二次多项式拟合曲线图 图3.6-4 三次多项式拟合曲线图图3.6-5 五次多项式拟合曲线图4.6数值积分的神经网络算法前面介绍了数值积分的Newton-Cotes方法、Romberg方法、Gauss方法等，其中Newton-Cotes方法是一种利用插值多项式来构造数值积分的常用方法，但是高阶的Newton-Cotes方法的收敛性没有保证，因此，在实际计算中很少使用高阶的Newton-Cotes公式；Romberg方法收敛速度快、计算精度较高，但是计算量较大；Gauss方法积分精度高、数值稳定、收敛速度较快，但是节点与系数的计算较麻烦，而且高斯点往往是无理数，在计算时常有舍入误差。本节介绍的数值积分计算方法，其基本思想是训练余弦基函数神经网络来逼近被积函数以实现定积分的数值计算。由于余弦基函数易于积分，因此，任意被积函数的积分计算可近似表示为余弦基函数的积分计算。研究结果表明，基于神经网络算法的数值积分方法有高的计算精度、收敛速度快、数值稳定，而且不需要给定被积函数，因此能有效解决建模困难的系统或未知系统的求积分问题，在工程实际中有较大的应用价值。下面详细讨论余弦基函数神经网络算法及其与数值积分的关系，并给出计算实例。4.6.1余弦基函数神经网络模型余弦基函数神经网络模型如图4.6-1所示，其中为神经网络权值，为隐层神经元激励函数，且。设权值矩阵为：，激励矩阵为：，则有：神经网络输出： (4.6-1)误差函数： (4.6-2)其中为样本点数，为被积函数。设误差矩阵为：，则有：性能指标： (4.6-3)其中为Euclidean 范数的平方。权值调整： (4.6-4) (4.6-5)其中为学习率，且。【定理4. 5】设为学习率，则当时，神经网络算法是收敛的，其中是隐层神经元个数。【证明】：为了节省篇幅，此略。【定理4. 6】设a、b为积分上下限，且，为神经网络权值，则有：【证明】： ，证毕。由定理4.6可以得出以下三条推论：【推论1】当时，【推论2】当时，【推论3】当时，4.6.2数值积分实例【例4.7】用梯形法和辛甫生方法在0，2积分区间分别计算被积函数：等六个函数的积分，结果如表4.7所示。在神经网络算法中，取神经网络结构为：，性能指标：，学习率为：，训练样本集为：。表4.74精确值2.6676.4001.0992.9581.4166.389梯形法结果4.00016.0001.3333.3260.9098.389Simpson方法结果2.6676.6671.1112.9641.4256.421神经网络算法结果2.6676.3971.0992.9581.4166.388例4.7的源程序(ex4p7.m)如下：clear,format longtol=input(tol=?);N=input(N=?);x=linspace(0,2,N);w=rand(N,1);f=x.2; % 在0，2区间计算被积函数的积分for n=0:N-1 c(n+1,:)=cos(n.*x);endnorm2=norm(c,2).2;eta=1.35./norm2;for n=1:100 y=c*w; e=f-y; J=0.5*norm(e,2).2; w=w+eta*c*e; if J=tol n=,n,break endendp=1:N-1; p=p;I=w(1)*2 % 输出积分值【例4.8】计算积分：用Romberg方法计算该积分时遇到了困难，用Composite Simpsons rule计算时，将积分区间0，48等分成100个子区间，计算结果为58.47082。考虑到被积函数是以为周期的函数，且48=15+0.8761，因此，神经网络计算方法如下：=15+根据定理4.6中的推论1和推论2有：取神经网络结构为：，性能指标：，学习率为：，训练样本集为：，计算结果为：。图4.6-1 余弦基函数神经网络模型例4.8的源程序（ex4p8.m）如下：clear, format longtol=input(tol=?);N=input(N=?);x=linspace(0,pi,N);w=rand(N,1);f=sqrt(1+cos(x).2);for n=0:N-1 c(n+1,:)=cos(n.*x);endeta=1.2./30;for n=1:100 y=c*w; e=f-y; J=0.5*norm(e,2).2; w=w+eta*c*e; if J=tol n=,n,break endendp=1:N-1;p=p;I=15*pi*w(1)+0.8761*w(1)+sum(w(2:N)./p.*sin(0.8761*p)第5章 求解非线性方程的现代方法5.7 求解非线性方程的神经网络算法求解非线性方程一直是国内外广大学者关心的问题之一。前面介绍了求解非线性方程的区间分半法、不动点迭代法、Newton迭代法和弦截法（割线法）等，其中Newton迭代法是解非线性方程的最著名和最有效的数值方法之一，但是，Newton迭代法要求只有当初始值充分接近方程的根时，才能收敛较快，否则难以保证算法的收敛性，为此提出了牛顿下山发以修正迭代过程中出现的严重偏差；区间分半法（也称二分法）是一种简单实用的方法，但收敛速度慢，在实际应用中，通常用该法来确定根的初始近似值；弦截法与Newton迭代法类似，它用差商代替Newton法的导数，在导数计算比较费事或不可能的情况下，显示了它的优点，但一般来说，其收敛速度比Newton法慢；而不动点迭代法（或称Picard法）需要根据给定的非线性方程来构造合适的迭代函数，以实现不动点的迭代计算：。众所周知，迭代函数的构造形式与收敛速度有很大关系，合适的迭代函数可能会使收敛速度很快，而不合适的迭代函数会使收敛速度很慢，甚至根本不收敛，因此，迭代函数的构造十分重要，但是，在实际应用中，如何构造合适的或最佳的迭代函数却不易掌握。此外还有许多其它的改进算法，如半隐式迭代法3等，半隐式迭代法要求在上非零，且具有任意阶连续导数，因此条件较苛刻。基于以上各种算法存在的局限性，编者提出了一种不依赖于初始值与迭代函数选择的新方法，即求解一元非线性方程的神经网络算法。该算法可在指定根的存在区间随机选择初始值，经过神经网络算法训练来获取精度可控的数值解。下面详细介绍求解一元非线性方程的神经网络算法。5.7.1 求解一元非线性方程的神经网络算法（1） 解一元非线性方程的神经网络算法描述设某一元非线性方程为：，且，以为训练权值，训练样本集为：，为隐层神经元激励函数，则神经网络模型如图5.7-1所示。神经网络输出为： (5.7-1)误差函数为： (5.7-2)定义性能指标： (5.7-3)权值调整量（采用最速下降法）：权值调整： (5.7-4)其中为学习率，且。图5.7-1 神经网络模型5.7.2 神经网络算法收敛性研究学习率的大小与算法的收敛性有十分重要的关系，太小了会使收敛速度很慢，增加了计算量；太大了又会使算法出现震荡现象而无法收敛，因此，为了保证算法快速收敛，必须要合理选择学习率的大小。下面给出了学习率的取值定理，为学习率的合理选取提供了理论依据。【定理5.4】设为学习率，在区间上连续，且存在一阶导数，则当学习率为：时，算法是收敛的。【证明】：取Lyapunov函数为：，则有。因为，且，于是有： (5.7-5)由式（5.7-1）知，所以，由式（5.7-5）知，要使神经网络算法收敛，必须有下式成立：因，所以有 ，如果的微分存在且有界，则有： ，所以为了使神经网络算法绝对收敛，取学习率为：，证毕。根据经验，通常取学习率为：时，算法收敛速度最快。5.7.3 神经网络算法步骤step1. 在区间上随机产生初始权值，给定无穷小正数Tol，令J=0，计算，确定学习率；step2. 计算神经网络输出： 计算误差函数：计算性能指标：step3. 权值调整：step4. 如果，返回step2重复上述训练过程；否则结束训练，输出神经网络权值即非线性方程的解：。5.7.4 算例为了检验神经网络算法的有效性，下面给出一个算例，并与其它算法结果作了比较分析。【例5.10】求解非线性方程：在0，1区间的解。【解】：因为，所以有：取学习率为，计算结果如表5.4所示。表5.4 例5.10 计算结果 神经网络算法在0，1区间随机产生初始值迭代次数计算结果误差30.35822156.865436360.35823221.6299892130.35823227.2164497初始值固定为：0.420.35823126.800036630.35823214.5268018120.35823222.7755576半隐式迭代法3未给出0.35167150.3584288-0.027863757723978.386282305697934Newton迭代法4不稳定例5.10的源程序(ex5p10.m)如下：Clear; format longeta=0.05;tol=input(tol=?);xk=rand;for n=1:20 fxk=exp(xk)-1-cos(pi*xk); fdot= exp(xk)+pi*sin(pi*xk); ek=-fxk; J=0.5*ek.2; xk_1=xk; xk=xk+eta*ek*fdot; if J=tol xk=,xk, break endend文献3未给出迭代次数，无法了解算法的收敛速度，而且计算结果依赖初始给定值，不同的初始给定值，其计算结果明显不同，因而算法的稳定性欠佳。由例5.10的计算结果可知，与Newton法和半隐式迭代法相比，神经网络算法具有收敛速度快、计算精度高、计算精度可控、收敛性不依赖初始值的选择等特点，是非线性方程求解的一种有效方法。5.7.5 算法改进由式(5.7-4)可知，神经网络权值训练过程中涉及信息，如果不易求得时，可用差商：代替微分，因此，式(5.7-4)改写为： (5.7-6)当结束训练。此外，神经网络学习率改写为： (5.7-7)【例5.11】用改进神经网络算法求解非线性方程：在0，1区间的解。例5.11的源程序(ex5p11.m)如下：clear; format longtol=input(tol=?);xk=rand;xk_1=rand;for n=1:20 fxk=exp(xk)-1-cos(pi*xk); fxk_1=exp(xk_1)-1-cos(pi*xk_1); ek=-fxk; if xk=xk_1 xk=,xk, break else fdot=(fxk-fxk_1)./(xk-xk_1); %eta=1.2./fdot.2; eta=0.05; end J=0.5*ek.2; xk_1=xk; xk=xk+eta*ek*fdot; if J=tol xk=,xk, break endend (1) 如果使用式(5.7-7)确定学习率，经过15次神经网络训练，结果为:, 误差为(2) 如果学习率取固定值，即，则经过11次神经网络训练后，结果为:，误差为;经过17次神经网络训练后，结果为: ，误差为：。【例5.12】同例5.9，计算下列方程【见式(5.1-1)】的根。【问题背景】：1000是指当年某城市的初始人口1000万，435是指当年迁入本城市的人口435万，指该城市人口出生率，1564指该城市当年末的人口总数1564万。【解】：因为，所以方程在0，1内有解。由于较复杂，在这里近似用差商来代替，则学习率为，式(5.7-7)的权值调整公式改写为：其中。例5.12的源程序(ex5p12.m)如下：(1) 函数定义（population2.m）function y=population2(x)n=length(x);for j=1:n if x(j)=0 y(j)=-129; else y(j)=1000*exp(x(j)+435*(exp(x(j)-1)./x(j)-1564; endend(2) 求方程的解的源程序(ex5p12.m)clear;format longx0=input(x0=?); x1=input(x1=?);N=input(N=?); Tol=input(Tol=?);n=0; flag=0;y0=population2(x0); y1=population2(x1);while n=N L=(y1-y0)./(x1-x0);eta=1.0./L.2;ek=-y1;x=x1+eta*ek*L; if abs(y1)=Tol flag=1;break end n=n+1; x0=x1; x1=x; y0=y1; y1=population2(x1);endif flag=0 The method failed after N iterations N=,Nelse x=,x,n=,n,y1=,y1end取初值，计算结果如表5.5所示。表5.5 例5.12计算结果神经网络算法初始值,迭代次数 0.2, 0.550.100997929685755.58884494239e-1060.100997929685754.54747350886e-13110.100997929685750.00000000000(exact value)0.3, 0.5100.100997929685750.00000000000(exact value)0.3, 0.490.100997929685750.00000000000(exact value)0.4, 0.560.100997929685750.00000000000(exact value)第6章 求解线性方程组的现代方法6.7求解线性方程组的共轭斜量(梯度)法设线性方程组则有，给定初始近似，取 (6.7-1)对计算： (6.7-2) (6.7-3) (6.7-4)且 (6.7-5) (6.7-6)【几点说明】：(1） 应用共轭梯度法求解实对称正定线性方程组，不必事先估计方程组的系数矩阵的特征值的上、下界，不需要选取任何迭代参数，使用方便。而且，它的收敛速度较快，对于非坏条件问题，一般来说，所需的迭代次数远小于矩阵的阶数（如果计算过程中没有舍入误差，理论上，最多迭代步便可得到方程组的准确解。如果矩阵只有个相异的特征值，则共轭梯度法至多进行步迭代，便可得到方程组的准确解）。因此，共轭梯度法目前已得到较为普遍应用。特别对于解大型稀疏矩阵的线性方程组，它是一个有效的方法。(2） 共轭梯度法适用于为对称正定矩阵的情形。对于为非对称矩阵的情形，可将线性方程组化为，此时，设为非奇异，则为对称正定矩阵。(3） 当矩阵的条件数较大时，各种计算方法的收敛速度都会很慢。因为求解线性方程组的收敛速度通常与矩阵的条件数有关，即：，其中是近似解，是精确解，是剩余向量。特别是当矩阵的条件数很大时（很怀条件），却会发生这样的情况：相当精确的解，其剩余向量并不小；不精确的解，反而其剩余向量却很小。由此可见，对于很坏条件线性方程组而言，从剩余向量的大小不能说明相应的近似解精确与否。为此，国外学者提出了改善矩阵条件数的方法。6.8基于梯度下降法的神经网络算法在自然科学和工程应用领域，许多数学问题的求解最终转化为线性方程组的求解问题，如样条函数插值、解非线性方程组、求解偏微分方程的差分方法及有限元方法、Wiener滤波器的优化设计、IIR数字滤波器和FIR数字滤波器的优化设计、经济学中的投入产出问题等等，所有这些问题均涉及到线性方程组的求解，因此，线性方程组的求解问题是一个非常实用又特别重要的问题。许多年来，国内外众多学者为此做了大量研究，取得了许多优秀成果，其中最典型的是众所周知的克莱姆法则和高斯消元法则等直接方法。前者适用于线性方程组有唯一解的情况，即系数矩阵A可逆，但计算量十分巨大，不适合大线性方程组的求解；后者的计算量也很大，特别是当线性方程组的阶数较高或系数矩阵是大稀疏的情况时，高斯消元法不适应。另一类是迭代法和最小二乘法，比较常用的有Jacobi迭代、Gauss-Seidel迭代、超松弛迭代等1-3和共轭梯度法1,4-6。用迭代法求解线性方程组，比较适合大稀疏的情况，但一般来说对线性方程组的系数矩阵有些条件限制，比如对于对称正定的系数矩阵A，Gauss-Seidel迭代法收敛但Jacobi迭代法却不收敛；同样情况，Jacobi迭代法收敛时，Gauss-Seidel迭代法不一定收敛，因此，这两种迭代法均受不同程度的限制。一般而言，Jacobi迭代法和Gauss-Seidel迭代法收敛太慢，在工程实际中很少采用；超松弛迭代法（SOR）是Gauss-Seidel迭代法的一种改进方法，只要合理选择松弛因子，其收敛速度明显加快，但并不是对任意线性方程组，SOR迭代法都收敛。共轭梯度法1-4是目前较为普遍使用的方法，特别对于解大型稀疏矩阵的线性方程组，它是一种有效方法，对于非坏条件问题，其收敛速度很快，但对于矩阵A的条件数较大时，其收敛速度很慢。为了解决各种算法的局限性，编者提出并详细研究了求解线性方程组的神经网络算法。6.8.1 基于梯度下降法(Gradient-descent method)的神经网络算法(NN-GDM)(1) 求解线性方程组的神经网络算法模型设线性方程组为： (6.8-1)其中：，式(6.8-1)可改写为：，() (6.8-2)其中。设 (6.8-3)以为神经网络输入，以为神经网络权值向量，以为神经网络输出，以为神经网络训练样本，则神经网络模型如图6.8-1所示。图6.8-1 神经网络模型由图6.8-1可得神经网络输出为： (6.8-4)误差为： (6.8-5)定义性能指标为： (6.8-6)则权值调整量为（采用最速下降法）： (6.8-7)因此，权值调整公式为：，() (6.8-8)或写为矩阵形式为： (6.8-9)其中为学习率，且，表示系数矩阵的第行元素构成的列向量，即：。(2) 求解线性方程组的神经网络算法收敛性研究众所周知，学习率的大小与算法收敛速度以及是否收敛有很大关系，取大了，算法因振荡而无法收敛；取小了，尽管能保证算法收敛，但收敛速度很慢。为了合理选择学习率的大小，下面给出并证明了学习率的取值定理，为学习率的合理选择提供了理论依据。【定理6.8】设为学习率，当且仅当学习率满足： (6.8-10)时，神经网络算法是收敛的，其中为系数矩阵第行向量的欧氏范数的平方。【证明】：取Lyapunov函数为： (6.8-11)则有 (6.8-12)因为 (6.8-13)且 (6.8-14)于是有： (6.8-15)由式(6.8-3)知，所以， (6.8-16)由式(6.8-15)知，要使神经网络算法收敛，必须有下式成立：因，所以有 。为了是神经网络算法绝对收敛，取学习率为：，证毕。根据经验，通常取学习率为时，算法收敛速度最快。(3) 求解线性方程组的神经网络算法步骤Step1. 随机产生权值列向量，给定任意小正数，令，准备好系数矩阵和向量，确定学习率；Step2. 计算神经网络输出：计算误差函数: 计算性能指标： ；Step3. 权值调整：Step4. 如果样本集未训练完毕，返回Step2 重复上述步骤，否则，判断性能指标是否满足给定任意小正数？如果，令，返回Step2 重复上述步骤，否则结束训练，输出神经网络训练结果，此即为满足给定指标要求的线性方程组的最佳解。6.8.2应用实例为了检验基于梯度下降法的神经网络计算方法(NN-GDM)求解线性方程组的有效性，下面给出了几个算例：【例6.10】文献1（第6、9页）中例1和例2，方程组有精确解为=1，1，1T。分别用Jacobi迭代法、Gauss-Seidel迭代法和本文算法求解线性方程组，计算结果如表6.1所示。 由表6.1计算结果可以看出，与Jacobi迭代法和Gauss-Seidel法相比，NN-GDM计算方法不仅具有更高的计算精度，而且还能得到精确解，这是Jacobi迭代法和Gauss-Seidel迭代法无法比拟的。表6.1 【例6.10】计算结果计算方法迭代次数计算结果最大绝对误差：Jacobi迭代法60.9997570.9992710.999028Gauss-Seidel法40.9998299000.9998469100.999938764NN-GDM初始值为0,0,0T51.000226080119961.000143900253420.9999852521876970.999999833028010.999997893861730.99999967173180201.000000000000001.000000000000001.000000000000000.00（精确解）例6.10源程序(ex6p10.m)如下：clear, format longA=10,-1,0;-1,10,-2;0,-4,10;b=9,7,6;n=length(b);for j=1:n normA(j)=norm(A(j,:),2);end max_norm=max(normA); eta=1.2./max_norm.2; x=zeros(n,1); J=0; tol=input(tol=?); for m=1:30 for k=1:n yk=A(k,:)*x; ek=b(k)-yk; J=J+0.5*ek.2; x=x+eta*ek*A(k,:); end if J=tol x=,x,m=,m,break end J=0; endxd=1,1,1;maxaer=norm(xd-x,inf)【例6.11】文献2（第200201页）中的例5.10，方程组有唯一解= 1，1 T 。由文献2知Jacobi迭代法和Gauss-Seidel迭代法均不收敛。用本文算法求解线性方程组，计算结果如表6.2所示。 表6.2 例6.11计算结果计算方法迭代次数计算结果最大绝对误差：Jacobi迭代法1发散Gauss-Seidel法1发散NN-GDM初始值为0,0T100.999999799922560.99999923449896190.999999999999830.99999999999933261.000000000000001.000000000000000.00（精确解）例6.11源程序(ex6p11.m)如下：clear, format longA=3,-10;9,-4;b=-7,5;n=length(b);for j=1:n normA(j)=norm(A(j,:),2);end max_norm=max(normA); eta=1.2./max_norm.2; x=zeros(n,1); J=0; tol=input(tol=?); for m=1:30 for k=1:n yk=A(k,:)*x; ek=b(k)-yk; J=J+0.5*ek.2; x=x+eta*ek*A(k,:); end if J=tol x=,x,m=,m,break end J=0; endxd=1,1;maxaer=norm(xd-x,inf) 【例6.12】文献1（第41页）中的例1，方程组有唯一解= 1，1，1 T。分别用共轭斜量法和神经网络算法进行计算，计算结果如表6.3所示。 表6.3 例6.12计算结果计算方法迭代次数计算结果最大绝对误差：共轭斜量法初始值为0,0,0T60.999999999985780.999999999985780.99999999998578101.000000000000001.000000000000001.000000000000000.00（精确解）NN-GDM初始值为0,0,0T361.000000000000100.999990433698200.99999999999998831.000000000000000.999999999897881.000000000000002001.000000000000001.000000000000001.00000000000000(1) 采用共轭斜量法的源程序(6p12_1.m)如下：clear, format longA=2,0,1;0,1,0;1,0,2;b=3,1,3;n=length(b);x=zeros(n,1);xd=1,1,1;tol=input(tol=?);for m=1:20 x=cgm(A,b,x); maer=norm(xd-x,inf); if maer=tol x=,x,m=,m,maer=,maer,break endend其中自定义的共轭斜量法函数cgm( )为：function x=cgm(A,b,x)m=length(b);for j=1:m pk=b-A*x; rk=-pk; papk=pk*A*pk; if papk=0 alfak=-rk*pk./papk; x=x+alfak*pk; rk=rk+alfak*A*pk; betak=rk*A*pk./papk; pk=-rk+betak*pk; else break endend(2) 采用神经网络算法的源程序(6p12_2.m)如下：clear, format longA=2,0,1;0,1,0;1,0,2;b=3,1,3;n=length(b);for j=1:n normA(j)=norm(A(j,:),2);end max_norm=max(normA); eta=1.2./max_norm.2; x=rand(n,1); J=0; tol=input(tol=?); for m=1:200 for k=1:n yk=A(k,:)*x; ek=b(k)-yk; J=J+0.5*ek.2; x=x+eta*ek*A(k,:); end if J=tol x=,x,m=,m,break end J=0; endxd=1,1,1;maxaer=norm(xd-x,inf)【例6.13】已知某线性系统方程有精确解：，。由于条件数，显然是一个病态矩阵方程，各类直接方法和迭代法都难以得到高精度解。同样，采用NN-GDM方法求解该方程组也无法得到高精度解。【例6.14】已知某线性系统方程有精确解：，。由于条件数，显然是一个病态矩阵方程，各类直接方法和迭代法以及NN-GDM方法都难以得到高精度解。【讨论】由例6.10和例6.11的计算结果可知，与Gauss-Seidel算法和Jacobi方法相比，基于梯度下降法的神经网络算法(NN-GDM)具有更高的计算精度。特别是针对例6.11的问题，Gauss-Seidel算法和Jacobi方法均发散，而NN-GDM方法却能得到精确解，显示了该算法的优越性。但是，NN-GDM也存在一些问题，比如例6.12的系数矩阵为对称矩阵时，共轭斜量法能得到精确解，而NN-GDM方法却难以得到精确解，表明了该方法的不足。例6.13和例6.14都是病态线性方程组，针对病态线性方程组的求解问题，NN-GDM方法和其它方法一样都无法得到高精度解。为了有效解决病态线性方程组的求解问题，作者提出了基于递推最小二乘算法的神经网络计算方法(简称NN-RLS)。6.9 基于递推最小二乘算法的神经网络计算方法(NN-RLS) 将式(6.8-3)重写如下： (6.9-1)神经网络模型如图6.8-1所示。采用递推最小二乘法(RLS)调整神经网络权值，算法描述如下：(1) 神经网络输出： (6.9-2)(2) 误差函数： (6.9-3)(3) 性能指标： (6.9-4)(4) 使的权值向量，即为所求的神经网络权值向量，这是一个多变量线性优化问题，为此，由可得递推最小二乘法（RLS）： (6.9-5) (6.9-6) (6.9-7)随机产生初始权值向量，设（是足够大的正数（一般取），是单位矩阵），通过对样本数据训练，即可获得神经网络权值向量，此即为线性方程组(6.8-1)的解。【几点说明】：(1) 当时，式(6.9-5)式(6.9-7)为基本递推最小二乘算法；(2) 的选择范围一般在：。当参数变化快时，取小点；变化慢时，取大点。为了验证基于递推最需二乘法的神经网络计算方法的有效性，下面给出几个实例。【例6.15】同例6.10，已知方程组有精确解为=1，1，1T。分别用Jacobi迭代法、Gauss-Seidel迭代法、NN-GDM方法和NN-RLS方法求解线性方程组，计算结果如表6.4所示。在NN-RLS方法中，,。 表6.4 例6.15计算结果计算方法迭代次数计算结果最大绝对误差：Jacobi迭代法160.9997570.9992710.999028Gauss-Seidel法140.9998299000.9998469100.999938764NN-GDM初始值：0,0,0T51.000226080119961.000143900253420.9999852521876970.999999833028010.999997893861730.99999967173180201.000000000000001.000000000000001.000000000000000.00（精确解）NN-RLS初始值：0,0,0T,10.999999986120030.999999977683860.9999999777768431.000000000000001.000000000000001.000000000000000.00（精确解）例6.15的源程序(ex6p15.m)如下：clear, format longA=10,-1,0;-1,10,-2;0,-4,10;b=9,7,6;n=length(b);x=zeros(n,1);xd=1,1,1;tol=input(tol=?);alfa=1e+6;nda=1;for k=1:10 x=rls(A,b,x,alfa,nda); maxaer=norm(xd-x,inf); %无穷范数 if maxaer=tol x=,x,maxaer=,maxaer,k=,k,break endend其中，自定义函数rls( )如下：function x=rls(A,b,x,alfa,nda)m=length(b);I=eye(m);pk=alfa*I;for k=1:m ek=b(k)-A(k,:)*x; qk=pk*A(k,:)./(nda+A(k,:)*pk*A(k,:); x=x+qk*ek; pk=(pk-qk*A(k,:)*pk)./nda;end【例6.16】同例6.11，已知方程组有唯一解= 1，1 T 。由文献2知Jacobi迭代法和Gauss-Seidel迭代法均不收敛。用NN-GDM 方法和NN-RLS算法求解线性方程组，计算结果如表6.5所示。在NN-RLS算法中，,。源程序与例6p15.m类似。 表6.5 例6.16计算结果计算方法迭代次数计算结果最大绝对误差：Jacobi迭代法1发散Gauss-Seidel法1发散NN-GDM初始值为0,0T100.999999799922560.99999923449896190.999999999999830.99999999999933261.000000000000001.000000000000000.00（精确解）NN-RLS初始值：0,0T,21.000000000000001.0000000000000031.000000000000001.000000000000000.00（精确解）【例6.17】同例6.12，已知方程组有唯一解= 1，1，1 T。分别用共轭斜量法和神经网络算法进行计算，计算结果如表6.6所示。 表6.6 例6.17计算结果计算方法迭代次数计算结果最大绝对误差：共轭斜量法初始值为0,0,0T60.999999999985780.999999999985780.99999999998578101.000000000000001.000000000000001.000000000000000.00（精确解）NN-GDM初始值为0,0,0T361.000000000000100.999990433698200.99999999999998831.000000000000000.999999999897881.000000000000002001.000000000000001.000000000000001.00000000000000NN-RLS初始值：0,0,0T,20.999999999999990.999999999999000.99999999999999201.000000000000001.000000000000001.00000000000000NN-RLS初始值：0,0,0T,10.999999999988500.999999999907840.9999999999910131.000000000000001.000000000000001.000000000000000.00（精确解）【例6.18】同例6.13，已知某线性系统方程有精确解：，。由于条件数，显然是一个病态矩阵方程，各类直接方法和迭代法都难以得到高精度解。同样，采用NN-GDM方法求解该方程组也无法得到高精度解。采用NN-RLS计算方法，计算结果如表6.7所示。表6.7 例6.18计算结果计算方法迭代次数计算结果最大绝对误差：NN-RLS初始值：0,0T,10.999999159105481.0000004204690521.000000000000120.99999999999994101.000000000000120.99999999999994【例6.19】同例6.14，已知某线性系统方程有精确解：，。由于条件数，显然是一个病态矩阵方程，各类直接方法和迭代法以及NN-GDM方法都难以得到高精度解。采用NN-RLS计算方法，计算结果如表6.8所示。表6.7 例6.19计算结果计算方法迭代次数计算结果最大绝对误差：NN-RLS初始值：0,0T,20.99999999393620-0.9999999915989930.99999999996747-0.99999999995493100.99999999993622-0.99999999991163第7章微分方程初值问题求解7.10 求解常微分方程初值问题的神经网络算法在许多科学研究和工程技术中所遇到的微分方程，只有少数是较简单和典型的常微分方程，比如线性常系数微分方程能够用初等方法求得它们的解析解。对于变系数微分方程的求解就有困难，更不用说一般的非线性微分方程了。因此，在多数情况下，常微分方程只能用近似方法求解。近似方法有两类：一类称为近似解析方法，如级数解法、逐次逼近法等；另一类近似方法称为数值解法，它可以给出解在一些离散点上的近似值。怎样利用计算机实现常微分方程初值问题的近似解就是本章将要分析讨论的数值计算方法问题。由于高阶常微分方程可以转化为一阶常微分方程组，因此，为了不失一般性，本章主要介绍一类一阶常微分方程初值问题 (7.10-1)在区间上的数值计算方法。假定函数足够光滑，且满足李普希兹（lipshitz）条件，以保证初值问题（7.10-1）的解存在且唯一。所谓初