- 1、本文档共13页,可阅读全部内容。
- 2、原创力文档(book118)网站文档一经付费(服务费),不意味着购买了该文档的版权,仅供个人/单位学习、研究之用,不得用于商业用途,未经授权,严禁复制、发行、汇编、翻译或者网络传播等,侵权必究。
- 3、本站所有内容均由合作方或网友上传,本站不对文档的完整性、权威性及其观点立场正确性做任何保证或承诺!文档内容仅供研究参考,付费前请自行鉴别。如您付费,意味着您自己接受本站规则且自行承担风险,本站不退款、不进行额外附加服务;查看《如何避免下载的几个坑》。如果您已付费下载过本站文档,您可以点击 这里二次下载。
- 4、如文档侵犯商业秘密、侵犯著作权、侵犯人身权等,请点击“版权申诉”(推荐),也可以打举报电话:400-050-0827(电话支持时间:9:00-18:30)。
查看更多
拟牛顿法求解非线性方程组 投稿:冯蕏蕐
拟牛顿法求解非线性方程组 投稿:冯蕏蕐
//拟牛顿法求解非线性方程组
#include
#include
#include
void n_newton(n, eps, x, f)//拟牛顿法
int n;
double eps, x[];
void (*f)(double [],double []);
{
int i,j,k=0;
double am, z, beta, d;
double t=0.1, h=0.1;
double *a, *b, *y;
int gauss(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{
b[i] = y[i];
z=fabs(y[i]);
if (zam) am=z;
}
if (am=eps)
{
k = k + 1;
if (k 500)
{
printf(“迭代超过500次,可能不收敛!\n”);
free(a); free(b); free(y); return;
}
for (j=0; j{
z=x[j]; x[j]=x[j]+h;
(*f)(x,y);
for (i=0; ix[j]=z;
}
if (gauss(n,a,b) == 0)
{
printf(“系数矩阵奇异,Gauss消去法失败!\n”);
free(a); free(b); free(y); return;
}
beta=1.0;
for (i=0; iif (fabs(beta)+1.0==1.0)
{
printf(“Beta=0. 程序工作失败!\n”);
free(a); free(b); free(y); return;
}
d=h/beta;
for (i=0; ih=t*h;
}
}
printf(“实际迭代次数为: %d\n”,k);
free(a); free(b); free(y);
return;
}
int gauss(int n, double a[], double b[])//全选主元Gauss消去法
{
int *js, k, i, j, is;
double d, t;
js = (int *)malloc(n*sizeof(int));
for (k=0; k{
d=0.0;
for (i=k;ifor (j=k;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{
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{
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;jb[k]=b[k]/d;
for (i=k+1;i{
for (j=k+1;jb[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;jb[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 ()
{
int k;
double x[3] = {1.0, 1.0, 1.0};
void f(double [],double []);
n_newton(3, 0.0000001, x, f);
for (k=0; kreturn;
}
void f(double x[], double y[]) //计算方程组各左端函数值
{
y[0] = x[0]*x[0] + x[1]*x[1
文档评论(0)