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

//   nodes - узлы, между которыми должны группироваться точки

//   npoints - кол-во элементов в nodes

//   alpha - коэффициент сглаживания

//  Выходные

//   sy - сглаженный вариант кривой

//

// Коды возврата

//  0 - не удалось построить сплайн, система вырождена

//  1 - все в порядке

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

double *sy, char* IntMem );

#endif

#include <math.h>

#include <mem.h>

#include "misc.h"

const double TinyPower = 1e-50;

int LUSolve( int n, double *A, double *F, double *X )

{

int i,k,r;

double sum;

double tiny = fabs(A[0]);

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

tiny = tiny*TinyPower;

if( fabs(A[0]) < tiny )return 0;

for( i = 1; i < n; i++ )A[i] /= A[0];

for( r = 1; r < n; r++ )

{

sum = 0;

for( k = 0; k < r; k++ )sum += A[r*n+k]*A[k*n+r];

A[r*n+r] -= sum;

if( fabs(A[r*n+r]) < tiny )

return 0;

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

{

sum = 0;

for( k = 0; k < r; k++ )sum += A[i*n+k]*A[k*n+r];

A[i*n+r] -= sum;

sum = 0;

for( k = 0; k < r; k++ )sum += A[r*n+k]*A[k*n+i];

A[r*n+i] = (A[r*n+i] - sum)/A[r*n+r];

}

}

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

{

sum = 0;

for( i = 0; i < r; i++ )sum += A[r*n+i]*X[i];

X[r] = (F[r]-sum)/A[r*n+r];

}

for( r = n - 1; r >= 0; r-- )

{

sum = 0;

for( i = r + 1; i < n; i++ )sum += A[r*n+i]*X[i];

X[r] = X[r]-sum;

}

return 1;

}

int LLSolve( int n, double *A, double *F, double *X )

{

int i,j,k;

double root,sum;

// Теперь массив начинается с индекса 1

A--; F--; X--;

// -----------------------------------double tiny=fabs(A[1]);

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

tiny = tiny*TinyPower;

// Count L (store in A ) A=L*L'

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

{

sum = 0;

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

sum += A[i*(i-1)/2+k]*A[i*(i-1)/2+k];

root = A[i*(i+1)/2] - sum;

if( root < tiny )

return 0;

A[i*(i+1)/2] = sqrt(root);

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

{

sum = 0;

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

sum += A[i*(i-1)/2+k]*A[j*(j-1)/2+k];

A[j*(j-1)/2+i] = (A[j*(j-1)/2+i] - sum)/A[i*(i+1)/2];

}

}

// Count  Y (store in X)

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

{

sum = 0;

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

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

X[i]=(F[i]-sum)/A[i*(i+1)/2];

}

// Count X

for( i = n; i >= 1; i-- )

{

sum = 0;

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

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

X[i]=(X[i]-sum)/A[i*(i+1)/2];

}

return 1;

}

int CholeskySR( int n, double *A, double *betta, int nx, double *L )

{

int i,j,k;

double root,sum;

if( L == NULL )L = A;

double tiny=fabs(A[0]);

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

tiny = tiny*TinyPower;

// Count L A=L*L'

// First column

root = A[0] + betta[0];

if( root < tiny )return 0;

L[0] = sqrt(root);

L[1] = (A[1] - betta[0])/L[0];

for( j = 2; j < n; j++ )

L[j*(j+1)/2] = A[j*(j+1)/2]/L[0];

// Others

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

{

sum = 0;

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

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

root = A[i*(i+3)/2] + (betta[i-1] + betta[i]) - sum;

if( root < tiny )

return 0;

L[i*(i+3)/2] = sqrt(root);

if( i == n-1 )break;

sum = 0;

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

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

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

for( j = i+2; j < n; j++ )

{

sum = 0;

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

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

L[j*(j+1)/2+i] = (A[j*(j+1)/2+i] - sum)/L[i*(i+3)/2];

}

}

return 1;

}

int CholeSolveS( int n, double *L, double *F, int nx, double *X )

{

int i,k;

double sum;

// Count  Y (store in X)

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

{

sum = 0;

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

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

X[i]=(F[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;

}

int CholeSolveSR( int n, double *L, double *betta, double *F, double *rf, int nx, double *X )