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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 专业学科区 » 信息科学 » 程序算法 » Levenberg-Marquardt 算法的例子代码
41/1返回列表
查看: 1659  |  回复: 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的回帖

stockline

铁虫 (初入文坛)
不错
回复此楼
2楼2011-10-16 07:18:24
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lrs2008

木虫 (正式写手)

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

lingshuning

木虫 (职业作家)
★
小木虫(金币+0.5):给个红包，谢谢回帖
很完整，要是再有个程序注释就完美了，不过还是多谢楼主的分享！
赞一下(1人) 回复此楼
不管未来如何，我必竭尽所能
4楼2011-12-03 19:18:50
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 wsbl 的主题更新
41/1返回列表
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[硕博家园] 深圳大学硕士招生（2026秋，传感器方向，仅录取第一志愿） +4 xujiaoszu 2026-03-11 9/450 2026-03-17 10:29 by xujiaoszu
[考研] 本人考085602 化学工程 专硕 +15 不知道叫什么！ 2026-03-15 17/850 2026-03-17 10:13 by Sammy2
[考研] 333求调剂 +3 文思客 2026-03-16 7/350 2026-03-16 18:21 by 文思客
[考研] 312求调剂 +3 陌宸希 2026-03-16 4/200 2026-03-16 15:06 by peike
[考研] 297一志愿上交085600求调剂 +5 指尖八千里 2026-03-14 5/250 2026-03-14 17:26 by a不易
[考研] 311求调剂 +5 牛乳糖的卡卡 2026-03-10 5/250 2026-03-14 00:05 by JourneyLucky
[考研] 341求调剂 +4 番茄头--- 2026-03-10 4/200 2026-03-13 23:12 by JourneyLucky
[考研] 304求调剂 +6 Mochaaaa 2026-03-12 7/350 2026-03-13 22:18 by 星空星月
[考研] 求材料调剂 085600英一数二总分302 前三科235 精通机器学习 一志愿哈工大 +4 林yaxin 2026-03-12 4/200 2026-03-13 22:04 by 星空星月
[考研] 求调剂 +7 18880831720 2026-03-11 7/350 2026-03-13 16:10 by JourneyLucky
[考研] 求调剂 +3 程雨杭 2026-03-12 3/150 2026-03-13 15:06 by JourneyLucky
[考研] 308求调剂 +3 是Lupa啊 2026-03-12 3/150 2026-03-13 14:30 by 求调剂zz
[论文投稿] 投稿问题 5+4 星光灿烂xt 2026-03-12 6/300 2026-03-13 14:17 by god_tian
[考研] 0856化工原理 +6 z2839474511 2026-03-10 6/300 2026-03-13 10:41 by houyaoxu
[考研] 化工学硕306求调剂 +9 42838695 2026-03-12 9/450 2026-03-13 10:16 by houyaoxu
[考研] 270求调剂 085600材料与化工专硕 +3 YXCT 2026-03-11 3/150 2026-03-13 10:13 by houyaoxu
[考研] 321求调剂（食品/专硕） +3 xc321 2026-03-12 6/300 2026-03-13 08:45 by xc321
[考研] 333求调剂 +3 152697 2026-03-12 4/200 2026-03-13 07:08 by Iveryant
[基金申请] 提交后的基金本子，已让学校撤回了，可否换口子提交 +3 dut_pfx 2026-03-10 3/150 2026-03-11 08:38 by kudofaye
[考研] 收调剂 +7 调剂的考研学生 2026-03-10 7/350 2026-03-10 17:57 by 麦茶汤圆
信息提示
关闭
请填处理意见
关闭