应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 专业学科区 » 信息科学 » 程序算法 » Levenberg-Marquardt 算法的例子代码
51/1返回列表
查看: 1662  |  回复: 3
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

wsbl

金虫 (小有名气)

[交流] 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();
        }  
}
回复此楼

» 收录本帖的淘帖专辑推荐

软件学习

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼 2011-10-14 19:26:34
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lrs2008

木虫 (正式写手)

观天
★
小木虫(金币+0.5):给个红包，谢谢回帖
乖乖   这么长 ...
回复此楼
残
3楼2011-10-16 09:21:30
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 4 个回答

lingshuning

木虫 (职业作家)
★
小木虫(金币+0.5):给个红包，谢谢回帖
很完整，要是再有个程序注释就完美了，不过还是多谢楼主的分享！
赞一下(1人) 回复此楼
不管未来如何，我必竭尽所能
4楼2011-12-03 19:18:50
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 4 个回答
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 317求调剂 +3 申子申申 2026-03-19 5/250 2026-03-19 13:47 by houyaoxu
[考研] 一志愿福大288有机化学，求调剂 +3 小木虫200408204 2026-03-18 3/150 2026-03-19 13:31 by houyaoxu
[考研] 一志愿西安交通大学材料工程专业 282分求调剂 +5 枫桥ZL 2026-03-18 6/300 2026-03-19 13:24 by 枫桥ZL
[考研] 一志愿武理材料305分求调剂 +5 想上岸的鲤鱼 2026-03-18 6/300 2026-03-18 17:53 by 无际的草原
[考研] 311求调剂 +6 26研0 2026-03-15 6/300 2026-03-18 14:43 by haxia
[考研] 材料与化工一志愿南昌大学327求调剂推荐 +8 Ncdx123456 2026-03-13 9/450 2026-03-18 14:40 by haxia
[考研] 307求调剂 +3 冷笙123 2026-03-17 3/150 2026-03-18 09:55 by macy2011
[考研] 环境工程调剂 +8 大可digkids 2026-03-16 8/400 2026-03-18 09:36 by zhukairuo
[考研] 265求调剂 +3 梁梁校校 2026-03-17 3/150 2026-03-18 09:12 by zhukairuo
[考研] 材料与化工专硕调剂 +5 heming3743 2026-03-16 5/250 2026-03-17 14:03 by 勇敢太监王公公
[考研] 东南大学364求调剂 +5 JasonYuiui 2026-03-15 5/250 2026-03-16 21:28 by 木瓜膏
[考研] 机械专硕325，寻找调剂院校 +3 y9999 2026-03-15 5/250 2026-03-16 19:58 by y9999
[基金申请] 今年的国基金是打分制吗？ 50+3 zhanghaozhu 2026-03-14 3/150 2026-03-16 17:07 by 北京莱茵润色
[考研] 318求调剂 +3 Yanyali 2026-03-15 3/150 2026-03-16 16:41 by houyaoxu
[考研] 070300化学学硕求调剂 +6 太想进步了0608 2026-03-16 6/300 2026-03-16 16:13 by kykm678
[考研] 285求调剂 +6 ytter 2026-03-12 6/300 2026-03-16 15:05 by njzyff
[考研] 求老师收留调剂 +4 jiang姜66 2026-03-14 5/250 2026-03-15 20:11 by Winj1e
[考研] 材料工程调剂 +9 咪咪空空 2026-03-12 9/450 2026-03-13 22:05 by 星空星月
[考研] 085600材料与化工 309分请求调剂 +7 dtdxzxx 2026-03-12 8/400 2026-03-13 14:43 by jxchenghu
[考研] 321求调剂（食品/专硕） +3 mxcz321 2026-03-12 6/300 2026-03-13 08:45 by xc321
信息提示
关闭
请填处理意见
关闭