前面的函数编码已经调试过均没有问题,就是每次的运行结果不一样,有时候对有时候错,错的时候多,调整*alamda的大小运行结果也不一样,初步认为是函数传递的问题,那位高人指点下啊,多谢了!!
当*lamda=1(赋值)时
第一次运行结果:
n=41,total=0.000000
0.83857
2.35192
2.87649
2.22281
0.98173
2.03835
1.147371e-08
Null pointer assignment
第二次运行结果:
n=41,total=0.000000
0.83857
2.35192
2.87649
2.22281
0.98173
2.03835
1.147371e-08
第三次运行结果:
502.29750
-36.82217
6.15297
0.98418
3.82171
4.53858
9.321532e+00
第四次以后:
Null pointer assignment
Floating point error: Domain.
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <stddef.h>
#include <malloc.h>
double *vector(long nl, long nh)
void free_vector(double *v, long nl, long nh)
double **matrix(long nrl, long nrh, long ncl, long nch)
void free_matrix(double **m, long nrl, long nrh, long ncl, long nch)
int *ivector(long nl, long nh)
void gaussj(double **a, int n, double **b, int m)
void covsrt(double **covar, int ma, int ia[], int mfit)
void funcs(double x, double a[], double *y, double dyda[], int na)
void mrqcof(double x[], double y[], double sig[], int ndata,double a[], int ia[], int ma, double **alpha, double beta[], double *chisq)
main()
{
int i,j,k,l;
long int n;
static double x[7]={1,2,3,4,5,6,7};
static double y[7]={2.895,2.558,1.631,0.852,0.405,0.173,0.062};
static double sig[7]={1,1,1,1,1,1,1};
int ndata=7,ma=6;
double a[6]={2,3,4,3,2,3},a2[6];
double total=0;
int ia[6]={1,1,1,1,1,1};
double *alamda, *chisq;
double **alpha,tr;
double **covar;
static int mfit;
static double ochisq,*atry,*da,**oneda,*beta;
double *p1[6],*p2[6];
*alamda=-0.1;
*chisq=10000;
covar=matrix(1,ma,1,ma);
alpha=matrix(1,ma,1,ma);
for(j=0;j<ma;j++)
for(k=0;k<ma;k++){ covar[j][k]=0.0;alpha[j][k]=0;}
for(n=0;n<100;n++)
{ total=0;
for(i=0;i<ma;i++) a2[i]=a[i];
if(*alamda<0.0) {
atry=vector(1,ma);
beta=vector(1,ma);
da=vector(1,ma);
for(mfit=0,j=0;j<ma;j++)
if(ia[j]) mfit++;
oneda=matrix(1,mfit,1,1);
*alamda=1;
mrqcof(x,y,sig,ndata,a,ia,ma,alpha,beta,chisq);
ochisq=(*chisq);
for(j=0;j<ma;j++) atry[j]=a[j];
}
for(j=0;j< mfit;j++)
for(k=0;k< mfit;k++) covar[j][k]=alpha[j][k];
for(j=0;j< mfit;j++) {
covar[j][j]=alpha[j][j]*(1.0+(*alamda));
oneda[j][0]=beta[j];
}
p1[0]=&covar[0][0];
p2[0]=&oneda[0][0];
for(i=1;i<6;i++)
{ p1[i]=covar[i];
p2[i]=oneda[i];
}
gaussj(p1,mfit,p2,1);
covar[0]=&p1[0][0];
oneda[0]=&p2[0][0];
for(i=1;i<6;i++)
{ covar[i]=p1[i];
oneda[i]=p2[i];
}
for(j=0;j<mfit;j++) da[j]=oneda[j][0];
for(j=-1,l=0;l<ma;l++)
if(ia[l]) atry[l]=a[l]+da[++j];
mrqcof(x,y,sig,ndata,atry,ia,ma,covar,da,chisq);
if(*chisq<ochisq) {
*alamda*=0.1;
ochisq=(*chisq);
for(j=0;j<mfit;j++) {
for(k=0;k<mfit;k++) alpha[j][k]=covar[j][k];
beta[j]=da[j];
}
for(j=0;j<ma;j++) a[j]=atry[j];
for(j=0;j<ma;j++) {tr=fabs(a2[j]-a[j]);total=total+tr;}
fprintf(fout,"total=%f\n",total);
if(total<1e-6)
{printf("n=%d,total=%f\n",n,total); break;}
} else {
*alamda*=10.0;
*chisq=ochisq;
}
}
for(i=0;i<ma;i++)
printf("%9.5f\n",a[i]);
printf("%13.7e\n",*chisq);
if(*alamda==0) {
covsrt(covar,ma,ia,mfit);
covsrt(alpha,ma,ia,mfit);
}
free_matrix(alpha,1,ma,1,ma);
free_matrix(covar,1,ma,1,ma);
free_matrix(oneda,1,mfit,1,1);
free_vector(da,1,ma);
free_vector(beta,1,ma);
free_vector(atry,1,ma);
}