牛顿法求解非线性方程组的解C程序

牛顿法求解非线性方程组的解C程序(经调试)

/*############################################################################*/
/*#### 本函数用用牛顿法求解非线性方程组的解 #########*/
/*#### 用户在子函数function中输入方程组的信息，即依次输入fi(x) #########*/
/*#### 用户在子函数Jacobi中输入Jacobi矩阵的信息，依次输入行与列 #########*/
/*#### 用户按提示在程序执行的过程中输入方程的维数n，以及初始向量X0#########*/
/*#### 本程序中迭代精度设为10^(-7),用户可在主程序中改动此值 #########*/
/*#### 程序输出为满足精度要求的方程组的解 #########*/
/*#### 本程序已经调试，作者：沈欢，邮箱:holyshen@http://doc.guandang.net #########*/
/*############################################################################*/

#include <stdio.h>
#include <math.h>
void function(double *x,double *f);//本子函数用于计算向量f(x)的值
//入口参数为方程x指针向量,f(x)向量指针

void Jacobi(double *x,double **J);//本子函数用于计算Jacobi矩阵的值
//入口参数为向量x的指针和J矩阵的指针

void LU(int n,double *b,double ** A,double *detax);//本子函数用于计算线性方程组Ax=b的解，用LU分解法实现
//入口参数为方程维数n,右端项b,系数矩阵A
//detax用于存储方程组的解
void main()
{
int n,i,j,breakindex;
double num;
double error=0.0000001;//定义精度为10^(-7)


//***********信息输入***************//
printf("本函数用用牛顿法求解非线性方程组的解\n");
printf("用户在子函数function中输入方程组的信息，即依次输入fi(x)\n");
printf("用户在子函数Jacobi中输入Jacobi矩阵的信息，依次输入行与列\n");
printf("本程序中迭代精度设为10^(-7),用户可在主程序中改动此值\n");
printf("现在请输入非线性方程的维数n(方程个数):(按Enter键继续)\n");
scanf("%d",&n);

double *x=new double [n+1];//开辟解的存储空间，用于存储xk向量

printf("\n");
printf("请依次输入初始迭代向量X0的各元素值，并按Enter键继续\n");
for(i=1;i<=n;i++)//输入A矩阵
{
scanf("%lf",&num);
x[i]=num;
}

//**************************************//

//**************迭代求解****************//
double *f=new double [n+1];//开辟空间用于存储f(xk)向量
double *detax=new doub
le [n+1];//开辟空间用于存储△xk向量

double **J = new double * [n+1];//为Jacobi矩阵开辟空间
for(i=0;i<n+1;i++)J[i]=new double [n+1];


for(i=1;;i++)//进行迭

牛顿法求解非线性方程组的解C程序(经调试)

代
{
function(x,f);//计算此时的f(x)向量

breakindex=0;//设置跳出标志
for(j=1;j<=n;j++)//判断迭代是否结束
{
if(abs(f[j])<error)breakindex=breakindex+1;
}
if(breakindex==n)break;

Jacobi(x,J);//计算此时的Jacobi矩阵
LU(n,f,J,detax);//用LU分解法求△x

for(j=1;j<=n;j++)//xk+1=xk+△x
{
x[j]=x[j]-detax[j];
}

}

//***************************************//



//*************信息输出*****************//
printf("经过 %d 次迭代，此非线性方程组的解如下：\n",i);
for(i=1;i<=n;i++)
{
printf("%lf\n",x[i]);
}
printf("欢迎您下次使用\n");
//**************************************//



}









void function(double *x,double *f)//本子函数用于计算向量f(x)的,入口参数为x向量指针及f(x)向量指针
{
f[1]=pow(x[1],2)+pow(x[2],2)+pow(x[3],2)-1;
f[2]=2*pow(x[1],2)+pow(x[2],2)-4*x[3];
f[3]=3*pow(x[1],2)-4*x[2]+pow(x[3],2);
}

void Jacobi(double *x,double **J)//本子函数用于计算Jacobi矩阵的值，入口参数为向量x的指针和J矩阵的指针
{
J[1][1]=2*x[1];
J[1][2]=2*x[2];
J[1][3]=2*x[3];
J[2][1]=4*x[1];
J[2][2]=2*x[2];
J[2][3]=-4;
J[3][1]=6*x[1];
J[3][2]=-4;
J[3][3]=2*x[3];
}




void LU(int n,double *b,double ** A,double *detax)//本子函数用于计算线性方程组Ax=b的解，用LU分解实现，入口参数为方程维数n,右端项b,系数矩阵A，其中detax用于存储方程组的解
{
int i,j,k;
double sum;


double *y=new double [n+1];//过渡向量


double **L = new double * [n+1];//为L矩阵开辟空间
for(i=0;i<n+1;i++)L[i]=new double [n+1];

double **U = new double * [n+1];//为U矩阵开辟空间
for(i=0;i<n+1;i++)U[i]=new double [n+1];



//****对A进行LU分解****//
for(i=1;i<=n;i++)
{
L[i][i]=1;
}
for(j=1;j<=n;j++)
{
U[1][j]=A[1][j];
}
for(i=2;i<=n;i++)
{
L[i][1]=A[i][1]/A[1][1];
}

for(k=2;k<=n;k++)
{
for(j=k;j<=n;j++)//计算U[k][j]
{
sum=0;
for(i=1;i<=k-1;i++)
{
sum=sum+L[k][i]*U[i][j];
}
U[k][j]=A[k][j]-sum;
}

if(k<n)//在k<n时计算L[i][k]
{
for(i=k+1;i<=n;i++)
{
sum=0;
for(j=1;j<=k-1;j++)
{
sum=sum+L[i][j]*U[j][k];
}
L[i][k]=(A[i][k]-sum)/U[k][k];
}

}


}

//************************//


//********回带求解*******//

y[1]=b[1];//计算y
for(i=2;i<=n;i++)
{
j=1;
sum=0;
while(j<i)
{
sum=sum+L[i][j]*y[j];
j++;
}
y[i]=b[i]-sum;
}

detax[n]=y[n]/U[n][n];//计算△x
for(i=n-1;i>=1;i--)
{
j=n;

sum=0;
while(j>i)
{
sum=sum+U[i][j]*detax[j];
j--;
}
detax[i]=(y[i]-sum)/U[i][i];

}


//************************//



}