牛顿法解非线性方程组

 

// 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 J[N][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 X[N],double J[N][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,J,Wuchaxian);
 }
 F(x,X);
 for(i=0;i<N;i++)//输出方程解
 {
  printf("x(%d)=%e F(%d)=%e\n",i,x[i],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,j,k,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 x[N],double J[N][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;//如果有一个误差不符合要求即退出
 }
}