24小时热门版块排行榜    

北京石油化工学院2026年研究生招生接收调剂公告
查看: 1667  |  回复: 3

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();
        }  
}
回复此楼
已阅   回复此楼   关注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):给个红包,谢谢回帖
很完整,要是再有个程序注释就完美了,不过还是多谢楼主的分享!
不管未来如何,我必竭尽所能
4楼2011-12-03 19:18:50
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 wsbl 的主题更新
普通表情 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 合肥区域性重点一本招收调剂 +4 6266jl 2026-03-30 4/200 2026-03-30 21:48 by zhuangyan123
[考研] 277跪求调剂 +8 1915668 2026-03-27 12/600 2026-03-30 21:01 by dophin1985
[考研] 11408总分309,一志愿东南大学求调剂,不挑专业 +5 天赋带到THU 2026-03-29 6/300 2026-03-30 20:49 by dick_runner
[考研] 285求调剂 +4 AZMK 2026-03-30 7/350 2026-03-30 20:24 by laoshidan
[考研] 08工科,295,接受跨专业调剂 +3 lmnlzy 2026-03-30 3/150 2026-03-30 17:49 by wangjy2002
[考研] 322求调剂:一志愿湖南大学 材料与化工(085600),已过六级。 +9 XX小邓 2026-03-29 9/450 2026-03-30 17:18 by limeifeng
[考研] 342求调剂 +4 加油a李zs 2026-03-26 4/200 2026-03-30 16:39 by 晶体之美
[考研] 085701求调剂初试286分 +5 secret0328 2026-03-28 5/250 2026-03-30 12:54 by fangnagu
[考研] 求调剂 +10 张zz111 2026-03-27 11/550 2026-03-30 09:17 by 无际的草原
[考研] 085600,材料与化工321分求调剂 +10 大馋小子 2026-03-28 10/500 2026-03-29 23:35 by 飞行日记西
[考研] 2026年华南师范大学欢迎化学,化工,生物,生医工等专业优秀学子加入! +3 llss0711 2026-03-28 6/300 2026-03-29 10:26 by llss0711
[考研] 356求调剂 +3 gysy?s?a 2026-03-28 3/150 2026-03-29 00:33 by 544594351
[考研] 315求调剂 +4 akie... 2026-03-28 5/250 2026-03-28 21:05 by zhq0425
[考研] 本科新能源科学与工程,一志愿华理能动285求调剂 +3 AZMK 2026-03-27 5/250 2026-03-28 16:19 by xxxsssccc
[考研] 347求调剂 +3 山顶见α 2026-03-25 3/150 2026-03-28 14:13 by 唐沐儿
[考研] 材料与化工(0856)304求B区调剂 +8 邱gl 2026-03-27 8/400 2026-03-28 12:42 by 唐沐儿
[考研] 086000调剂 +3 7901117076 2026-03-26 3/150 2026-03-27 21:34 by Jianing_Mi
[考研] 一志愿郑大085600,310分求调剂 +5 李潇可 2026-03-26 5/250 2026-03-27 11:14 by 不吃魚的貓
[考研] 一志愿 南京邮电大学 288分 材料考研 求调剂 +3 jl0720 2026-03-26 3/150 2026-03-26 13:39 by zzll406
[考研] 求b区院校调剂 +4 周56 2026-03-24 5/250 2026-03-25 17:12 by yishunmin
信息提示
请填处理意见