// 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;//如果有一个误差不符合要求即退出
}
}