24小时热门版块排行榜

>论坛更新日志 (4346)
>文献求助 (471)
>虫友互识 (402)
>导师招生 (377)
>考博 (206)
>博后之家 (183)
>硕博家园 (163)
>休闲灌水 (162)
>论文投稿 (136)
>基金申请 (96)
>考研 (74)
>招聘信息布告栏 (70)
>教师之家 (60)
>绿色求助(高悬赏) (47)
>公派出国 (46)
>找工作 (42)

返回列表

当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

wsbl

金虫 (小有名气)

应助: 1 (幼儿园)
金币: 1845.9
散金: 50
红花: 1
帖子: 114
在线: 36.4小时
虫号: 1147591
注册: 2010-11-15
性别: GG
专业: 无机材料化学

[交流] Levenberg-Marquardt 算法的例子代码已有3人参与

#include "stdio.h"
#include
#include
void out(double*p,int n,int m)//输出矩阵
{int i,j;

for (i=0;i    for (j=0;j {
   if ((m*i+j)%m==0)
   printf("\n"

;
printf("%16.8f",*(p+m*i+j));

}
printf("\n\n"

;
}

void LM(double *c, double *sigma, int N, double sigma_H2O,double r[3])
{
int i,j,maxtime=100; //如果不收敛，允许的最大迭代次数
double s,snext,det,temp0;


double x[]={.5,.5},xnext[2],ep[2];

double f[6];
double f_x1[6];
double f_x2[6];

double m[4];  //矩阵

double temp1[4];//逆矩阵
double temp2[2];
   double temp3[2];

double alph,beita;

j=0;
   for (i=0;i {f=(sigma-(sigma_H2O-sigma_H2O*x[0]*log(1+c/x[1])));
   f_x1=sigma_H2O*log(1+c/x[1]);
   f_x2=-sigma_H2O*x[0]*c/(x[1]*x[1])/(1+c/x[1]);
}

  s=.0;
  for (i=0;i s=s+f*f;

  temp0=.0;
  for (i=0;i   temp0=temp0+f_x1*f_x1;
  m[0]=temp0;

  temp0=.0;
  for (i=0;i   temp0=temp0+f_x1*f_x2;
  m[1]=temp0;
  m[2]=temp0;

  temp0=.0;
  for (i=0;i   temp0=temp0+f_x2*f_x2;
  m[3]=temp0;

//out(m,2,2);//到此步正确。
  alph=(fabs(*m)+fabs(*(m+3)))/2;
  //printf("%8.5e\n",alph);
  beita=2.0;

m[0]=m[0]+alph;
m[3]=m[3]+alph;
//out(m,2,2);//到此步正确。

  det=1.0/(m[0]*m[3]-m[1]*m[2]);
//printf("迭%15.10f\n",det);
  temp1[0]=det*m[3];
  temp1[1]=-det*m[1];
  temp1[2]=-det*m[1];
  temp1[3]=det*m[0];

//out(temp1,2,2);//到此步正确。
//temp2=mult(temp,f,2,N,1);
  temp0=.0;
  for (i=0;i   temp0=temp0+f*f_x1;
  temp2[0]=temp0;

  temp0=.0;
  for (i=0;i   temp0=temp0+f*f_x2;
  temp2[1]=temp0;

//out(temp2,2,1);//到此步正确。
//temp3=mult(temp1,temp2,2,2,1);
  temp0=.0;
  for (i=0;i<2;i++)
{ temp0=temp0+temp1*temp2;
}
temp3[0]=temp0;

  temp0=.0;
  for (i=0;i<2;i++)
{ temp0=temp0+temp1[2+i]*temp2;
}
temp3[1]=temp0;

//out(temp3,2,1);//到此步正确。

  xnext[0]=x[0]-temp3[0];
  xnext[1]=x[1]-temp3[1];
//out(xnext,2,1);//到此步正确。

  snext=0;
  for (i=0;i   { f=(sigma-(sigma_H2O-sigma_H2O*xnext[0]*log(1+c/xnext[1])));
   snext=snext+f*f;
  }
//printf("%8.5e",snext);//到此步正确。

  ep[0]=xnext[0]-x[0];
  ep[1]=xnext[1]-x[1];

//out(ep,2,1);//到此步正确。

while (((fabs(ep[0])+fabs(ep[1]))/2)>1.0e-6 && j {j=j+1;
if (snext>=s)
   alph=alph*beita;
else
   alph=alph/beita;
   //printf("%8.3f\n",(float)j);

   x[0]=xnext[0];
   x[1]=xnext[1];
   s=snext;

   for (i=0;i {f=(sigma-(sigma_H2O-sigma_H2O*x[0]*log(1+c/x[1])));
   f_x1=sigma_H2O*log(1+c/x[1]);
   f_x2=-sigma_H2O*x[0]*c/(x[1]*x[1])/(1+c/x[1]);
}

s=.0;
for (i=0;i s=s+f*f;

temp0=.0;
  for (i=0;i   temp0=temp0+f_x1*f_x1;
  m[0]=temp0;

  temp0=.0;
  for (i=0;i   temp0=temp0+f_x1*f_x2;
  m[1]=temp0;
  m[2]=temp0;

  temp0=.0;
  for (i=0;i   temp0=temp0+f_x2*f_x2;
  m[3]=temp0;

m[0]=m[0]+alph;
m[3]=m[3]+alph;
//out(m,2,2);//到此步正确。

  det=1.0/(m[0]*m[3]-m[1]*m[2]);
//printf("迭%15.10f\n",det);
  temp1[0]=det*m[3];
  temp1[1]=-det*m[1];
  temp1[2]=-det*m[1];
  temp1[3]=det*m[0];

//out(temp1,2,2);//到此步正确。
//temp2=mult(temp,f,2,N,1);
temp0=.0;
  for (i=0;i   temp0=temp0+f*f_x1;
  temp2[0]=temp0;

  temp0=.0;
  for (i=0;i   temp0=temp0+f*f_x2;
  temp2[1]=temp0;

//out(temp2,2,1);//到此步正确。
//temp3=mult(temp1,temp2,2,2,1);
  temp0=.0;
  for (i=0;i<2;i++)
{ temp0=temp0+temp1*temp2;
}
temp3[0]=temp0;

  temp0=.0;
  for (i=0;i<2;i++)
{ temp0=temp0+temp1[2+i]*temp2;
}
temp3[1]=temp0;

//out(temp3,2,1);//到此步正确。

xnext[0]=x[0]-temp3[0];
xnext[1]=x[1]-temp3[1];
//out(xnext,2,1);//到此步正确。

snext=0;
for (i=0;i { f=(sigma-(sigma_H2O-sigma_H2O*xnext[0]*log(1+c/xnext[1])));
   snext=snext+f*f;
}

ep[0]=xnext[0]-x[0];
ep[1]=xnext[1]-x[1];

}
r[0]=xnext[0];
r[1]=xnext[1];
r[2]=(double)j;
}

void main()
{
//double p[] = {0.02,0.04,0.08,0.12,0.16,0.2,0.24,0.28};//mol/L
//double q[] = {0.06810,0.06442,0.0601,0.05678,0.05467,0.05237,0.05108,0.05025};//N/m

double p[] = {0.02,0.04,0.08,0.12,0.16,0.2};
double q[] = {0.06810,0.0644,0.0601,0.0568,0.0547,0.0524};
double r[3];
LM(p,q,6,0.07197,r);

printf("迭代次数：%5.0f\n",r[2]);

if (r[2]<100)
{
printf("\n希什科夫模型sigma_H2O-sigma_H2O*a*ln(1+x/b)中\n\n参数a的最优值为:%9.9f\n",r[0]);
printf("\n参数b的最优值为:%9.9f\n\n",r[1]);
}
  else
printf("\n用Levenberg法失败\n\n";

  for(int i=0;i<1;i++)//延时几秒。
{
printf("按回车键结束程序:";
getchar();
}
}