📄 d9r4.cpp
字号:
#include "iostream.h"
#include "math.h"
#include "stdlib.h"
#include "string.h"
int sgn (double f)
{
if (f>0)
{
return 1;
}
if (f<0)
{
return -1;
}
return 0;
}
void svbksb(double u[], double w[], double v[], int m, int n, double b[], double x[])
{
int i,j,jj;
double tmp[101];
double s;
for (j = 1; j<=n; j++)
{
s = 0.0;
if (w[j] != 0.0)
{
for (i =1; i<=m; i++)
{
s = s + u[(i-1)*n+j] * b[i];
}
s = s / w[j];
}
tmp[j] = s;
}
for (j = 1; j<=n; j++)
{
s = 0.0;
for (jj = 1; jj<=n; jj++)
{
s = s + v[(j-1)*n+jj] * tmp[jj];
}
x[j] = s;
}
}
void svdcmp(double a[], int m, int n, double w[], double v[])
{
int i,j,k,l,its,nm,aaaaa,pp;
double rv1[100];
if (m < n)
{
cout<<"you must augment a with extra zero rows."<<endl;
}
double x,y,z,c,f,h,s,g = 0.0;
double scale1 = 0.0;
double anorm = 0.0;
for (i = 1; i<=n;i++)
{
l = i + 1;
rv1[i] = scale1 * g;
g = 0.0;
s = 0.0;
scale1 = 0.0;
if (i <= m)
{
for (k = i; k<=m; k++)
{
scale1 = scale1 + fabs(a[(k-1)*n+i]);
}
if (scale1 != 0.0)
{
for (k = i; k<=m; k++)
{
a[(k-1)*n+i] = a[(k-1)*n+i] / scale1;
s = s + a[(k-1)*n+i] * a[(k-1)*n+i];
}
f = a[(i-1)*n+i];
g = -sqrt(s) * sgn(f);
h = f * g - s;
a[(i-1)*n+i] = f - g;
if (i != n)
{
for (j = l; j<=n; j++)
{
s = 0.0;
for (k = i; k<=m; k++)
{
s = s + a[(k-1)*n+i] * a[(k-1)*n+j];
}
f = s / h;
for (k = i; k<=m; k++)
{
a[(k-1)*n+j] = a[(k-1)*n+j] + f * a[(k-1)*n+i];
}
}
}
for (k = i; k<=m; k++)
{
a[(k-1)*n+i] = scale1 * a[(k-1)*n+i];
}
}
}
w[i] = scale1 * g;
g = 0.0;
s = 0.0;
scale1 = 0.0;
if ((i <= m) && (i != n))
{
for (k = l; k<=n; k++)
{
scale1 = scale1 + fabs(a[(i-1)*n+k]);
}
if (scale1 != 0.0)
{
for (k = l; k<=n; k++)
{
a[(i-1)*n+k] = a[(i-1)*n+k] / scale1;
s = s + a[(i-1)*n+k] * a[(i-1)*n+k];
}
f = a[(i-1)*n+l];
g = -sqrt(s) * sgn(f);
h = f * g - s;
a[(i-1)*n+l] = f - g;
for (k = l; k<=n; k++)
{
rv1[k] = a[(i-1)*n+k] / h;
}
if (i != m)
{
for (j = l; j<=m; j++)
{
s = 0.0;
for (k = l; k<=n; k++)
{
s = s + a[(j-1)*n+k] * a[(i-1)*n+k];
}
for (k = l; k<=n; k++)
{
a[(j-1)*n+k] = a[(j-1)*n+k] + s * rv1[k];
}
}
}
for (k = l; k<=n; k++)
{
a[(i-1)*n+k] = scale1 * a[(i-1)*n+k];
}
}
}
if (anorm > (fabs(w[i]) + fabs(rv1[i])))
{
anorm = anorm;
}
else
{
anorm = fabs(w[i]) + fabs(rv1[i]);
}
}
for (i=n; i>=1; i--)
{
if (i < n)
{
if (g != 0.0)
{
for (j = l; j<=n; j++)
{
v[(j-1)*n+i] = (a[(i-1)*n+j] / a[(i-1)*n+l]) / g;
}
for (j = l; j<=n; j++)
{
s = 0.0;
for (k = l; k<=n; k++)
{
s = s + a[(i-1)*n+k] * v[(k-1)*n+j];
}
for (k = l; k<=n; k++)
{
v[(k-1)*n+j] = v[(k-1)*n+j] + s * v[(k-1)*n+i];
}
}
}
for (j = l; j<=n; j++)
{
v[(i-1)*n+j] = 0.0;
v[(j-1)*n+i] = 0.0;
}
}
v[(i-1)*n+i] = 1.0;
g = rv1[i];
l = i;
}
for (i = n; i>=1; i--)
{
l = i + 1;
g = w[i];
if (i < n)
{
for (j = l; j<=n; j++)
{
a[(i-1)*n+j] = 0.0;
}
}
if (g != 0.0)
{
g = 1.0 / g;
if (i != n)
{
for (j = l; j<=n; j++)
{
s = 0.0;
for (k = l; k<=m; k++)
{
s = s + a[(k-1)*n+i] * a[(k-1)*n+j];
}
f = (s / a[(i-1)*n+i]) * g;
for (k = i; k<=m; k++)
{
a[(k-1)*n+j] = a[(k-1)*n+j] + f * a[(k-1)*n+i];
}
}
}
for (j= i; j<=m; j++)
{
a[(j-1)*n+i] = a[(j-1)*n+i] * g;
}
}
else
{
for (j= i; j<=m; j++)
{
a[(j-1)*n+i] = 0.0;
}
}
a[(i-1)*n+i] = a[(i-1)*n+i] + 1.0;
}
for (k = n; k>=1; k--)
{
for (its = 1; its<=30; its++)
{
for (l = k; l>=1; l--)
{
nm = l - 1;
if (fabs(rv1[l]) + anorm == anorm)
{
goto r2;
}
if (fabs(w[nm]) + anorm == anorm)
{
goto r1;
}
}
r1: c = 0.0;
s = 1.0;
for (i = l; i<=k; i++)
{
f = s * rv1[i];
if (fabs(f) + anorm != anorm)
{
g = w[i];
h = sqrt(f * f + g * g);
w[i] = h;
h = 1.0 / h;
c = (g * h);
s = -(f * h);
for (j = 1; j<=m; j++)
{
y = a[(j-1)*n+nm];
z = a[(j-1)*n+i];
a[(j-1)*n+nm] = (y * c) + (z * s);
a[(j-1)*n+i] = -(y * s) + (z * c);
}
}
}
r2: z = w[k];
if (l == k)
{
if (z < 0.0)
{
w[k] = -z;
for (j = 1; j<=n; j++)
{
v[(j-1)*n+k] = -v[(j-1)*n+k];
}
}
goto r3;
}
if (its == 30)
{
cout<< "no convergence in 30 iterations"<<endl;
}
x = w[l];
nm = k - 1;
y = w[nm];
g = rv1[nm];
h = rv1[k];
f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
g = sqrt(f * f + 1.0);
f = ((x - z) * (x + z) + h * ((y / (f + fabs(g) * sgn(f))) - h)) / x;
c = 1.0;
s = 1.0;
for (j = l; j<=nm; j++)
{
i = j + 1;
g = rv1[i];
y = w[i];
h = s * g;
g = g * c;
z = sqrt(f * f + h * h);
rv1[j] = z;
c = f / z;
s = h / z;
f = (x * c) + (g * s);
g = -(x * s) + (g * c);
h = y * s;
y = y * c;
for (pp = 1; pp<=n; pp++)
{
x = v[(pp-1)*n+j];
z = v[(pp-1)*n+i];
v[(pp-1)*n+j] = (x * c) + (z * s);
v[(pp-1)*n+i] = -(x * s) + (z * c);
}
z = sqrt(f * f + h * h);
w[j] = z;
if (z != 0.0)
{
z = 1.0 / z;
c = f * z;
s = h * z;
}
f = (c * g) + (s * y);
x = -(s * g) + (c * y);
for (pp = 1; pp<=m; pp++)
{
y = a[(pp-1)*n+j];
z = a[(pp-1)*n+i];
a[(pp-1)*n+j] = (y * c) + (z * s);
a[(pp-1)*n+i] = -(y * s) + (z * c);
}
}
rv1[l] = 0.0;
rv1[k] = f;
w[k] = x;
}
r3: aaaaa = 1;
}
}
void fleg(double x, double pl[], int nl)
{
pl[1] = 1.0;
pl[2] = x;
if (nl > 2)
{
double twox = 2.0 * x;
double f2 = x;
double d = 1.0;
for (int j = 3; j<=nl; j++)
{
double f1 = d;
f2 = f2 + twox;
d = d + 1.0;
pl[j] = (f2 * pl[j - 1] - f1 * pl[j - 2]) / d;
}
}
}
void fpoly(double x, double p[], int np)
{
p[1] = 1.0;
for (int j = 2; j<=np; j++)
{
p[j] = p[j - 1] * x;
}
}
void svdfit(double x[], double y[], double sig[], int ndata, double a[],
int ma, double u[], double v[], double w[], int mp, int np,
double& chisq, char* funcs)
{
int i, j;
double thresh,tmp,tol = 0.00001;
double b[1001], afunc[51];
for (i = 1; i<=ndata; i++)
{
if (!strcmp(funcs,"fpoly"))
{
fpoly(x[i], afunc, ma);
}
if (!strcmp(funcs,"fleg"))
{
fleg(x[i], afunc, ma);
}
tmp = 1.0 / sig[i];
for (j = 1; j<=ma; j++)
{
u[(i-1)*ma+j] = afunc[j] * tmp;
}
b[i] = y[i] * tmp;
}
svdcmp(u, ndata, ma, w, v);
double wmax = 0.0;
for (j = 1; j<=ma; j++)
{
if (w[j] > wmax)
{
wmax = w[j];
}
}
thresh = tol * wmax;
for (j = 1; j<=ma; j++)
{
if (w[j] < thresh)
{
w[j] = 0.0;
}
}
svbksb(u, w, v, ndata, ma, b, a);
chisq = 0.0;
for (i = 1; i<=ndata; i++)
{
if (!strcmp(funcs,"fpoly"))
{
fpoly(x[i], afunc, ma);
}
if (!strcmp(funcs,"fleg"))
{
fleg(x[i], afunc, ma);
}
double sum1 = 0.0;
for (j = 1; j<=ma; j++)
{
sum1 = sum1 + a[j] * afunc[j];
}
chisq = chisq + ((y[i] - sum1) / sig[i]) * ((y[i] - sum1) / sig[i]);
}
}
void svdvar(double v[], int ma, int np, double w[], double cvm[], int ncvm)
{
double wti[21];
int i,j;
for (i = 1; i<=ma; i++)
{
wti[i] = 0.0;
if (w[i] != 0.0)
{
wti[i] = 1.0 / (w[i] * w[i]);
}
}
for (i = 1; i<=ma; i++)
{
for (j = 1; j<=i; j++)
{
double sum1 = 0.0;
for (int k = 1; k<=ma; k++)
{
sum1 = sum1 + v[(i-1)*np+k] * v[(j-1)*np+k] * wti[k];
}
cvm[(i-1)*ma+j] = sum1;
cvm[(j-1)*ma+i] = sum1;
}
}
}
double ran1(long& idum)
{
int j,iff=-1;
static long ix1,ix2,ix3;
static double r[98];
long m1 = 259200; long m2 = 134456; long m3 = 243000;
long ia1 = 7141; long ia2 = 8121; long ia3 = 4561;
long ic1 = 54773; long ic2 = 28411; long ic3 = 51349;
double rm1 = 0.0000038580247; double rm2 = 0.0000074373773;
if (idum < 0 || iff == 0)
{
iff = 1;
ix1 = (ic1 - idum) % m1;
ix1 = (ia1 * ix1 + ic1) % m1;
ix2 = ix1 % m2;
ix1 = (ia1 * ix1 + ic1) % m1;
ix3 = ix1 % m3;
for (j = 1; j<=97; j++)
{
ix1 = (ia1 * ix1 + ic1) % m1;
ix2 = (ia2 * ix2 + ic2) % m2;
r[j] = (double(ix1) + double(ix2) * rm2) * rm1;
}
idum = 1;
}
ix1 = (ia1 * ix1 + ic1) % m1;
ix2 = (ia2 * ix2 + ic2) % m2;
ix3 = (ia3 * ix3 + ic3) % m3;
j = 1 + int((97 * ix3) / m3);
if (j > 97 || j < 1)
{
cout<<"abnormal exit in ran1"<<endl;
exit(1);
}
double temp=r[j];
r[j] = (double(ix1) + double(ix2) * rm2) * rm1;
return temp;
}
double gasdev(long& idum)
{
static int iset;
static double gset;
double v1,v2,r,fac;
if (iset == 0)
{
do
{
v1 = 2.0 * ran1(idum) - 1.0;
v2 = 2.0 * ran1(idum) - 1.0;
r = v1 * v1 + v2 * v2;
}while (r >= 1.0 || r == 0);
fac = sqrt(-2.0 * log(r) / r);
gset = v1 * fac;
iset = 1;
return v2 * fac;
}
else
{
iset = 0;
return gset;
}
}
void main()
{
//program d9r4
//driver for routine svdfit
int i,mp,np,npt = 100;
double chisq;
double spread = 0.02;
int npol = 5;
double x[101], y[101], sig[101], a[6], cvm[26];
double u[501], v[26], w[6];
//polynomial fit
long idum = -911;
mp = npt;
np = npol;
for (i = 1; i<=npt; i++)
{
x[i] = 0.02 * i;
y[i] = 1.0 + x[i] * (2.0 + x[i] * (3.0 + x[i] * (4.0 + x[i] * 5.0)));
y[i] = y[i] + spread * gasdev(idum);
sig[i] = spread;
}
svdfit(x, y, sig, npt, a, npol, u, v, w, mp, np, chisq, "fpoly");
svdvar(v, npol, np, w, cvm, npol);
cout<<endl;
cout<<"Polynomial fit"<<endl;
cout<<endl;
cout.setf(ios::fixed|ios::right);
cout.precision(6);
for (i = 1; i<=npol; i++)
{
cout.width(12);
cout<<a[i]<<" +-";
cout.width(12);
cout<<sqrt(cvm[(i-1)*npol+i])<<endl;
}
cout<<endl;
cout<<"Chi-squared ";
cout<<chisq<<endl;
svdfit(x, y, sig, npt, a, npol, u, v, w, mp, np, chisq, "fleg");
svdvar(v, npol, np, w, cvm, npol);
cout<<endl;
cout<<"Legendre polynomial fit"<<endl;
for (i = 1; i<=npol; i++)
{
cout.width(12);
cout<<a[i]<<" +-";
cout.width(12);
cout<<sqrt(cvm[(i-1)*npol+i])<<endl;
}
cout<<endl;
cout<<"Chi-squared "<<chisq<<endl;
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -