(完整word版)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的初始化，撤销下面两行的注释符，

就可以由键盘向x0读入初始近似解向量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<

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<

}

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<

}