拟牛顿法求解非线性方程组

//ÎļþÃûNNEWTON.C//ÄâÅ£¶Ù·¨Çó½â·ÇÏßÐÔ·½³Ì×é#includestdio.h#includestdlib.h#includemath.hvoidn_newton(n,eps,x,f)//ÄâÅ£¶Ù·¨intn;doubleeps,x[];void(*f)(double[],double[]);{inti,j,k=0;doubleam,z,beta,d;doublet=0.1,h=0.1;double*a,*b,*y;intgauss(int,double[],double[]);a=(double*)malloc(n*n*sizeof(double));b=(double*)malloc(n*sizeof(double));y=(double*)malloc(n*sizeof(double));am=1.0+eps;while(am=eps){(*f)(x,y);am=0.0;for(i=0;i=n-1;i++){b[i]=y[i];z=fabs(y[i]);if(zam)am=z;}if(am=eps){k=k+1;if(k500){printf(µü´ú³¬¹ý500´Î,¿ÉÄܲ»ÊÕÁ²!\n);free(a);free(b);free(y);return;}for(j=0;j=n-1;j++){z=x[j];x[j]=x[j]+h;(*f)(x,y);for(i=0;i=n-1;i++)a[i*n+j]=y[i];x[j]=z;}if(gauss(n,a,b)==0){printf(ϵÊý¾ØÕóÆæÒì,GaussÏûÈ¥·¨Ê§°Ü!\n);free(a);free(b);free(y);return;}beta=1.0;for(i=0;i=n-1;i++)beta=beta-b[i];if(fabs(beta)+1.0==1.0){printf(Beta=0.³ÌÐò¹¤×÷ʧ°Ü£¡\n);free(a);free(b);free(y);return;}d=h/beta;for(i=0;i=n-1;i++)x[i]=x[i]-d*b[i];h=t*h;}}printf(ʵ¼Êµü´ú´ÎÊýΪ:%d\n,k);free(a);free(b);free(y);return;}intgauss(intn,doublea[],doubleb[])//È«Ñ¡Ö÷ÔªGaussÏûÈ¥·¨{int*js,k,i,j,is;doubled,t;js=(int*)malloc(n*sizeof(int));for(k=0;k=n-2;k++){d=0.0;for(i=k;i=n-1;i++)for(j=k;j=n-1;j++){t=fabs(a[i*n+j]);if(td){d=t;js[k]=j;is=i;}}if(d+1.0==1.0){free(js);return(0);}if(js[k]!=k)for(i=0;i=n-1;i++){t=a[i*n+k];a[i*n+k]=a[i*n+js[k]];a[i*n+js[k]]=t;}if(is!=k){for(j=k;j=n-1;j++){t=a[k*n+j];a[k*n+j]=a[is*n+j];a[is*n+j]=t;}t=b[k];b[k]=b[is];b[is]=t;}d=a[k*n+k];for(j=k+1;j=n-1;j++)a[k*n+j]=a[k*n+j]/d;b[k]=b[k]/d;for(i=k+1;i=n-1;i++){for(j=k+1;j=n-1;j++)a[i*n+j]=a[i*n+j]-a[i*n+k]*a[k*n+j];b[i]=b[i]-a[i*n+k]*b[k];}}d=a[n*n-1];if(fabs(d)+1.0==1.0){free(js);return(0);}b[n-1]=b[n-1]/d;for(i=n-2;i=0;i--){t=0.0;for(j=i+1;j=n-1;j++)t=t+a[i*n+j]*b[j];b[i]=b[i]-t;}js[n-1]=n-1;for(k=n-1;k=0;k--)if(js[k]!=k){t=b[k];b[k]=b[js[k]];b[js[k]]=t;}free(js);return(1);}//°üº¬ÔÚÎļþNNEWTON.CÖÐmain(){intk;doublex[3]={1.0,1.0,1.0};voidf(double[],double[]);n_newton(3,0.0000001,x,f);for(k=0;k=2;k++)printf(x(%d)=%f\n,k,x[k]);return;}voidf(doublex[],doubley[])//¼ÆËã·½³Ì×é¸÷×ó¶Ëº¯ÊýÖµ{y[0]=x[0]*x[0]+x[1]*x[1]+x[2]*x[2]-1.0;y[1]=2*x[0]*x[0]+x[1]*x[1]-4*x[2];y[2]=3*x[0]*x[0]-4*x[1]+x[2]*x[2];}