C++实现 牛顿迭代 解非线性方程组
C++实现牛顿迭代解非线性方程组（二元二次为例）
求解{0=x*x-2*x-y+0.5;
0=x*x+4*y*y-4;
}的方程
#include
#include
#define N 2 // 非线性方程组中方程个数、未知量个数#define Epsilon 0.0001 // 差向量1范数的上限
#define Max 100 // 最大迭代次数
using namespace std;
const int N2=2*N;
int main()
{
void ff(float xx[N],float yy[N]); //计算向量函数的因变量向量yy[N]
void ffjacobian(float xx[N],float yy[N][N]); //计算雅克比矩阵yy[N][N]
void inv_jacobian(float yy[N][N],float inv[N][N]); //计算雅克比矩阵的逆矩阵inv
void newdundiedai(float x0[N], float inv[N][N],float y0[N],float x1[N]);
//由近似解向量x0 计算近似解向量x1
float x0[N]={2.0,0.25},y0[N],jacobian[N][N],invjacobian[N][N],x1[N],errornorm;
int i,j,iter=0;
//如果取消对x0的初始化，撤销下面两行的注释符，就可以由键盘x读入初始近似解向量for( i=0;i cin>>x0[i]; cout<<"初始近似解向量："< for (i=0;i cout< cout< do { iter=iter+1; cout<<"第"< ff(x0,y0); //计算雅克比矩阵jacobian ffjacobian(x0,jacobian); //计算雅克比矩阵的逆矩阵invjacobian inv_jacobian(jacobian,invjacobian); //由近似解向量x0 计算近似解向量x1 newdundiedai(x0, invjacobian,y0,x1); //计算差向量的1范数errornorm errornorm=0; for (i=0;i errornorm=errornorm+fabs(x1[i]-x0[i]); if (errornorm for (i=0;i x0[i]=x1[i]; } while (iter return 0; } void ff(float xx[N],float yy[N]) //调用函数 {float x,y; int i; x=xx[0]; y=xx[1]; yy[0]=x*x-2*x-y+0.5; yy[1]=x*x+4*y*y-4; //计算初值位置的值 cout<<"向量函数的因变量向量是："< for( i=0;i cout< cout< cout< } void ffjacobian(float xx[N],float yy[N][N]) { float x,y; int i,j; x=xx[0]; y=xx[1]; //jacobian have n*n element //计算函数雅克比的值yy[0][0]=2*x-2; yy[0][1]=-1; yy[1][0]=2*x; yy[1][1]=8*y; cout<<"雅克比矩阵是："< for( i=0;i {for(j=0;j cout< cout< } cout< } void inv_jacobian(float yy[N][N],float inv[N][N]) {float aug[N][N2],L; int i,j,k; cout<<"开始计算雅克比矩阵的逆矩阵："< { for(j=0;j aug[i][j]=yy[i][j]; for(j=N;j if(j==i+N) aug[i][j]=1; else aug[i][j]=0; } for (i=0;i { for(j=0;j cout< cout< } cout< for (i=0;i { for (k=i+1;k {L=-aug[k][i]/aug[i][i]; for(j=i;j aug[k][j]=aug[k][j]+L*aug[i][j]; } } for (i=0;i { for(j=0;j cout< cout< cout< for (i=N-1;i>0;i--) { for (k=i-1;k>=0;k--) {L=-aug[k][i]/aug[i][i]; for(j=N2-1;j>=0;j--) aug[k][j]=aug[k][j]+L*aug[i][j]; } } for (i=0;i { for(j=0;j cout< cout< } cout< for (i=N-1;i>=0;i--) for(j=N2-1;j>=0;j--) aug[i][j]=aug[i][j]/aug[i][i]; for (i=0;i { for(j=0;j cout< cout< for(j=N;j inv[i][j-N]=aug[i][j]; } cout< cout<<"雅克比矩阵的逆矩阵："< { for(j=0;j cout< cout< cout< } void newdundiedai(float x0[N], float inv[N][N],float y0[N],float x1[N]) { int i,j; float sum=0; for(i=0;i { sum=0; for(j=0;j sum=sum+inv[i][j]*y0[j]; x1[i]=x0[i]-sum; } cout<<"近似解向量："< for (i=0;i cout< cout< } 解线性方程组的迭代法 1.rs里查森迭代法求线性方程组Ax=b的解 function[x,n]=rs(A,b,x0,eps,M) if(nargin==3) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值elseif(nargin==4) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1; %迭代过程 while(tol>eps) x=(I-A)*x0+b; n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x; if(n>=M) disp('Warning:迭代次数太多，可能不收敛！'); return; end end 2.crs里查森参数迭代法求线性方程组Ax=b的解 function[x,n]=crs(A,b,x0,w,eps,M) if(nargin==4) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值 elseif(nargin==5) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1; %迭代过程 while(tol>eps) x=(I-w*A)*x0+w*b; n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x; if(n>=M) disp('Warning:迭代次数太多，可能不收敛！'); return; end end 3.grs里查森迭代法求线性方程组Ax=b的解 function[x,n]=grs(A,b,x0,W,eps,M) if(nargin==4) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值 elseif(nargin==5) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1;%前后两次迭代结果误差 %迭代过程 while(tol>eps) x=(I-W*A)*x0+W*b;%迭代公式 n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x; if(n>=M) disp('Warning:迭代次数太多，可能不收敛！'); return; end end 4.jacobi雅可比迭代法求线性方程组Ax=b的解 function[x,n]=jacobi(A,b,x0,eps,varargin) if nargin==3 eps=1.0e-6; M=200; elseif nargin<3 error return elseif nargin==5 M=varargin{1}; end D=diag(diag(A));%求A的对角矩阵 L=-tril(A,-1);%求A的下三角阵 学科前沿讲座论文 班级：工程力学13-1班姓名：陆树飞 学号：02130827 牛顿法求非线性方程的根 一 实验目的 （1）用牛顿迭代法求解方程的根 （2）了解迭代法的原理，了解迭代速度跟什么有关 题目：用Newton 法计算下列方程 （1） 013=--x x ， 初值分别为10=x ，7.00=x ，5.00=x ； （2） 32943892940x x x +-+= 其三个根分别为1,3,98-。当选择初值02x =时 给出结果并分析现象,当6510ε-=?，迭代停止。 二 数学原理 对于方程f(x)=0，如果f(x)是线性函数，则它的求根是很容易的。牛顿迭代法实质上是一种线性化方法，其基本思想是将非线性方程f(x)=0逐步归结为某种线性方程来求解。 设已知方程f(x)=0有近似根x k (假定k f'(x )0≠) ，将函数f(x)在点x k 进行泰勒展开，有 k k k f(x)f(x )+f'(x )(x-x )+≈??? 于是方程f(x)=0可近似的表示为 k k k f(x )+f'(x )(x-x )=0 这是个线性方程，记其根为x k+1，则x k+1的计算公式为 k+1k ()x =x -'() k k f x f x ，k=0,1,2,… 这就是牛顿迭代法。 三 程序设计 （1）对于310x x --=，按照上述数学原理，编制的程序如下 program newton implicit none real :: x(0:50),fx(0:50),f1x(0:50)!分别为自变量x ，函数f(x)和一阶导数f1(x) integer :: k write(*,*) "x(0)=" read(*,*) x(0) !输入变量：初始值x(0) open(10,file='1.txt') do k=1,50,1 fx(k)=x(k-1)**3-x(k-1)-1 f1x(k)=3*x(k-1)**2-1 x(k)=x(k-1)-fx(k)/f1x(k) !牛顿法 write(*,'(I3,1x,f11.6)') k,x(k) !输出变量：迭代次数k 及x 的值 write(10,'(I3,1x,f11.6)') k,x(k) if(abs(x(k)-x(k-1))<1e-6) exit !终止迭代条件 end do stop end （2）对于32943892940x x x +-+=，按照上述数学原理，编制的程序如下 program newton implicit none C++实现牛顿迭代解非线性方程组（二元二次为例） 求解{0=x*x-2*x-y+0.5; 0=x*x+4*y*y-4; }的方程 #include 《MATLAB 程序设计实践》课程考核 ---第37-38页 题１ ： 编程实现以下科学计算算法，并举一例应用之。（参考书籍《精 通ＭＡT ＬＡＢ科学计算》，王正林等著，电子工业出版社，２００９ 年） “牛顿法非线性方程求解” 弦截法本质是一种割线法，它从两端向中间逐渐逼近方程的根；牛顿法本质上是一种切线法，它从一端向一个方向逼近方程的根，其递推公式为： - =+n n x x 1) ()(' n n x f x f 初始值可以取)('a f 和)('b f 的较大者，这样可以加快收敛速度。 和牛顿法有关的还有简化牛顿法和牛顿下山法。 在MATLAB 中编程实现的牛顿法的函数为：NewtonRoot 。 功能：用牛顿法求函数在某个区间上的一个零点。 调用格式：root=NewtonRoot )(```eps b a f 其中，f 为函数名； a 为区间左端点； b 为区间右端点 eps 为根的精度； root 为求出的函数零点。 ， 牛顿法的matlab程序代码如下： function root=NewtonRoot(f,a,b,eps) %牛顿法求函数f在区间[a,b]上的一个零点%函数名:f %区间左端点:a %区间右端点:b %根的精度:eps %求出的函数零点:root if(nargin==3) eps=1.0e-4; end f1=subs(sym(f),findsym(sym(f)),a); f2=subs(sym(f),findsym(sym(f)),b); if (f1==0) root=a; end if (f2==0) root=b; end if (f1*f2>0) disp('两端点函数值乘积大于0 ！'); return; else tol=1; fun=diff(sym(f)); %求导数 fa=subs(sym(f),findsym(sym(f)),a); fb=subs(sym(f),findsym(sym(f)),b); dfa=subs(sym(fun),findsym(sym(fun)),a); dfb=subs(sym(fun),findsym(sym(fun)),b); if(dfa>dfb) %初始值取两端点导数较大者 root=a-fa/dfa; else root=b-fb/dfb; end while(tol>eps) r1=root; fx=subs(sym(f),findsym(sym(f)),r1); dfx=subs(sym(fun),findsym(sym(fun)),r1); %求该点的导数值 root=r1-fx/dfx; %迭代的核心公式 tol=abs(root-r1); end end 例：求方程3x^2－exp（x）＝0的一根 解：在MATLAB命令窗口输入： >> r=NewtonRoot('3*x^2-exp(x)',3,4) 输出结果： X=3.7331 作业六：分别编写用Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组Ax=B的标准程序，并求下列方程组的解。 可取初始向量 X(0) =(0，0，0)’； 迭代终止条件||x(k+1)-x(k)||<=10e-6 （1） = （2） = Jacobi迭代法： 流程图 开 始 判断b中的最大值 有没有比误差大 给x赋初值 进行迭代 求出x，弱到100次还没到，警告不收 结束 程序 clear;clc; A=[8,-1,1;2,10,01;1,1,-5]; b=[1;4;3]; e=1e-6; x0=[0;0;0]'; n=length(A); x=zeros(n,1); k=0; r=max(abs(b)); while r>e for i=1:n d=A(i,i); if abs(d) 程序结果（1） （2） Gauss-Seidel迭代法： 程序 clear;clc; %A=[8,-1,1;2,10,01;1,1,-5]; %b=[1;4;3]; A=[5,2,1;-1,4,2;2,-3,10]; b=[-12;20;3]; m=size(A); if m(1)~=m(2) error('矩阵A不是方阵'); end n=length(b); %初始化 N=0;%迭代次数 L=zeros(n);%分解A=D+L+U,D是对角阵，L是下三角阵，U是上三角阵U=zeros(n); D=zeros(n); G=zeros(n);%G=-inv(D+L)*U d=zeros(n,1);%d=inv(D+L)*b x=zeros(n,1); for i=1:n%初始化L和U for j=1:n if i 中南大学MATLAB程序设计实践学长有爱奉献，下载填上信息即可上交，没有下载券的自行百度。所需m文件照本文档做即可，即新建（FILE）→脚本（NEW-Sscript）→复制本文档代码→运行（会跳出保存界面，文件名默认不要修改，保存）→结果。第一题需要把数据文本文档和m文件放在一起。全部测试无误，放心使用。本文档针对做牛顿法求非线性函数题目的同学，当然第一题都一样，所有人都可以用。←记得删掉这段话 班级： ？ 学号： 姓名： 一、《MATLAB程序设计实践》Matlab基础 表示多晶体材料织构的三维取向分布函数（f＝f（φ1，φ，φ2））是一个非常复杂的函数，难以精确的用解析函数表达，通常采用离散 空间函数值来表示取向分布函数，是三维取向分布函数的一个实例。 由于数据量非常大，不便于分析，需要借助图形来分析。请你编写一 个matlab程序画出如下的几种图形来分析其取向分布特征： （1）用Slice函数给出其整体分布特征； " ~ （2）用pcolor或contour函数分别给出(φ2＝0, 5, 10, 15, 20, 25, 30, 35 … 90)切面上f分布情况(需要用到subplot函数)； (3) 用plot函数给出沿α取向线(φ1=0~90，φ＝45，φ2＝0)的f分布情况。 ( 备注：数据格式说明 解： （1）（ （2）将文件内的数据按照要求读取到矩阵f(phi1,phi,phi2)中，代码如 下： fid=fopen(''); for i=1:18 tline=fgetl(fid); end phi1=1;phi=1;phi2=1;line=0; f=zeros(19,19,19); [ while ~feof(fid) tline=fgetl(fid); data=str2num(tline); line=line+1;数据说明部分，与 作图无关此方向表示f随着 φ1从0,5,10,15, 20 …到90的变化而 变化 此方向表示f随着φ 从0,5,10,15, 20 … 到90的变化而变化 表示以下数据为φ2＝0的数据，即f（φ1，φ，0） 求解线性方程组——超松弛迭代法 #include cout<<"输入松弛因子w (1 (1)二分法求解非线性方程： #include Gauss-Seidel迭代法求解线性方程组 一． 问题描述 用Gauss-Seidel 迭代法求解线性方程组 由Jacobi 迭代法中，每一次的迭代只用到前一次的迭代值。使用了两倍的存储空间，浪费了存储空间。若每一次迭代充分利用当前最新的迭代值，即在计算第i 个分量 ) 1(+k i x 时，用最新分量 ) 1(1 +k x ， ???+) 1(2 k x ) 1(1 -+k i x 代替旧分量 ) (1 k x ， ???) (2 k x ) (1 -k i x ，可以起 到节省存储空间的作用。这样就得到所谓解方程组的Gauss-Seidel 迭代法。 二． 算法设计 将A 分解成U D L A --=，则b x =A 等价于b x =--U)D (L 则Gauss-Seidel 迭代过程 ) ()1()1(k k k Ux Lx b Dx ++=++ 故 ) ()1()(k k Ux b x L D +=-+ 若设1 )(--L D 存在，则 b L D Ux L D x k k 1)(1)1()()(--+-+-= 令 b L D f U L D G 11)()(---=-=， 则Gauss-Seidel 迭代公式的矩阵形式为 f Gx x k k +=+) () 1( 其迭代格式为 T n x x x x ) ()0()0(2)0(1)0(，，，???= (初始向量)， ) (1 1 1 1 1 )()1()1(∑∑-=-+=++--=i j i i j k j ij k j ij i ii i i x a x a b a x )210i 210(n k ???=???=，，，；，，， 或者 ?? ???--=???=???==?+=∑∑-=-+=+++) (1)210i 210(111 1)()1()1()()1(i j i i j k j ij k j ij i ii i i i k i k i x a x a b a x n k k x x x ，，，；，，， 三． 程序框图 牛顿迭代法求解非线性方程组 非线性方程组如下： 221122121210801080 x x x x x x x ?-++=??+-+=?? 给定初值()00.0T x =，要求求解精度达到0.00001 1.首先建立函数()F X ,方程编程如下，将F.m 保存到工作路径中： function f=F(x) f(1)=x(1)^2-10*x(1)+x(2)^2+8; f(2)=x(1)*x(2)^2+x(1)-10*x(2)+8; f=[f(1),f(2)] ; 2.建立函数()DF X ,用于求方程的jacobi 矩阵,将DF.m 保存到工作路径中： function df=DF(x) df=[2*x(1)-10,2*x(2);x(2)^2+1,2*x(1)*x(2)-10]; %jacobi 矩阵是一阶偏导数以一定方式排列成的矩阵。 3.编程牛顿迭代法解非线性方程组，将newton.m 保存在工作路径中： clear,clc; x=[0,0]'; f=F(x); df=DF(x); fprintf('%d %.7f %.7f\n',0,x(1),x(2)); N=4; for i=1:N y=df\f'; x=x-y; f=F(x); df=DF(x); fprintf('%d %.7f %.7f\n',i,x(1),x(2)); if norm(y)<0.0000001 break; else end end ezplot('x^2-10*x+y^2+8',[-6,6,-6,6]); hold on ezplot('x*y^2+x-10*y+8',[-6,6,-6,6]); 运行结果如下： 0 0.0000000 0.0000000 1 0.8000000 0.8800000 2 0.9917872 0.9917117 迭代法解线性方程组作业 沈欢00986096 北京大学工学院，北京100871 2011年10月12日 摘要 由所给矩阵生成系数矩阵A和右端项b，分析系数矩阵A，并用Jacobi迭代法、GS迭代法、SOR(逐步松弛迭代法)解方程组Ax=b 1生成系数矩阵A、右端项b,并分析矩阵A 由文件”gr900900c rg.mm”得到了以.mm格式描述的系数矩阵A。A矩阵是900?900的大型稀 疏对称矩阵。于是，在matlaB中，使用”A=zeros(900,900)”语句生成900?900的零矩阵。再 按照.mm文件中的描述，分别对第i行、第j列的元素赋对应的值，就生成了系数矩阵A，并 将A存为.mat文件以便之后应用。 由于右端项是全为1的列向量，所以由语句”b=ones(900,1)”生成。 得到了矩阵A后，求其行列式，使用函数”det(A)”,求得结果为”Inf”,证明行列式太大，matlaB无法显示。由此证明，矩阵A可逆，线性方程组 Ax=b 有唯一解。 接着，判断A矩阵是否是对称矩阵(其实，这步是没有必要的，因为A矩阵本身是对称矩阵，是.mm格式中的矩阵按对称阵生成的)。如果A是对称矩阵，那么 A?A T=0 。于是，令B=A?A T,并对B求∞范数。结果显示： B ∞=0,所以，B是零矩阵，也就是：A是对称矩阵。 然后，求A的三个条件数： Cond(A)= A ? A?1 所求结果是，对应于1范数的条件数为：377.2334;对应于2范数的条件数为：194.5739;对应 于3范数的条件数为：377.2334; 1 从以上结果我们看出，A是可逆矩阵，但是A的条件数很大，所以，Ax=b有唯一解并且矩阵A相对不稳定。所以，我们可以用迭代方法来求解该线性方程组，但是由于A的条件数太大迭代次数一般而言会比较多。 2Jacobi迭代法 Jacobi迭代方法的程序流程图如图所示： 图1:Jacobi迭代方法程序流程图 在上述流程中，取x0=[1,1,...,1]T将精度设为accuracy=10?3，需要误差满足： error= x k+1?x k x k+1 基于Matlab 实现牛顿迭代法解非线性方程组 已知非线性方程组如下 2211221212 10801080x x x x x x x ?-++=??+-+=?? 给定初值0(0,0)T x =，要求求解精度达到0.00001 首先建立函数F(x)，方程组编程如下，将F.m 保存到工作路径中： function f=F(x) f(1)=x(1)^2-10*x(1)+x(2)^2+8; f(2)=x(1)*x(2)^2+x(1)-10*x(2)+8; f=[f(1) f(2)]; 建立函数DF(x),用于求方程组的Jacobi 矩阵，将DF.m 保存到工作路径中： function df=DF(x) df=[2*x(1)-10,2*x(2);x(2)^2+1,2*x(1)*x(2)-10]; 编程牛顿迭代法解非线性方程组，将newton.m 保存到工作路径中： clear; clc x=[0,0]'; f=F(x); df=DF(x); fprintf('%d %.7f %.7f\n',0,x(1),x(2)); N=4; for i=1:N y=df\f'; x=x-y; f=F(x); df=DF(x); fprintf('%d %.7f %.7f\n',i,x(1),x(2)); if norm(y)<0.0000001 break ; else end end 运行结果如下： 0 0.0000000 0.0000000 1 0.8000000 0.8800000 2 0.9917872 0.9917117 3 0.9999752 0.9999685 4 1.0000000 1.0000000 数值分析实验报告 0 a 12 K a 1,n 1 K a 2,n 1 U O M 则有： 第一步： Jacobi 迭代法 a 1n a 2n M , 则有： A D L U a n 1,n Ax b A A x D b L U (D L U)x b Dx (L U)x b x D (L U)x D b 令 J D (L U) 则称 J 为雅克比迭代矩阵 f D b 由此可得雅克比迭代的迭代格式如下： x (0) , 初始向量 x (k 1) Jx (k) f ,k 0,1,2,L 第二步 Gauss-Seidel 迭代法 Ax b (D L U )x b (D L)x Ux b x (D L) Ux (D L) b A D L U a 11 a 12 L a 1n a 11 A a 21 a 22 L a 2n a 22 M MM MO a n1 a n2 L a nn a 11 得到 D a 22 O a nn 由 a 21 0 M M O a n 1,1 a n 1,2 L 0 a nn a n1 a n2 L a n,n a 21 L M M O a n 1,1 a n 1,2 L a n1 a n2 L a n,n 1 a 12 K a 1,n 1 a 1n 0 K a 2,n 1 a 2n O M M a n 1,n 10 令 G (D L) U ，则称G 为Gauss-Seidel 迭代矩阵 f (D L) b 由此可得 Gauss-Seidel 迭代的迭代格式如下： x (0) , 初始向量 第三步 SOR 迭代法 w0 AD L U 1 ( D 1 wL ((1 w)D wU )) (D 1 wL) ((1 w)D wU ) w w w 令M w 1 (D wL), N 1 ((1 w)D wU )则有：A MN w w Ax b AM L W N M (M N )x b Mx Nx b x M Nx M b N M, 令W f Mb 带入 N 的值可有 L W ((1 w)D wU) (D wL) 1((1 w)D wU) (D wL) f 1 b w 1(D wL) 1b 1 (D wL) w 称 L W 为 SOR 迭代矩阵，由此可得 SOR 迭代的迭代格式如下： x (0) ,初始向量 二、算法程序 Jacobi 迭代法的 M 文件： function [y,n]=Jacobi(A,b,x0,eps) %************************************************* %函数名称 Jacobi 雅克比迭代函数 %参数解释 A 系数矩阵 % b 常数项 % x0 估计解向量 x (k 1) Gx (k) f ,k 0,1,2,L (k 1) f,k 0,1,2,L 一． 问题描述 用Gauss-Seidel 迭代法求解线性方程组 由Jacobi 迭代法中，每一次的迭代只用到前一次的迭代值。使用了两倍的存储空间，浪 费了存储空间。若每一次迭代充分利用当前最新的迭代值，即在计算第i 个分量) 1(+k i x 时， 用最新分量) 1(1 +k x ，???+) 1(2 k x ) 1(1 -+k i x 代替旧分量)(1k x ，???) (2 k x ) (1-k i x ，可以起到节省存储 空间的作用。这样就得到所谓解方程组的Gauss-Seidel 迭代法。 二． 算法设计 将A 分解成U D L A --=，则b x =A 等价于b x =--U)D (L 则Gauss-Seidel 迭代过程 ) ()1()1(k k k Ux Lx b Dx ++=++ 故 )()1()(k k Ux b x L D +=-+ 若设1 )(--L D 存在，则 b L D Ux L D x k k 1)(1)1()()(--+-+-= 令 b L D f U L D G 11)()(---=-=， 则Gauss-Seidel 迭代公式的矩阵形式为 f Gx x k k +=+)()1( 其迭代格式为 T n x x x x )()0()0(2)0(1)0(，，，???= (初始向量)， )(1111 1 )() 1()1(∑∑-=-+=++--=i j i i j k j ij k j ij i ii i i x a x a b a x )210i 210(n k ???=???=，，，；，，， 或者 ?? ???--=???=???==?+=∑∑-=-+=+++) (1)210i 210(111 1)() 1()1()()1(i j i i j k j ij k j ij i ii i i i k i k i x a x a b a x n k k x x x ，，，；，，， 三． 程序框图 《数值计算方法》实验报告 实验名称：实验1 非线性方程的简单迭代法和Steffensen 迭代法 实验题目：分别用简单迭代法和Steffensen 迭代法求方程 010423=-+x x 在 [1, 2] 内的一个实根. 实验目的：理解并掌握简单迭代法和Steffensen 迭代法 基础理论：简单迭代法和Steffensen 迭代法 1）.简单迭代法的原理：将一元非线性方程：0)(=x f 改写成等价方程：)(x x ρ= ，对此，从某个初始值x0开始，对应式)(x x ρ= 构成迭代公式 ,...1,0),(1==+k x x k k ρ ，这样就可以确定序列 {}k x （k=0，1，2…）。如果 {}k x 有极限 *lim x x k k =∞→ ，由式 ,...1,0),(1==+k x x k k ρ 两边取极限可得 )(**x x ρ= ，可知 * x 为方程0)(=x f 的近似解。 2）Steffensen 迭代法的原理： 通过把改进的Aitken 方法应用于根据不动点迭代所得到的线性收敛序列，将收敛速度加速到二阶。 ()???? ?????+---===+k k k k k k k k k k k x y z x y x x y z x y 2) ()(21ρρ []x x x x x x x +---=)(2)(()()(2ρρρρψ 实验环境：操作系统：Windows 7； 实验平台：Turbo C++ 实验过程：写出算法→编写程序→调试运行程序→计算结果 1）简单迭代法的算法： Input:初始近似值x0,精度要求del,最大迭代次数N Output:近似解x 或失败信息 1. n ←1 2. While n ≤N do; 3. x ←f(x0); 4. if | x-x0| Newton迭代法求解非 线性方程 一、 Newton 迭代法概述 构造迭代函数的一条重要途径是用近似方程来代替原方程去求根。因此，如果能将非线性方程f (x )=0用线性方程去代替，那么，求近似根问题就很容易解决，而且十分方便。牛顿(Newton)法就是一种将非线性方程线化的一种方法。 设k x 是方程f (x )=0的一个近似根，把如果)(x f 在k x 处作一阶Taylor 展开，即： )x x )(x ('f )x (f )x (f k k k -+≈ (1-1) 于是我们得到如下近似方程： 0)x x )(x ('f )x (f k k k =-+ (1-2) 设0)('≠k x f ，则方程的解为： x ?=x k +f (x k ) f (x k )? (1-3) 取x ~作为原方程的新近似根1+k x ，即令： ) x ('f ) x (f x x k k k 1k -=+, k=0,1,2,… (1-4) 上式称为牛顿迭代格式。用牛顿迭代格式求方程的根的方法就称为牛顿迭代法，简称牛顿法。 牛顿法具有明显的几何意义。方程： )x x )(x ('f )x (f y k k k -+= (1-5) 是曲线)x (f y =上点))x (f ,x (k k 处的切线方程。迭代格式(1-4)就是用切线式(1-5)的零点来代替曲线的零点。正因为如此，牛顿法也称为切线法。 牛顿迭代法对单根至少是二阶局部收敛的，而对于重根是一阶局部收敛的。一般来说，牛顿法对初值0x 的要求较高，初值足够靠近*x 时才能保证收敛。若 要保证初值在较大范围内收敛，则需对)x (f 加一些条件。如果所加的条件不满足，而导致牛顿法不收敛时，则需对牛顿法作一些改时，即可以采用下面的迭代格式： ) x ('f ) x (f x x k k k 1k λ -=+, ?=,2,1,0k (1-6) 上式中，10<λ<，称为下山因子。因此，用这种方法求方程的根，也称为牛顿下山法。 牛顿法对单根收敛速度快，但每迭代一次，除需计算)x (f k 之外，还要计算 )x ('f k 的值。如果)x (f 比较复杂，计算)x ('f k 的工作量就可能比较大。为了避免计算导数值，我们可用差商来代替导数。通常用如下几种方法： 1. 割线法 如果用 1 k k 1k k x x ) x (f )x (f ----代替)x ('f k ，则得到割线法的迭代格式为： )x (f ) x (f )x (f x x x x k 1k k 1 k k k 1k --+---= (1-7) 2. 拟牛顿法 如果用 ) x (f )) x (f x (f )x (f k 1k k k ---代替)x ('f k ，则得到拟牛顿法的迭代格式为： )) x (f x (f )x (f ) x (f x x 1k k k k 2k 1k -+--- = (1-8) 3. Steffenson 法 如果用 ) x (f ) x (f ))x (f x (f k k k k -+代替)x ('f k ，则得到拟牛顿法的迭代格式为： ) x (f ))x (f x (f ) x (f x x k k k k 2k 1 k -+- =+ 不动点迭代法求解非线性方程组 摘要：一般非线性方程组可以写成()0F x =的形式，其中:n m F R R →是定义在区域 n D R ?上的向量函数。把方程组()0F x =改写成与之等价的形式：(x G x =)。因此，求 方程组()0F x =的解就转化为求函数的(G x )的不动点。本文首先介绍了多变量函数()F x 的微积分性质，接着介绍了用不动点迭代法求解非线性方程组。 关键词：多变量函数；微积分；不动点 Fixed Point Iteration Method For Solving Nonlinear Equations Abstract ：General nonlinear equations can be written in the form of ()F x θ=, where the vector function :n m F R R →is defined on the region n D R ?. Transform the equations ()0F x = into its equivalent form: (x G x =).Therefore, we can get the solution of ()0F x = by finding the fixed point of (G x ).In this paper, we first introduce some knowledge about multivariable calculus, then introduce the fixed point iteration method for solving nonlinear equations. Key words: multi-variable function; calculus; fixed point 1 引言 一般非线性方程组及其向量表示法:含有n 个方程的n 元非线性方程组的一般形式为 11221212(,...,)0(,...,)0......(,...,)0 n n m n f x x x f x x x f x x x =??=?? ? ?=? （1） 其中，()1,2,...,i f i n =是定义在n D R ?上的n 元实值函数，且i f 中至少有一个是非线性 函数。 用Gauss-Seidel迭代法解线性方程组的C语言源代码：#include ptr2->link=ptr1; ptr2=ptr1; } } j++; } if(ptr2!=NULL){ ptr2->link=NULL; ptr2=NULL; } if(Rhead!=NULL){ p1=(struct Line*)malloc(sizeof(Line)); if(p1==NULL){ printf("内存分配错误!\n"); exit(1); } p1->L=i; p1->head=Rhead; if(p2==NULL){ Lhead=p1; p2=p1; } else{ p2->next=p1; p2=p1; } } i++; Rhead=NULL; j=1; } if(p2!=NULL) p2->next=NULL; return Lhead; } struct Line *Change(struct Line*Lhead,int n){ struct Line*p1,*p2,*p3,*p; struct Row*ptr; int i=1,k,j; float max,t; if(Lhead==NULL){ printf("链表为空!\n");MATLAB代码 解线性方程组的迭代法
牛顿法求非线性方程的根
C++实现 牛顿迭代 解非线性方程组
牛顿法非线性方程求解
数值分析5-用Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组
matlab程序设计实践-牛顿法解非线性方程
求解线性方程组——超松弛迭代法(c)
二分法和牛顿法求解非线性方程(C语言)
Gauss-Seidel迭代法求解线性方程组
牛顿迭代法求解非线性方程组的代码
迭代法解线性方程组
基于Matlab的牛顿迭代法解非线性方程组
实验解线性方程组的基本迭代法实验
Gauss-Seidel迭代法求解线性方程组
非线性方程的简单迭代法和Steffensen迭代法
Newton迭代法求解非线性方程
不动点迭代法求解非线性方程组
迭代法解线性方程组(C语言描述)