cpp4.cpp

牛顿法解非线性方程组输入系数可直接运行的

#include <iostream.h>
#include <stdlib.h>
#include <math.h>
#include <conio.h>
#define f0(x1,x2) (x1+2*x2-3)
#define f1(x1,x2) (2*x1*x1+x2*x2-5)
#define x_ 0.000001
#define matrixNum 2
double *matrixF2(double *x);

int y=0;

void main()
{
    int i,j;
    double p,x[2];
    double *b;
	double *matrixF;     //矩阵F
    double *matrixF_;   //矩阵F的雅可比矩阵的逆 
    b=(double *)malloc(matrixNum);
	matrixF=(double *)malloc(matrixNum);
    matrixF_=(double *)malloc(matrixNum*matrixNum);
    cout<<"请输入初值:";
	for(i=0;i<matrixNum;i++)
        cin>>*(x+i);
    do
	{
	    p=0.0;
		for(i=0;i<matrixNum;i++)
            *(b+i)=0;
		*matrixF=f0(*x,*(x+1));
        *(matrixF+1)=f1(*x,*(x+1));
        matrixF_=matrixF2(x);
		for(i=0;i<matrixNum;i++)
		{
			for(j=0;j<matrixNum;j++)
				*(b+i)+=*(matrixF_+i*matrixNum+j)*(*(matrixF+j));//矩阵F和矩阵F的雅可比矩阵的逆相乘
            *(x+i)=*(x+i)-*(b+i);       //                         //*(matrixF+j)为矩阵F的列数
			cout<<*(x+i)<<" ";
		}
	    cout<<endl;
        for (i=0;i<matrixNum;i++)
            p+=pow(*(b+i),2);   //pow(x,y)=x.^y
		y++;
	}while(sqrt(p)>x_);     //while为0时结束循环
	cout<<"停止迭代，最终迭代结果为"<<*x<<','<<*(x+1)<<""<<endl;
    free(matrixF);
    free(matrixF_);
	//delete [] matrixF;
    //delete [] matrixF_;
    getch();
}


double *matrixF2(double *x)
{
    int i,j;
    double t; 
	double *matrixF1;   //矩阵F的雅可比矩阵
    double *matrixF2;   //矩阵F的雅可比矩阵的逆 
    matrixF1=(double *)malloc(matrixNum*matrixNum);
    matrixF2=(double *)malloc(matrixNum*matrixNum);
    for(i=0;i<matrixNum;i++)
        for(j=0;j<matrixNum;j++)
			if(i==j)
              *(matrixF2+i*matrixNum+j)=1;
			else *(matrixF2+i*matrixNum+j)=0;
    *matrixF1=(f0((*x+x_),*(x+1))-f0(*x,*(x+1)))/x_;
    *(matrixF1+1)=(f0(*x,(*(x+1)+x_))-f0(*x,*(x+1)))/x_;
    *(matrixF1+2)=(f1((*x+x_),*(x+1))-f1(*x,*(x+1)))/x_;
    *(matrixF1+3)=(f1(*x,(*(x+1)+x_))-f1(*x,*(x+1)))/x_;
    //for(i=0;i<matrixNum;i++)
    //   cout<<*(x+i)<<endl;  
	for(i=0;i<matrixNum;i++)    
	{
		for(j=0;j<matrixNum;j++)
	       cout<<*(matrixF1+i*matrixNum+j)<<" ";
        cout<<endl;
	}
	//求矩阵F的雅可比矩阵的逆 
	t=*matrixF1;
    for(i=0,j=0;j<matrixNum;j++)           
	{
		*(matrixF1+i*matrixNum+j)/=t;
        *(matrixF2+i*matrixNum+j)/=t;
	}
    t=*(matrixF1+1*matrixNum);
	for(i=1,j=0;j<matrixNum;j++)
	{
		*(matrixF1+i*matrixNum+j)-=*(matrixF1+j)*t;
        *(matrixF2+i*matrixNum+j)-=*(matrixF2+j)*t;
	}
	t=*(matrixF1+1*matrixNum+1);
    for(i=1,j=0;j<matrixNum;j++)
	 {
		 *(matrixF1+i*matrixNum+j)/=t;
         *(matrixF2+i*matrixNum+j)/=t;
	}  
	t=*(matrixF1+1);
	for(i=i,j=0;j<matrixNum;j++)
	{	 
	    *(matrixF1+j)-=*(matrixF1+i*matrixNum+j)*t;
        *(matrixF2+j)-=*(matrixF2+i*matrixNum+j)*t;
	}
    for(i=0;i<matrixNum;i++)           
	{
		for(j=0;j<matrixNum;j++)
	       cout<<*(matrixF1+i*matrixNum+j)<<" ";
        cout<<endl;
	}
    for(i=0;i<matrixNum;i++)           
	{
		for(j=0;j<matrixNum;j++)
	       cout<<*(matrixF2+i*matrixNum+j)<<" ";
        cout<<endl;
	}
	cout<<"第"<<y<<"次迭代结果为"<<*x<<','<<*(x+1)<<""<<endl;
	getch();
    return matrixF2;
     free( matrixF1 );
     free( matrixF2 );
	//delete [] matrixF1;
    //delete [] matrixF2;
}