牛顿迭代法解三元二次方程组C

1、牛顿迭代法解三元二次方程组（C+版） *方程组为：Y1=-X1+0.3X2*X3-X32+0.6=0Y2=-0.1X12-X2 +0.8X1*X3-X3+0.4=0Y3=0.3X1-0.5X22 +0.7X1*X2-X3+0.5=0 * #include<fstream>#include<iostream>#include <iomanip>#include<cmath>#include <cstdlib>#define N 3 / 非线性方程组中变量个数及方程个数const int N2=N*N; / jacobi矩阵的元素个数#d

2、efine eps 0.00001 / 收敛精度#define Max 2000 / 最大迭代次数using namespace std;double main()double x0N,yN,x1N,es,esmax,jacobiNN,aij;int i,j,k,it=0,iter=0; ofstream fpout1("d:Program FilesMicrosoft Visual StudioMyProjectsxieFCZx0.txt",ios:out); if(!fpout1)cout<<"fpout open fail!"<&

3、lt;endl;return 0; /打开文件失败则结束运行ifstream fpin("d:Program FilesMicrosoft Visual StudioMyProjectsxieFCZX0.txt",ios:in);if(!fpin)cout<<"fpin open fail!"<<endl;return 0; /打开文件失败则结束运行 ofstream fpout2("d:Program FilesMicrosoft Visual StudioMyProjectsxieFCZanswer.txt"

4、;,ios:out); if(!fpout2)cout<<"fpout open fail!"<<endl;return 0; /打开文件失败则结束运行cout<<"* Please input X0 *"<<endl; / 技巧一：为保证jacobi矩阵尽快收敛，初值最好在 0 附近for(i=0;i<N;i+)cin>>x0i;cout<<"输入的初值为："<<endl;for(i=0;i<N;i+)cout<<x0i<

5、;<"t"cout<<endl;fpout1<<"输入的初值为："<<endl;do it+;for(i=0;i<N;i+)fpout1<<x0i<<"t"fpout1<<endl;/jacobi have N2 element /计算jacobi矩阵的值jacobi00=0;jacobi01=0.3*x02;jacobi02=0.3*x01-2*x02;jacobi10=-0.2*x00+0.8*x02;jacobi11=0;jacobi12=0.8

6、*x00-1;jacobi20=0.3+0.7*x01;jacobi21=-x01+0.7*x00;jacobi22=0;fpout2<<"第 "<<it<<" 次jacobi矩阵是："<<endl;for(i=0;i<N;i+) for(j=0;j<N;j+) fpout2.setf(cout.showpoint); /固定小数位数输出fpout2<<setprecision(8)<<jacobiij<<"t"fpout2<<

7、endl;fpout2<<endl;aij=0;for(i=0;i<N;i+) /计算jacobi矩阵元素的平方和aijfor(j=0;j<N;j+)aij=aij+jacobiij*jacobiij;if(aij>1)cout<<"Sorry, aij>1"<<endl;k=rand()%3; x0k=(double)rand()/RAND_MAX;cout<<endl<<"由计算机第 "<<it<<" 次随机产生初值："&l

8、t;<endl;fpout1<<endl<<"由计算机第 "<<it<<" 次随机产生初值："<<endl;cout<<" "for(i=0;i<N;i+) cout<<x0i<<"t" cout<<endl; while(aij>1); /保证aij<1for(i=0;i<N;i+)x1i=x0i;fpout2<<"* 牛顿迭代如下 *"<

9、<endl;fpout2<<"迭代次数 x1tt x2tt x3"<<endl<<endl;doiter=iter+1;fpout2<<" "<<iter<<" "<<"t"for(i=0;i<N;i+)fpout2<<x1i<<"t"fpout2<<endl;/迭代格式y0=0.3*x11*x12-x12*x12+0.6; y1=-0.1*x10*x10+0.8*

10、x10*x12-x12+0.4;y2=0.3*x10-0.5*x11*x11+0.7*x10*x11+0.5;/技巧二：由于要使最终的aij小于0，可让方程同除某个数使其系数小于1esmax=0.0;for(i=0;i<N;i+)es=yi-x1i;if(fabs(es)>fabs(esmax)esmax=es;if(fabs(esmax)<eps)fpout2<<endl<<endl<<"方程组的解为 : "<<endl;cout<<endl<<endl<<"方程组的解为 : "<<endl;for(i=0;i<N;i+)cout<<x1i<<"t"fpout2<<x1i<<"t"cout<<endl<<endl;fpout2<<endl;break;for(i=0;i<N;i+)x1i=yi;while(iter<Max);fpin.close();fpo