牛顿法解非线性方程组

先把方程组给定到一个矩阵，带入初始值求解，再求解雅可比矩阵，带值求解，
最后用公式求解值，直到计算结束。
牛顿法解非线性方程组源代码
//#include
//#include
//#include
//#define N 2 // 非线性方程组中方程个数、未知量个数
//#define E 0.0001 // 差向量1范数的上限
//#define Max 100 //最大迭代次数
using namespace std;
const int N2 = 2 * N;
int main()
{
void ff(float xx[N], float yy[N]);
void ffjacobian(float xx[N], float yy[N][N]);
void inv_jacobian(float yy[N][N], float inv[N][N]);
void newdundiedai(float x0[N], float inv[N][N], float y0[N], float x1[N]);
float x0[N] = { 2.0, 0.25 }, y0[N], jacobian[N][N], invjacobian[N][N], x1[N], errornorm;
int i, count = 0;
//for( i=0;i < N;i++)
// cin>>x0[i];
cout << “初始解向量是：” << endl;
for (i = 0; i < N; i++)
cout << x0[i] << ” “;
cout << endl; cout << endl;
do
{
count = count + 1;
cout << “第 ” << count << ” 次迭代开始” << endl << endl;
//计算向量函数值—因变量向量存放在一维数组y0中，y0做函数ff的实参，ff(x0,y0)的调用将传出因变量向量
ff(x0, y0); //计算雅克比矩阵，二维数组jacobian作为ffjacobian的实参;调用ffjacobian(x0,jacobian)，
//按照地址传送方式，二维数组jacobian将传出雅克比矩阵的计算结果
ffjacobian(x0, jacobian); //调用函数inv_jacobian计算雅克比矩阵的逆矩阵，计算结果存放在实参invjacobian中
inv_jacobian(jacobian, invjacobian); //由近似解向量x0迭代计算近似解向量x1
cout << “第 ” << count << ” 次迭代得到”;
newdundiedai(x0, invjacobian, y0, x1); //计算x1与x0的差向量的1范数
errornorm = 0;
for (i = 0; i < N; i++)
errornorm = errornorm + fabs(x1[i] - x0[i]);
if (errornorm < E)
break; //为下次循环计算准备新的x0
for (i = 0; i < N; i++)
x0[i] = x1[i];
}
while (count < Max);
cin >> count;
return 0;
}
void ff(float xx[N], float yy[N])
{
float x, y;
int i;
x = xx[0];
y = xx[1];
yy[0] = x*x - 2 * x - y + 1/2;
yy[1] = x*x + 4 * y*y - 4;
cout << “当前自变量向量为：” << endl;
for (i = 0; i < N; i++)
{
cout << setw(10) << setiosflags(ios::right) << setprecision(5) << setiosflags(ios::fixed) << xx[i] << ” “;
cout << endl;
cout << endl;
cout << “当前因变量向量为：” << endl;
}
for (i = 0; i < N; i++)
{
cout << setw(10) << setiosflags(ios::right) << setprecision(5) << setiosflags(ios::fixed) << yy[i] << ” “;
cout << endl;
cout << endl;
}
}
void ffjacobian(float xx[N], float yy[N][N])
{
float x, y;
int i, j;
x = xx[0];
y = xx[1];
yy[0][0] = 2 * x - 2;
yy[0][1] = -1;
yy[1][0] = 2 * x;
yy[1][1] = 8 * y;
cout << “当前的雅克比矩阵是： ” << endl;
for (i = 0; i < N; i++)
{
for (j = 0; j < N; j++)
cout << setw(10) << setiosflags(ios::right) << setprecision(5) << setiosflags(ios::fixed) << yy[i][j] << ” “;
cout << endl;
}
cout << endl;
}
void inv_jacobian(float yy[N][N], float inv[N][N])
{
float aug[N][N2], L;
int i, j, k;
cout << “计算雅克比矩阵的逆矩阵” << endl;

for (i = 0; i<N; i++)    //构造aug{    for (j = 0; j<N; j++)        aug[i][j] = yy[i][j];    for (j = N; j<N2; j++)    if (j == i + N) aug[i][j] = 1;    else  aug[i][j] = 0;}for (i = 0; i<N; i++)                 //输出aug{    for (j = 0; j<N2; j++)        cout << setw(10) << setiosflags(ios::right) << setprecision(5) << setiosflags(ios::fixed) << aug[i][j] << "  ";        cout << endl;}cout << endl;for (i = 0; i<N; i++)              //在aug上实施初等行变换{    for (k = i + 1; k<N; k++)    {        L = -aug[k][i] / aug[i][i];        for (j = i; j<N2; j++)            aug[k][j] = aug[k][j] + L*aug[i][j];    }}for (i = 0; i<N; i++){    for (j = 0; j<N2; j++)        cout << setw(10) << setiosflags(ios::right) << setprecision(5) << setiosflags(ios::fixed) << aug[i][j] << "  ";        cout << endl;}cout << endl;for (i = N - 1; i>0; i--){    for (k = i - 1; k >= 0; k--)    {        L = -aug[k][i] / aug[i][i];        for (j = N2 - 1; j >= 0; j--)            aug[k][j] = aug[k][j] + L*aug[i][j];    }}for (i = 0; i<N; i++){    for (j = 0; j<N2; j++)        cout << setw(10) << setiosflags(ios::right) << setprecision(5) << setiosflags(ios::fixed) << aug[i][j] << "  ";    cout << endl;}cout << endl;for (i = N - 1; i >= 0; i--)      //对aug的第i行每个元素除aug[i][i]{    for (j = N2 - 1; j >= 0; j--)        aug[i][j] = aug[i][j] / aug[i][i];}for (i = 0; i<N; i++)     //从aug中提取逆矩阵{    for (j = 0; j<N2; j++)        cout << setw(10) << setiosflags(ios::right) << setprecision(5) << setiosflags(ios::fixed) << aug[i][j] << "  ";    cout << endl;    for (j = N; j<N2; j++)        inv[i][j - N] = aug[i][j];}cout << endl;cout << "雅克比矩阵的逆矩阵是： " << endl;for (i = 0; i<N; i++){    for (j = 0; j<N; j++)        cout << setw(10) << setiosflags(ios::right) << setprecision(5) << setiosflags(ios::fixed) << inv[i][j] << "  ";    cout << endl;}cout << endl;

}
void newdundiedai(float x0[N], float inv[N][N], float y0[N], float x1[N])
{
int i, j;
float sum = 0;
for (i = 0; i < N; i++)
{
sum = 0;
for (j = 0; j < N; j++)
sum = sum + inv[i][j] * y0[j];
x1[i] = x0[i] - sum;
}
cout << “近似解向量” << endl;
for (i = 0; < N; i++)
cout << x1[i] << ” “;
cout << endl; cout << endl;
}