Разложение импульса.Разложение ЭДС.Использование регуляризации, страница 20

{

int i,k;

double sum;

// Count  Y (store in X)

X[0]=(F[0]-betta[0]*(rf[1]-rf[0]))/L[0];

for( i = 1; i < nx; i++ )

{

sum = 0;

for( k = 0; k < i; k++ )

sum += L[i*(i+1)/2+k]*X[k];

X[i]=(F[i]+(betta[i-1]*(rf[i]-rf[i-1])-betta[i]*(rf[i+1]-rf[i]))-sum)/L[i*(i+3)/2];

}

// Count X

for( i = nx-1; i >= 0; i-- )

{

sum = 0;

for( k = i+1; k < n; k++ )

sum += L[k*(k+1)/2+i]*X[k];

X[i]=(X[i]-sum)/L[i*(i+3)/2];

}

return 1;

}

// smooth cubic spline

int LLSolveTape( int n, int m, double *A, double *F, double *X )

{

int i,j,k;

double sum,root;

int r, de = m-1;

double tiny=fabs(A[0]);

for( i = 1; i < n*m; i++ )tiny = tiny < fabs(A[i]) ? fabs(A[i]) : tiny;

tiny = tiny*TinyPower;

for( i = 0; i < n; i++ )

{

sum = 0;

for( k = 0; k < de; k++ )sum += A[i*m+k]*A[i*m+k];

root = A[i*m+de] - sum;

if( root < tiny )return 0;

A[i*m+de] = sqrt(root);

r = i+1;

for( j = de-1; j >= 0 && r < n; j-- )

{

sum = 0;

for( k = 0; k < j; k++ )

sum += A[r*m+k]*A[i*(m-1)+r+k];

A[r*m+j] = (A[r*m+j] - sum)/A[i*m+de];

r++;

}

}

// ----------------------------------for( i = 0; i < m; i++ )

{

sum = 0;

for( k = de-i; k < de; k++ )

sum += A[i*m+k]*X[i-de+k];

X[i] = (F[i] - sum)/A[i*m+de];

}

for( i = m; i < n; i++ )

{

sum = 0;

for( k = 0; k < de; k++ )

sum += A[i*m+k]*X[i-de+k];

X[i] = (F[i] - sum)/A[i*m+de];

}

// ---------------------------------X[n-1] /= A[n*m-1];

for( i = n-2; i > n-m; i-- )

{

sum = 0;

r = i+1;

for( k = de-1; k >= m+i-n; k-- )

{

sum += A[r*m+k]*X[r];

r++;

}

X[i] = (X[i] - sum)/A[i*m+de];

}

for( i = n-m; i >= 0; i-- )

{

sum = 0;

r = i+1;

for( k = de-1; k >= 0; k-- )

{

sum += A[r*m+k]*X[r];

r++;

}

X[i] = (X[i] - sum)/A[i*m+de];

}

return 1;

}

double fi1(double ksi) { double ksi2=ksi*ksi; return 1.0 - 3.0*ksi2 + 2.0*ksi*ksi2; }

double fi2(double ksi) { double ksi2=ksi*ksi; return ksi - 2.0*ksi2 + ksi*ksi2; }

double fi3(double ksi) { double ksi2=ksi*ksi; return 3.0*ksi2 - 2.0*ksi*ksi2; }

double fi4(double ksi) { double ksi2=ksi*ksi; return -ksi2 + ksi*ksi2; }

int SmoothSpline( double *x, double *y, int n, double *nodes, int npoints, double alfa,

double *sy, char* IntMem )

{

// make matrix

double cy,cx,ksi,d1,d2,d3,d4,hk;

int j,k;

int nslae = npoints*2;

int SelfAlloc = 0;

double *A, *F;

if( IntMem == NULL )

{

IntMem = new char[nslae*5*sizeof(double)];

SelfAlloc = 1;

}

A = (double *)&IntMem[0];

F = (double *)&IntMem[nslae*4*sizeof(double)];

memset( A, 0, nslae*4*sizeof(double) );

memset( F, 0, nslae*sizeof(double) );

for( k=0; k < n; k++)

{

cx = x[k];

if(cx<nodes[0] || cx>nodes[npoints-1]) continue;

for(j=1;j<npoints;j++) if(cx<=nodes[j]) break;

j--;

cy = y[k];

hk=nodes[j+1]-nodes[j];

ksi=(cx-nodes[j])/hk;

d1=fi1(ksi);

d2=fi2(ksi)*hk;

d3=fi3(ksi);

d4=fi4(ksi)*hk;

j*=2;

A[j*4+3] += d1*d1;

A[(j+1)*4+2] += d2*d1;

A[(j+1)*4+3] += d2*d2;

A[(j+2)*4+1] += d3*d1;

A[(j+2)*4+2] += d3*d2;

A[(j+2)*4+3] += d3*d3;

A[(j+3)*4] += d4*d1;

A[(j+3)*4+1] += d4*d2;

A[(j+3)*4+2] += d4*d3;

A[(j+3)*4+3] += d4*d4;

F[j] += d1*cy;

F[j+1] += d2*cy;

F[j+2] += d3*cy;

F[j+3] += d4*cy;

}

for( k=0; k<npoints-1; k++)

{

hk = nodes[k+1]-nodes[k];

double hk3=hk*hk*hk;

j=k*2;

A[j*4+3] += alfa*12.0/hk3;

A[(j+1)*4+2] += alfa*6.0/hk;

A[(j+1)*4+3] += alfa*4.0*hk;

A[(j+2)*4+1] += -alfa*12.0/hk3;

A[(j+2)*4+2] += -alfa*6.0/hk;

A[(j+2)*4+3] += alfa*12.0/hk3;

A[(j+3)*4] += alfa*6.0/hk;

A[(j+3)*4+1] += alfa*2.0*hk;

A[(j+3)*4+2] += -alfa*6.0/hk;

A[(j+3)*4+3] += alfa*4.0*hk;

}

int res = LLSolveTape( nslae, 4, A, F, F );

if( res )

{

for( k=0; k<n; k++ )

{

cx = x[k];

if( cx > nodes[npoints-1] || cx < nodes[0])

{

sy[k] = 0.0;

continue;

}

for(j=1;j<npoints;j++) if(cx<=nodes[j]) break;

j--;

hk = nodes[j+1] - nodes[j];

ksi = (cx - nodes[j])/hk;

j *= 2;

sy[k] = F[j]*fi1(ksi) +

F[j+1]*fi2(ksi)*hk +

F[j+2]*fi3(ksi) +

F[j+3]*fi4(ksi)*hk;

}

}

if( SelfAlloc )delete[] IntMem;

return res;

}