常微分方程组的参数估值（整整弄了一个月没解决！望高手帮帮忙）

[znpmg]

先用Treanor法解常微分方程的数值解，然后用DFP法（求解无约束n元函数的极值问题）寻优。

#include <stdio.h>
#include <math.h>
#include "stdlib.h"

#define MAX_N 10
#define MAX_ITERATE 3000
#define DEBUG 1

double y[2];
double result[2][8];

double DFP(int N,double XK[],double EPS,double(*f)(int ,double[]));
double Parab(int N,double X0[],double P[],double H,double EPS,double (*f)(int ,double[] ),double* TM,double X[]);
void gtnr2(double t,double h,int n,double y[],int k,double x[],double z[]);
/*----------------测试函数-----------------------------*/
double OBF_For_DFP(int N,double X[]);
void testTreanor();
double AimFunction(double z[2][8]);
double OBF_For_Treanor(int N,double x[]);


/*----For testing Parab-----*/
double OBF_For_Parab(int N,double X[])
{
double FX=0.0,X1;
int i;
for(i=0;i<N;i++)
{
X1=X[i]-(i+1);
FX=FX+X1*X1;
}
return FX;
}
void testParab()
{
int N;
double EPS,FX,H,TM;
double X0[MAX_N],P[MAX_N],X[MAX_N];
int i;
N=5;
for(i=0;i<N;i++)
{
X0[i]=1;
P[i]=(i+1)/10.0;
}
H=1;
EPS=1.0E-7;
TM=7;
FX=Parab(N,X0,P,H,EPS,OBF_For_Parab,&TM,X);

printf("\n\nThe test result of Parab as following:\n");
printf("TM=%f,FX=%f\n",TM,FX);
for(i=0;i<N;i++)
printf("X[%d]=%f\t",i+1,X[i]);
}
/*----For testing DFP-----*/
double OBF_For_DFP(int N,double X[])
{
double f1,f2,f3;
/*f1= X[0]-8;
f2= X[1]-1;//+5*sin(X[0]);
f3=X[2]-cos(X[2]);
//f3=1.0/2-cos(X[2]);
return f1*f1+f2*f2+f3*f3;
*/
f1=X[1]-X[0]*X[0];
f2=1-X[0];
return 100*f1*f1+f2*f2;
}
void testDFP()
{
int N;
double EPS,FX;
double XK[MAX_N];
int i;

N=2;
EPS=1.0E-10;
//for(i=0;i<N;i++)
// XK[i]=1;
XK[0]=-1.2;
XK[1]=1.5;
FX=DFP(N,XK,EPS,OBF_For_DFP);

printf("The test result of DFP as following:\n");
for(i=0;i<N;i++)
printf("XK[%d]=%1.30f\t",i+1,XK[i]);
printf("\nFX=%1.30f\n",FX);
}
/*------For test Treanor------------------*/
void testTreanor()
{
int i,j;
double x[8];
double f;
y[0]=0;
y[1]=0;

x[0]=0.1151;
x[1]=2.45;
x[2]=0.02;
x[3]=0.3047;
x[4]=2.75;
x[5]=0.008;
x[6]=0.85;
x[7]=0.51;
f=OBF_For_Treanor(8,x);
while(f>0.1)
{
printf("f=%f\n",f);
for(i=0;i<8;i++)
{
printf("result[0][%d]=%f\tresult[1][%d]=%f\n",i,result[0][i],i,result[1][i]);
}
DFP(8,x,1.0e-10,OBF_For_Treanor);
for(i=0;i<8;i++)
printf("x[%d]=%f\t",i,x[i]);
y[0]=0;
y[1]=0;
f=OBF_For_Treanor(8,x);
}
for(i=0;i<8;i++)
printf("x[i]=%f",x[i]);
}
double AimFunction(double z[2][8])
{
double A[2][8];
double f;
int i,j;
double x[8];
f=0.0;
A[0][0]=0.074;
A[0][1]=0.323;
A[0][2]=0.295;
A[0][3]=0.319;
A[0][4]=0.315;
A[0][5]=0.311;
A[0][6]=0.289;
A[0][7]=0.294;
A[1][0]=0.515;
A[1][1]=0.592;
A[1][2]=0.684;
A[1][3]=0.667;
A[1][4]=0.704;
A[1][5]=0.723;
A[1][6]=0.765;
A[1][7]=0.773;
for(i=0;i<8;i++)
{
f+=(z[0][i]-A[0][i])*(z[0][i]-A[0][i]);
f+=(z[1][i]-A[1][i])*(z[1][i]-A[1][i]);
}
return f;
}
double OBF_For_Treanor(int N,double x[])
{
int i;

double t,h,z[2][101];
double f;
t=0.0;
h=2;

gtnr2(t,h,2,y,101,x,z);
result[0][0]=z[0][0];
result[1][0]=z[1][0];
for(i=0;i<7;i++)
{
result[0][i+1]=z[0][i*15+5];
result[1][i+1]=z[1][i*15+5];
}
/*for(i=0;i<8;i++)
{
printf("result[0][%d]=%f\tresult[1][%d]=%f",i,result[0][i],i,result[1][i]);
}*/
return AimFunction(result);
}
void gtnr2f(t,y,n,x,d)
int n;
double t,y[],d[],x[];
{
d[0]=x[0]*pow((0.8947-y[0]-y[1]),x[1])*pow((0.726-y[0]-y[1]),x[2])-x[3]*pow(y[0],x[4]);
d[1]=x[5]*pow((0.8947-y[0]-y[1]),x[6])*pow((0.726-y[0]-y[1]),x[7]);
return;
}

int main()
{
// testParab(); /*用例子中的函数测试Parab函数*/
// testDFP(); /*用例子中的函数测试DFP函数*/
testTreanor();
// getch();

return 0;
}

[znpmg]

/*---------------------------------------------------------------*/
double DFP(int N,double XK[],double EPS,double(*f)(int ,double[]))
{
double P[MAX_N],G[MAX_N],Y[MAX_N],H[MAX_N][MAX_N],X[MAX_N];
double FX,TM,E,D1,D2,S;
int i,j,k;
long nTimes;

S=1.0E-5;

for(i=0;i<N;i++)
X[i]=XK[i];
/*最大循环次数MAX_ITERATE，若迭代MAX_ITERATE次后偏差还是大于EPS，则输出结果*/
for(nTimes=0;nTimes<MAX_ITERATE;nTimes++)
{
/*开始矩阵H[N,N]=I*/
for(i=0;i<N;i++)
{
for(j=0;j<N;j++)
H[i][j]=0;
H[i][i]=1;
}

/*用S为步长的中心差分计算F(X)的梯度，存放在G中*/
for(i=0;i<N;i++)
{
X[i]=XK[i]+S;
G[i]=f(N,X);
X[i]=XK[i]-S;
G[i]=0.5*(G[i]-f(N,X))/S;
X[i]=XK[i];
}
/*变尺度法一般每做N次共轭方向后再重新开始*/
for(i=0;i<N;i++)
{
/*利用公式P=-H*G求出搜索方向P*/
for(j=0;j<N;j++)
{
P[j]=0;
for(k=0;k<N;k++)
P[j]-=H[j][k]*G[k];
}
TM=1;
FX=Parab(N,XK,P,0.2,0.0001,f,&TM,X);
for(j=0;j<N;j++)
{
P[j]*=TM;
XK[j]+=P[j];
}
E=0;
for(j=0;j<N;j++)
{
X[j]=XK[j]+S;
Y[j]=f(N,X);
X[j]=XK[j]-S;
FX=f(N,X);
Y[j]=0.5*(Y[j]-FX)/S;
E=E+Y[j]*Y[j];
Y[j]-=G[j];
X[j]=XK[j];
}
/*如果梯度的模小于EPS，最有点就是XK,并返回最小值FX*/
if(E<EPS)
{
if(DEBUG)
printf("nTimes=%d\n",nTimes);
return FX;
}
D1=D2=0;
for(j=0;j<N;j++)
{
G[j]=0;
for(k=0;k<N;k++)
G[j]+=H[j][k]*Y[j];
D1+=P[j]*Y[j];
D2+=G[j]*Y[j];
}
for(j=0;j<N;j++)
for(k=0;k<N;k++)
H[j][k]+=P[j]*P[k]/D1-G[j]*G[k]/D2;
for(j=0;j<N;j++)
G[j]=Y[j];
}
}
return FX;
}

double Parab(int N,double X0[],double P[],double H,double EPS,double (*f)(int ,double[] ),double* TM,double X[])
{
double S,FM=0,DD,DU,T[3],Z[3],Y0[3],Y[3];
int i,M;

S=1.0E-5;

M=0;
Again: for(i=0;i<N;i++)
X[i]=X0[i]+*TM*P[i];
Y0[0]=f(N,X);
for(i=0;i<N;i++)
X[i]=X0[i]+(*TM+S)*P[i];
Y0[1]=f(N,X);
if(Y0[1]-Y0[0]<0)
{
T[0]=*TM;
Y[0]=Y0[0];
T[1]=*TM+S;
Y[1]=Y0[1];
while(1)
{
T[2]=T[1]+H;
for(i=0;i<N;i++)
X[i]=X0[i]+T[2]*P[i];
Y0[2]=f(N,X);
if(Y0[2]<Y[1])
{
T[0]=T[1];
Y[0]=Y[1];
T[1]=T[2];
Y[1]=Y0[2];
}
else
{
Y[2]=Y0[2];
break;
}
}
}
else if(Y0[1]-Y0[0]>0)
{
T[2]=*TM+S;
Y[2]=Y0[1];
T[1]=*TM;
Y[1]=Y0[0];
while(1)
{
T[0]=T[1]-H;
for(i=1;i<N;i++)
X[i]=X0[i]+T[0]*P[i];
Y0[2]=f(N,X);
if(Y0[2]<Y[1])
{
T[2]=T[1];
Y[2]=Y[1];
T[1]=T[0];
Y[1]=Y0[2];
}
else
{
Y[0]=Y0[2];
break;
}
}
}
else
{
*TM=*TM+S;
M=M+1;
if(M>=1000) /*如果相等的次数达到1000次，说明f可能是一个常函数*/
return Y0[1];
else
goto Again;
}
while(1)
{
Z[0]=Y[1]-Y[2];
Z[1]=Y[2]-Y[0];
Z[2]=Y[0]-Y[1];
DU=DD=0;
for(i=0;i<3;i++)
{
DU+=T[i]*T[i]*Z[i];
DD+=T[i]*Z[i];
}
if(DD+DU-DU==0.0)
{
*TM=T[1];
return FM;
}
else
*TM=0.5*DU/DD;
for(i=0;i<N;i++)
X[i]=X0[i]+(*TM)*P[i];
FM=f(N,X);
if(fabs(*TM-T[1])<=EPS)
return FM;
if(fabs(*TM-T[1])<0.001&&fabs(FM-Y[1])<=EPS)
return FM;

if(*TM<T[1]) /**TM在区间(T[0],T[1])之间*/
{
if(FM < Y[1])
{
T[2]=T[1];
Y[2]=Y[1];
T[1]=*TM;
Y[1]=FM;
continue;
}
else
{
T[0]=*TM;
Y[0]=FM;
continue;
}
}
else
{
if(FM<Y[1])
{
T[0]=T[1];
Y[0]=Y[1];
T[1]=*TM;
Y[1]=FM;
continue;
}
else
{
T[2]=*TM;
Y[2]=FM;
continue;
}
}
return FM;
}
return FM;
}

void gtnr2(double t,double h,int n,double y[],int k,double x[],double z[])
{
int i,j;
double a,s,aa,bb,dd,g,dy,dy1,pm,*d,*p,*w,*q,*r,*yn;
d=malloc(n*sizeof(double));
p=malloc(n*sizeof(double));
w=malloc(4*n*sizeof(double));
q=malloc(4*n*sizeof(double));
r=malloc(4*n*sizeof(double));
yn=malloc(n*sizeof(double));
a=t;
for (j=0; j<=n-1; j++)
{
yn[j]=y[j];
z[j*k]=yn[j];
}
for (i=1; i<=k-1; i++)
{
t=a+(i-1)*h;
for (j=0; j<=n-1; j++)
w[j]=yn[j];
gtnr2f(t,yn,n,x,d);
for (j=0; j<=n-1; j++)
{
q[j]=d[j]; yn[j]=w[j]+h*d[j]/2.0;
w[n+j]=yn[j];
}
s=t+h/2.0;
gtnr2f(s,yn,n,x,d);
for (j=0; j<=n-1; j++)
{
q[n+j]=d[j];
yn[j]=w[j]+h*d[j]/2.0;
w[n+n+j]=yn[j];
}
gtnr2f(s,yn,n,x,d);
for (j=0; j<=n-1; j++)
q[n+n+j]=d[j];
pm=0.0;
for (j=0; j<=n-1; j++)
{
aa=q[n+n+j]-q[n+j];
bb=w[n+n+j]-w[n+j];
if(pm<-aa/bb)
pm=-aa/bb;
}
for (j=0; j<=n-1; j++)
{
if(pm>0.0)
{
dd=-pm*h;
r[j]=exp(dd);
r[n+j]=(r[j]-1.0)/dd;
r[n+n+j]=(r[n+j]-1.0)/dd;
r[3*n+j]=(r[n+n+j]-1.0)/dd;
}
else pm=0.0;
if (p[j]<=0.0)
g=q[n+n+j];
else
{
g=2.0*(q[n+n+j]-q[j])*r[n+n+j];
g=g+(q[j]-q[n+j])*r[n+j]+q[n+j];
}
w[3*n+j]=w[j]+g*h;
yn[j]=w[3*n+j];
}
s=t+h;
gtnr2f(s,yn,n,x,d);
for (j=0; j<=n-1; j++) q[3*n+j]=d[j];
for (j=0; j<=n-1; j++)
{
if (p[j]<=0.0)
{
dy=q[j]+2.0*(q[n+j]+q[n+n+j]);
dy=(dy+q[n+n+n+j])*h/6.0;
}
else
{
dy=-3.0*(q[j]+p[j]*w[j])+2.0*(q[n+j]
+p[j]*w[n+j]);
dy=dy+2.0*(q[n+n+j]+p[j]*w[n+n+j]);
dy=dy-(q[n+n+n+j]+p[j]*w[n+n+n+j]);
dy=dy*r[n+n+j]+q[j]*r[n+j];
dy1=q[j]-q[n+j]-q[n+n+j]+q[n+n+n+j];
dy1=dy1+(w[j]-w[n+j]-w[n+n+j]+w[n+n+n+j])*p[j];
dy=(dy+4.0*dy1*r[n+n+n+j])*h;
}
yn[j]=w[j]+dy;
z[j*k+i]=yn[j];
}
}
free(d); free(p); free(w); free(q); free(r);free(yn);
return;
}

[znpmg]

运算到第2步结果就溢出了，不知道是不是DFP算法部分写的有问题。

[dongpy]

不懂,学习,帮你顶一下

[bill_li]

不懂,学习,帮你顶一下

[江南书童]

我建议你把算法仔细想一下!猜想出最可能出错的地方,再设置断点,进行单步执行查看你的每个变量的变化有没有错的!

[nwpulipeng]

高手啊！！！

[greenteanet]

偶在网吧，没有编译器，不能帮你拉。。

[idler]

up一下，楼主以后最好加点注释。。。这个。。。数值算法最怕这个。。。

[laomai]

先mark,有空研究,呵呵

[daylove]

只能看看，友情up一下！

[melonliu]

有点晕……
虽然扔掉高数了，还是要copy下来研究研究