牛頓迭代法解三元二次方程組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 / 非線性方程組中變量個(gè)數(shù)及方程個(gè)數(shù)const int N2=N*N; / jacobi矩陣的元素個(gè)數(shù)#d

2、efine eps 0.00001 / 收斂精度#define Max 2000 / 最大迭代次數(shù)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; /打開文件失敗則結(jié)束運(yùn)行ifstream fpin("d:Program FilesMicrosoft Visual StudioMyProjectsxieFCZX0.txt",ios:in);if(!fpin)cout<<"fpin open fail!"<<endl;return 0; /打開文件失敗則結(jié)束運(yùn)行 ofstream fpout2("d:Program FilesMicrosoft Visual StudioMyProjectsxieFCZanswer.txt"

4、;,ios:out); if(!fpout2)cout<<"fpout open fail!"<<endl;return 0; /打開文件失敗則結(jié)束運(yùn)行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 /計(jì)算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); /固定小數(shù)位數(shù)輸出fpout2<<setprecision(8)<<jacobiij<<"t"fpout2<<

7、endl;fpout2<<endl;aij=0;for(i=0;i<N;i+) /計(jì)算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<<"由計(jì)算機(jī)第 "<<it<<" 次隨機(jī)產(chǎn)生初值："&l

8、t;<endl;fpout1<<endl<<"由計(jì)算機(jī)第 "<<it<<" 次隨機(jī)產(chǎn)生初值："<<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<<"迭代次數(shù) x1tt x2tt x3"<<endl<<endl;doiter=iter+1;fpout2<<" "<<iter<<" "<<"t"for(i=0;i<N;i+)fpout2<<x1i<<"t"fpout2<<endl;/迭代格式y(tǒng)0=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，可讓方程同除某個(gè)數(shù)使其系數(shù)小于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