牛顿法解非线性方程组

// Newton.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include<iostream.h>
#include<math.h>
#define N 3
void get_Chuzhi(double a[N],int n);
void F(double x[N],double X[N]);
void get_J(double x[N],double J[N][N]);
void get_J_1(double J[N][N]);
void Get_Flash_x(double X[N],double x[N],double J[N][N],double Wuchanxian);
void Newton(double x[N],double X[N],double Wuchaxian);
void main()
{
?double Wuchaxian=1.0,x[N],X[N];//定义定义：误差"Wuchaxian"、未知数x值"x[N]"、函数值"X[N]"、
?double J[N][N];//雅阁比矩阵"J[N][N]"(后面程序将雅阁比逆阵赋给J[N][N]
?Newton(x,X,J,Wuchaxian);//进入迭代过程
}

void Newton(double x[N],double Wuchaxian)////////////////解非线性方程组
{
?int i;//控制循环最高次数
?for(i=0;i<100&&Wuchaxian>0.0000001;i++)//迭代过程
?{
??if(i==0) {get_Chuzhi(x,N);}
??F(x,X);
??get_J(x,J);
??get_J_1(J);
??Get_Flash_x(X,x,Wuchaxian);
?}
?F(x,X);
?for(i=0;i<N;i++)//输出方程解
?{
??printf("x(%d)=%e F(%d)=%en",i,x[i],X[i]);
?}
}

void get_Chuzhi(double a[],int n)///////////////////////////////////////////////////录入方程初值
{
?int i;
?printf("请输入%d个初值n",N);
?for(i=0;i<N;i++)
?{
??cin>>a[i];
??if(i==0){cout<<"您输入的初值为："<<endl;}
??cout<<a[i]<<" ";
?}
?cout<<endl;
?cout<<endl;
}

void F(double x[N],double X[N])/////////////////////////////////////////////////////求方程函数值F(x)
{
?/////////////////////////
?/////////////////////////
?X[0]=x[0]*x[0]+x[1]*x[1]+x[2]*x[2]-2.0;//方程可替换
?X[1]=3.0*x[0]*x[0]+x[1]*x[1]-3.0*x[2];
?X[2]=2.0*x[0]*x[0]-3.0*x[1]+2.0*x[2]*x[2];
?/////////////////////////
?/////////////////////////
}

void get_J(double x[N],double J[N][N])//////////////////////////////////////////////求雅各比矩阵
{
?/////////////////////////
?/////////////////////////
?J[0][0]=2.0*x[0];J[0][1]=2.0*x[1];J[0][2]=2.0*x[2];//方程偏导函数可替换
?J[1][0]=6.0*x[0];J[1][1]=2.0*x[1];J[1][2]=-3.0;
?J[2][0]=4.0*x[0];J[2][1]=-3.0;J[2][2]=4.0*x[2];
?/////////////////////////
?/////////////////////////
}

void get_J_1(double D[N][N])////////////////////////////////////////////////////////求雅各比逆阵（通过建立单位增广矩阵） { ?int i,j,k,n,s[N];//i,n控制循环,s[N]存放矩阵行编号 ?double P,Q,T,S[N][N];//P,Q存放主元所在行数据，上消、下消、过程需要乘的系数； ?????? //T 数据互换时需要的中间变量；S[N][N]是单位阵 ?for(i=0;i<N;i++)//建立一个同阶单位阵 ??for(j=0;j<N;j++) ??{ ???s[i]=i; ???if(i==j) S[i][j]=1; ???else S[i][j]=0; ??} ?for(i=0;i<N;i++)//列主元消去 ?{ ??for(j=i;j<N;j++)//选列主元 换行过程 ??{ ???if(fabs(D[i][i])<fabs(D[j][i])) ????for(k=0;k<N;k++) ????{ ?????T=D[i][k];D[i][k]=D[j][k];D[j][k]=T; ?????T=S[i][k];S[i][k]=S[j][k];S[j][k]=T;//行所有数据的互换 ?????s[i]=j;s[j]=i;//记住行 ????} ??} ??for(j=i+1,n=i-1;;n--,j++)//消去过程 ??{ ???for(k=0;k<N;k++) ???{ ????if(k==0) P=D[j][i]/D[i][i];//P、Q是消去时 需要乘的系数 ????if(j<N) S[j][k]=S[j][k]-S[i][k]*P;//下消 ????if(k==0) Q=D[n][i]/D[i][i]; ????if(n>=0) S[n][k]=S[n][k]-S[i][k]*Q;//上消 ???} ???for(k=i;k<N;k++) ???{ ????if(k==0) P=D[j][i]/D[i][i]; ????if(j<N) D[j][k]=D[j][k]-D[i][k]*P; ????if(k==0) Q=D[n][i]/D[i][i]; ????if(n>=0) D[n][k]=D[n][k]-D[i][k]*Q;//逆阵做相同的消去 ???} ???if(n<0 && j>=N) break;//直到将上消下消完成后退出循环 ??} ?} ?for(i=0;i<N;i++)//行还原 ??for(j=0;j<N;j++) ??{ ???if(s[i]==j) ????for(k=0;k<n;k++) ????{ ?????T=D[i][k];D[i][k]=D[j][k];D[j][k]=T; ?????T=S[i][k];S[i][k]=S[j][k];S[j][k]=T;//行所有数据互换 ????} ??} ?for(i=0;i<N;i++)//逆阵每行除以原矩阵对角线上的数据 ??for(j=0;j<N;j++) ??{ ???S[i][j]=S[i][j]/D[i][i]; ??} ?for(i=0;i<N;i++)//将逆阵赋给D带回 ??for(j=0;j<N;j++) ??{ ???D[i][j]=S[i][j]; ??} } void Get_Flash_x(double X[N],double Wuchaxian)///////////刷新x数值 { ?double DX=0.0;//定义DX存放△X ?int i,j;//控制循环 ?for(i=0;i<N;i++) ?{ ??for(j=0;j<N;j++) ??{ ???DX=DX+J[i][j]*X[j];//计算DX数值 ??} ??x[i]=x[i]-DX; ?} ?for(i=0;i<N;i++) ?{ ??Wuchaxian=fabs(X[i]); ??if(Wuchaxian>0.0000001) break;//如果有一个误差不符合要求即退出 ?} }