newton3.cpp

数值分析中牛顿迭代算法的C++实现

#include<iostream.h>
#include<math.h>
#include<stdio.h>
#include <iomanip.h>
//选出数组中的最大值
double max(double a[],int n)
{int i;
double max;
max=fabs(a[0]);
for (i=0;i<=n-1;i++)
{if (fabs(a[i])>max)
    max=fabs(a[i]);
}
return (max);
}
void main()
{double X[4]={0,0,0,0},err=1,x,y,F1[4][4],F[4],xx[4],c[4][4],d[4],sum1,m[4];
int i,j,k=0,n=4,ik,i1,j1,flag=0;
for (i=0;i<=10;i++)
	for(j=0;j<=20;j++)
	{x=0.08*i;
	y=0.5+0.05*j;
	err=1;
	for (i1=0;i1<4;i1++)
	X[i1]=2;
	while (err>1e-12)
		{F[0]=-(0.5*cos(X[0])+X[1]+X[2]+X[3]-x-2.67);
		F[1]=-(X[0]+0.5*sin(X[1])+X[2]+X[3]-y-1.07);
		F[2]=-(0.5*X[0]+X[1]+cos(X[2])+X[3]-x-3.74);
		F[3]=-(X[0]+0.5*X[1]+X[2]+sin(X[3])-y-0.79);   /*对数组F赋值*/
		F1[0][0]=-0.5*sin(X[0]);
		F1[0][1]=1;
		F1[0][2]=1;
		F1[0][3]=1;
		F1[1][0]=1;
		F1[1][1]=0.5*cos(X[1]);
		F1[1][2]=1;
		F1[1][3]=1;
		F1[2][0]=0.5;
		F1[2][1]=1;
		F1[2][2]=-sin(X[2]);
		F1[2][3]=1;
		F1[3][0]=1;
		F1[3][1]=0.5;
		F1[3][2]=1;
		F1[3][3]=cos(X[3]);                /*对数组 F1 赋值*/
		for (k=0;k<=(n-2);k++)
			{ik=k;
			for(i1=k;i1<=(n-1);i1++)
			{if(fabs(F1[ik][k])<fabs(F1[i1][k]))    //选ik
				ik=i1;}
			for (j1=k;j1<=(n-1);j1++)                        //交换k行与ik行
				{c[k][j1]=F1[k][j1];
				F1[k][j1]=F1[ik][j1];
				F1[ik][j1]=c[k][j1];
				d[k]=F[k];
				F[k]=F[ik];
				F[ik]=d[k];
				}
			for (i1=k+1;i1<=(n-1);i1++)
				{m[i1]=F1[i1][k]/F1[k][k];
				for (j1=k+1;j1<=(n-1);j1++)
					{F1[i1][j1]=F1[i1][j1]-m[i1]*F1[k][j1];}       //计算a,b
				F[i1]=F[i1]-m[i1]*F[k];
				}
			}
			//回带过程
		xx[n-1]=F[n-1]/F1[n-1][n-1];
		for (k=(n-2);k>=0;k--)
			{sum1=0;
			for (j1=k+1;j1<=(n-1);j1++)
				{sum1=sum1+F1[k][j1]*xx[j1];}
				xx[k]=(F[k]-sum1)/F1[k][k];
			}
		err=max(xx,n)/max(X,n);
		if (err<=1e-12)
		{break;}
		for (i1=0;i1<n;i1++)
			{X[i1]=X[i1]+xx[i1];}
		flag=flag+1;			
	 }
		printf("(%10.6f,%10.6f):t=%10.6f,u=%10.6f\n",x,y,X[0],X[1]);
	}
	
getchar();
}