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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 专业学科区 » 信息科学 » 程序算法 » Levenberg-Marquardt 算法的例子代码
41/1返回列表
查看: 1658  |  回复: 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返回列表
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 326求调剂 +3 上岸的小葡 2026-03-15 4/200 2026-03-15 18:50 by 无际的草原
[考研] 化学工程321分求调剂 +5 大米饭！ 2026-03-15 5/250 2026-03-15 18:49 by a不易
[考研] 321求调剂 +3 大米饭！ 2026-03-15 3/150 2026-03-15 17:48 by 哈哈哈哈嘿嘿嘿
[考研] 材料与化工（0856）304求B区调剂 +7 邱gl 2026-03-10 11/550 2026-03-14 12:18 by 邱gl
[考研] 332分材料工程调剂 +3 莓好时光海苔 2026-03-09 3/150 2026-03-14 02:03 by JourneyLucky
[考研] 求调剂 +6 yfihxh 2026-03-09 6/300 2026-03-14 01:18 by JourneyLucky
[考研] b区环境工程求调剂 +4 Maps1 2026-03-10 6/300 2026-03-14 00:23 by JourneyLucky
[考研] 四川大学085601材料工程专硕 初试294求调剂 +4 祝我们好在冬天 2026-03-11 4/200 2026-03-13 21:39 by peike
[考研] 315求调剂 +9 小羊小羊_ 2026-03-11 10/500 2026-03-13 21:13 by SXNU李老师
[硕博家园] 深圳大学硕士招生（2026秋，传感器方向，仅录取第一志愿） +4 xujiaoszu 2026-03-11 7/350 2026-03-13 17:28 by xujiaoszu
[考研] 求b区学校调剂 +3 周56 2026-03-11 3/150 2026-03-13 16:20 by JourneyLucky
[考研] 310求调剂 +3 【上上签】 2026-03-11 3/150 2026-03-13 16:16 by JourneyLucky
[考研] 工科278分求调剂 +5 周慢热啊 2026-03-12 7/350 2026-03-13 15:49 by JourneyLucky
[考研] 314求调剂 +7 无懈可击的巨人 2026-03-12 7/350 2026-03-13 15:40 by JourneyLucky
[考研] 274求调剂 +3 S.H1 2026-03-12 3/150 2026-03-13 15:15 by JourneyLucky
[考研] 285求调剂 +4 ytter 2026-03-12 4/200 2026-03-13 14:48 by jxchenghu
[考研] 321求调剂（食品/专硕） +3 xc321 2026-03-12 6/300 2026-03-13 08:45 by xc321
[考研] 290求调剂 +3 柯淮然 2026-03-10 8/400 2026-03-11 13:48 by 柯淮然
[考研] 化工0817调剂 +8 灿若星晨 2026-03-10 8/400 2026-03-10 22:44 by 星空星月
[考研] 085602化工求调剂 +7 董boxing 2026-03-10 7/350 2026-03-10 17:07 by BruceLiu320
信息提示
关闭
请填处理意见
关闭