nonlinearequation.cpp

通过此程序能够得到非线性方程组的数值解

//拟牛顿(QuasiNewton)法求解非线性方程组一组实根
#include "Comm.h"		//公共部分头文件
#include "Matrix.h"		//矩阵类及相关函数等的定义
#include "LinearEquation.h"	//线性方程(组)求解头文件
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<string.h>
#include <iostream>		//输入输出流头文件

const int  N=2;         //方程组变量的个数
const int  M=20;        //迭代的最大次数
const double eps=1e-5;  //eps为精度

template <class _Ty> 
void FunctionValueRSN(valarray<_Ty>& x, valarray<_Ty>& y)
{//输入方程组
/*
	x^2+y^2-5=0
	(x+1)*y-(3*x+1)=0
*/
	y[0]=x[0]*x[0]+x[1]*x[1]-5;
	y[1]=(x[0]+1)*x[1]-(3*x[0]+1);
}


template <class _Ty>
inline int 
RootQuasiNewton(_Ty eps, _Ty t, _Ty h, valarray<_Ty>& x, int k)
{
	 int i, j, l;
    _Ty am, z, beta, d;
	
	int n = x.size();		//方程组中各方程个数，也是未知数个数

	matrix<_Ty> a(n, n);	//定义二维数组对象a
	valarray<_Ty> b(n);		//定义一维数组对象b
	valarray<_Ty> y(n);		//定义一维数组对象y
	    
	l = k;
	am = 1.0 + eps;
    while(am>=eps)
    {
		FunctionValueRSN(x, b);
        am = 0.0;
        for(i=0; i<n; i++)
        {
			z = Abs(b[i]);
            if(z>am) am = z;
        }
        if(am>=eps)
        {
			l--;
            if(l==0)
			{
				cout << "Fail(0)!" << endl;
				return(0);
			}

            for(j=0; j<n; j++)
            {
				z = x[j];
				x[j] += h;
                FunctionValueRSN(x, y);
                for(i=0; i<n; i++) 
					a(i,j) = y[i];
                x[j] = z;
            }

			if(LE_TotalChoiceGauss(a, b) == 0)
			{
				cout << "Fail(-1)!" << endl;
				return(-1);					//矩阵a奇异
			}

            beta = 1.0;

            for(i=0; i<n; i++) beta -= b[i];
            
			if(FloatEqual(beta,0))
			{
				cout << "Fail(-2)!" << endl;
				return(-2);					//beta=0
			}

            d = h / beta;
            
			for(i=0; i<n; i++) x[i]=x[i]-d*b[i];

            h = t * h;
        }
    }
    return(k-l);					//正常结束，返回迭代次数
}

void main()
{
	double a[N] = {1, 1};              //在(1，1)附近求方程的解,维数“2”是方程组变量的个数
	valarray<double> va(a, N);
	double t=0.01,h=0.001;  
	RootQuasiNewton(eps, t, h, va,  M);
	cout<< "方程组的解为 :  ";
	ValarrayPrint(va);
}