Си, метод минимальной невязки для решения СЛАУ

[ajaxtpm]

Форумчане и математики, прошу вашей помощи :) Почти что вторую неделю я пытаюсь разобрать, что в моей программе не работает. В начале по заданию не использовал массивы вообще и считал всё по сложным штукам, но потом забил и упростил, стало полегче, но он не химичит всё равно ><
Есть уравнения вида :

$X_1 = X_2$
$A(i)(i-1)X_{i-1} + A(i)(i)X_{i} + A(i)(i+1)X_{i+1} = Bi , i = 2..N-1$
$Xn = 10$

Где матрица А - трёхдиагональная (ну почти) и коэфициенты её я организовал в отдельной функции, а B - интегралы, зависящие от i. Всё круто расписано, организованно, но никак не желает меня радовать, помогите кто может, а ? :)
В качестве проверки по заданию положенно использовать простенький метод обратной матрицы, делал через матлаб, ответы матлаба и моей проги расходятся.

Вот ссылка на метод :)
http://www.intuit.ru/department/calculate/calcmathbase/2/9.html

Integr.m

function I = Integr( J , T )
I = -( log( 10 + J + cos( 1 + exp( T ) ) ) ) ; 
end

main.m

N = 5 ;
A = zeros( N , N ) ;
A(1,1)=1;
A(1,2)=-1;
for i = 2 : N-1
    A(i,i-1)=1+i/2
    A(i,i)=-(i*i+1)
    A(i,i+1)=2*i-1
    B(i,1)=quad(@(T)Integr(i,T),1,3)
end
A(N,N)=1
B(N,1)=10
X=A\B
%plot(X,'.-')
plot(X)

Ну и сама моя прога :

Используется синтаксис C

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

#define N 5
#define XLast 10. 
#define Accuracy 0.01
#define IntegralStart 1.
#define IntegralEnd 3. 

//double Norm = ( double ) ( N * N ) ;
double Norm = 1. ;
double Tau ;
double X[ N ] ;
double R[ N ] ;
double B[ N ] ;
double AR[ N ] ;

/* Calculate coefficient A0 , A1 , A2 */
/* ================================== */
double GetAElement( int I , int J ) {
       if ( I == 1 && J == 1 ) return 1. / Norm ;
       if ( I == 1 && J == 2 ) return -1. / Norm ;
       if ( I == N && J == N ) return 1. / Norm ;
       if ( I == N && J == N - 1 ) return 0. ;
       if ( J == I - 1 ) return ( ( double ) I / 2. + 1. ) / Norm ; 
       if ( J == I ) return -1. * ( ( double ) ( I * I ) + 1. ) / Norm ; 
       if ( J == I + 1 ) return ( 2. * ( double ) I - 1. ) / Norm ; 
       return 0. ;
} ;

void MakeNorm() {
     int I , J ;
     double Res = 0. ;
     for ( I = 0 ; I < N ; I++ ) 
         for ( J = 0 ; J < N ; J++ ) {
             Res += GetAElement( I + 1 , J + 1 ) * GetAElement( I + 1 , J + 1 ) ;
//             printf( " Element ( %i ; %i ) = %g      Res = %g\n" , I + 1 , J + 1 , GetAElement( I + 1 , J + 1 ) , Res ) ;
         } ;
     Norm = sqrt( Res ) ;
     printf( " Norm = %g\n" , Norm ) ; scanf( "%*c" ) ;
} ;
     
/* Calculate integral function at X point and I parametr */
/* ===================================================== */
double IFunc( double X, int I ) {
       return -1. * log( 10. + ( double ) I + cos( 1. + exp( X ) ) ) ;
} ;

/* Calculate Integral at I parametr and save it into B array */
/* ========================================================= */
void IntegralCalc() {
     int I ;
     double T = IntegralStart ;
     double Res = 0. ;
     for ( I = 1 ; I < N - 1 ; I++ ) {
         while ( T < IntegralEnd ) {
               Res += ( IFunc( T , I ) + IFunc( T + Accuracy , I ) ) * Accuracy / 2. ;             
               T += Accuracy ;
         } ;
         B[ I ] = Res / Norm ;
     } ;
     B[ 0 ] = 0. ; B[ N - 1 ] = XLast / Norm ;
} ;

void Initialize() { int I ; for ( I = 0 ; I < N ; I++ ) { X[ I ] = 1. ; R[ I ] = 0. ; } ; } ;

void CalcRArray() {
     int I ;
     
     R[ 0 ] = GetAElement( 1 , 1 ) * X[ 0 ] + GetAElement( 1 , 2 ) * X[ 1 ] ;
     for ( I = 1 ; I < N - 1 ; I++ ) R[ I ] = GetAElement( I + 1 , I ) * X[ I - 1 ] + GetAElement( I + 1 , I + 1 ) * X[ I ] + GetAElement( I + 1 , I + 2 ) * X[ I + 1 ] ;
     R[ N - 1 ] = GetAElement( N , N ) * X[ N - 1 ] ;
     
     for ( I = 0 ; I < N ; I++ ) R[ I ] -= B[ I ] ;
} ;

void CalcARArray() {
     int I ;
     
     AR[ 0 ] = GetAElement( 1 , 1 ) * R[ 0 ] + GetAElement( 1 , 2 ) * R[ 1 ] ;
     for ( I = 1 ; I < N - 1 ; I++ ) AR[ I ] = GetAElement( I + 1 , I ) * R[ I - 1 ] + GetAElement( I + 1 , I + 1 ) * R[ I ] + GetAElement( I + 1 , I + 2 ) * R[ I + 1 ] ;
     AR[ N - 1 ] = GetAElement( N , N ) * R[ N - 1 ] ;
     
} ;     

void CalcTau() {
     int I ;
     double Ch = 0. ;
     double Zn = 0. ;
     
     for ( I = 0 ; I < N ; I++ ) {
         Ch += AR[ I ] * R[ I ] ;
         Zn += AR[ I ] * AR[ I ] ;
     } ;

     Tau = Ch / Zn ;
     printf( " ## Tau = %g\n" , Tau ) ;
} ;

void CalcNewX() { int I ; for ( I = 0 ; I < N ; I++ ) X[ I ] -= Tau * R[ I ] ; } ;

void PrintX() { int I ; for ( I = 0 ; I < N ; I++ ) printf( "   %g\n" , X[ I ] ) ; printf( "\n\n" ) ; } ;


int main() {
    int I = 0 , J ;
    
    IntegralCalc() ;
    Initialize() ;
//    MakeNorm() ;

    while ( 1 == 1 ) {
        I++ ;
        CalcRArray() ;
        CalcARArray() ;
        CalcTau() ;
        CalcNewX() ;
        if ( I % 1000 == 0 ) { 
             PrintX() ; scanf( "%*c" ) ; 
        } ;
    } ;
/**/    
    
    scanf( "%*c" ) ;            
    return 0 ;
} ;
 

[Toucan]

ajaxtpm в сообщении #574383 писал(а):

Чуваки, это катастрофа

Не совсем правильно начинать общение на форуме с подобного обращения. Придумайте, пожалуйста, что-нибудь другое, а тема пока полежит в Карантине.

Для редактирования темы пользуйтесь кнопкой "Правка".

После того как исправите сообщение, сообщите об этом в теме Сообщение в карантине исправлено.

[Toucan]

Так лучше. Вернул.

[hurtsy]

ajaxtpm в сообщении #574383 писал(а):

Всё круто расписано, организованно, но никак не желает меня радовать, помогите кто может, а ? :)
В качестве проверки по заданию положенно использовать простенький метод обратной матрицы, делал через матлаб, ответы матлаба и моей проги расходятся.

Цитата:

Не надо печалиться ...

Ваш итерационный процесс стабилизировался, где то после десятой итерации. После двадцатой итерации машина трудится вхолостую. А как вы обращаете трех(ну почти)диагональную матрицу, интересно однако. С уважением,

[ajaxtpm]

Задание решено, пускай я буду тем, кто даст код на халяву
Кстати говоря, система из миллиона уравнений прекрасно решится минут за .. 5 )

Используется синтаксис C

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

#define N 10000
#define Accuracy 0.01


double Delta ;
double Tau = 1. ; // Tau parametr of system
double X[ N ] ; // Solution of system
double R[ N ] ; // R array
double AR[ N ] ;// A*R array
double B[ N ] ;  // Integral array


/* Function to Get needed A array value */
/* ==================================== */
double GetA ( int I , int J ) {
       double Res ;
       double Aver = -( ( double ) ( I * I ) + 1. ) ;
       
       if ( I == 1 && J == 1 ) Res = 0.5 ; // First A row
       else if ( I == 1 && J == 2 ) Res = -0.5 ;// First A row
       else if ( I == N && J == N ) Res = 1. ; // Last A row
       
       else if ( I < N && J == I - 1 ) Res = ( 1. + ( double ) I / 2. ) / ( -( ( double ) ( I * I ) + 1. ) ) ;
       else if ( I < N && J == I ) Res = 1. ;
       else if ( I < N && J == I + 1 ) Res = ( ( double ) I * 2. - 1. ) / ( -( ( double ) ( I * I ) + 1. ) ) ;
              
       else Res = 0. ;
       return Res ;
} ;


/* Calculate Integral Of Partition ( IOP ) size */
/* ==================================== */
double IOP( double P , int I ) {
       /* Calculate UnderIntegral Function ( UIF ) by I and T */
       /* =================================================== */
       double UIF( int I , double T ) { return log( 10. + ( double ) I + cos( 1. + exp( T ) ) ) ; } ;

       double Delta = 2. / P ;
       double Sum = 0. ;
       double T = 1. ;
       
       while ( T < 3. ) {
             Sum += Delta * ( UIF( I , T ) + UIF( I , T + Delta ) ) / 2. ;
             T += Delta ;
       } ;
       return Sum ;
} ;

/* Fill B array by values and integrals */
/* ==================================== */
void CalcB() {
     double Partition = 101. ;
     double Delta ;
     double Sum1 , Sum2 ;
     int I ; 
     
     B[ 0 ] = 0. ;
     B[ N - 1 ] = 10. ;
     for ( I = 1 ; I < N - 1 ; I++ ) {
         Sum1 = IOP( Partition - 1. , I ) ;
         Sum2 = IOP( Partition , I ) ;
         while ( fabs( Sum2 - Sum1 ) > Accuracy ) {
               Sum1 = Sum2 ;
               Partition += 1. ;
               Sum2 = IOP( Partition , I ) ;
         } ; 
         B[ I ] = Sum2 / ( -( ( double ) ( ( I + 1 ) * ( I + 1 ) ) + 1. ) ) ;
     } ;
} ;

void CalcR() {
     int I , J ;
     
     R[ 0 ] = GetA( 1 , 1 ) * X[ 0 ] + GetA( 1 , 2 ) * X[ 1 ] - B[ 0 ] ;
     R[ N - 1 ] = GetA( N , N ) * X[ N - 1 ] - B[ N - 1 ] ;
     for ( I = 1 ; I < N - 1 ; I++ ) {
         R[ I ] = 0. ;
         for ( J = I - 1 ; J < I + 2 ; J++ ) R[ I ] += GetA( I + 1 , J + 1 ) * X[ J ] ;
         R[ I ] -= B[ I ] ;
     } ;
} ;

void CalcAR() {
     int I , J ;

     AR[ 0 ] = GetA( 1 , 1 ) * R[ 0 ] + GetA( 1 , 2 ) * R[ 1 ] ;
     AR[ N - 1 ] = GetA( N , N ) * R[ N - 1 ] ;
     for ( I = 1 ; I < N - 1 ; I++ ) {
         AR[ I ] = 0. ;
         for ( J = I - 1 ; J < I + 2 ; J++ ) AR[ I ] += GetA( I + 1 , J + 1 ) * R[ J ] ;
     } ;

} ;

void CalcTau() {
     int I ;
     double Ch = 0. , Zn = 0. ;
     
     for ( I = 0 ; I < N ; I++ ) {
         Ch += AR[ I ] * R[ I ] ;
         Zn += AR[ I ] * AR[ I ] ;
     } ;
     
     Tau = Ch / Zn ;
} ;

void CalcX() {
     double NewX[ N ] ; 
     int I , J ; 
     NewX[ 0 ] = ( 1. - Tau * GetA( 1 , 1 ) ) * X[ 0 ] - Tau * GetA( 1 , 2 ) * X[ 1 ] + Tau * B[ 0 ] ;
     NewX[ N - 1 ] = ( 1. - Tau * GetA( N , N ) ) * X[ N - 1 ] + Tau * B[ N - 1 ] ;
     for ( I = 1 ; I < N - 1 ; I++ ) {
         NewX[ I ] = -Tau * GetA( I + 1 , I ) * X[ I - 1 ] + ( 1. - Tau * GetA( I + 1 , I + 1 ) ) * X[ I ] - Tau * GetA( I + 1 , I + 2 ) * X[ I + 1 ]  + Tau * B[ I ] ;  
     } ;
     Delta = 0. ;
     for ( I = 0 ; I < N ; I++ ) { if( fabs( NewX[ I ] - X[ I ] ) > Delta ) Delta = fabs( NewX[ I ] - X[ I ] ) ; X[ I ] = NewX[ I ] ; } ;
} ;

/* Functions to out Solution */
/* ======================= */
void PrintX() { int I ; for ( I = 0 ; I < N ; I++ ) printf( "  %16.14lf\n" , X[ I ] ) ; printf( "\n" ) ; } ;


int main() {
    int I = 0 ;

    CalcB() ;
    for ( I = 0 ; I < N ; I++ ) X[ I ] = 1. ;

    do {
        CalcR() ; CalcAR() ; CalcTau() ; CalcX() ;        
    } while ( Delta > 0.001 ) ;
    PrintX() ;
    
    scanf( "%*c" ) ;
    return 0 ;
} ;
 

[hurtsy]

ajaxtpm в сообщении #575446 писал(а):

Задание решено, пускай я буду тем, кто даст код на халяву
Кстати говоря, система из миллиона уравнений прекрасно решится минут за .. 5 )

Судя по косоватости и кривоватости ваших принтов(выводов на консоль), прогноз не верен. Надо поработать а потом одаряйте форум и человечество. С уважением,

[ajaxtpm]

Знаете ли, кем бы вы ни были, такое изречение от вас кажется мне несерьёзным.
Напиши поконкретней, т.к. программа выполняет все требуемые действия, либо не пишите о кривоватости вообще :)

[hurtsy]

ajaxtpm в сообщении #576262 писал(а):

Напиши поконкретней, т.к. программа выполняет все требуемые действия, либо не пишите о кривоватости вообще :)

Ну, если серьезно и только для вас. Выбрать вывод в один столбик там где строчек одновременно нельзя рассмотреть больше 300, и рассуждать о миллионах элементов. Задача просмотра даже тысячи значений требует какого-то программного решения. Такие подарки и даром не нужны. С уважением,

[Maslov]

Решение системы из тысяч уравнений само по себе редко когда интересно: обычно это часть какой-нибудь другой задачи. Вывод результатов в данном случае носит чисто вспомогательный характер и может легко быть переделан под конкретные нужды конкретной "объемлющей" программы.

(Оффтоп)

hurtsy в сообщении #576600 писал(а):

Такие подарки и даром не нужны.

Есть предложение делать менее категоричные заявления. Ну, например, "мне такие подарки и даром не нужны". А то от лица человечества -- как-то нескромно.

[ajaxtpm]

hurtsy в сообщении #576600 писал(а):

ajaxtpm в сообщении #576262 писал(а):

Напиши поконкретней, т.к. программа выполняет все требуемые действия, либо не пишите о кривоватости вообще :)

Ну, если серьезно и только для вас. Выбрать вывод в один столбик там где строчек одновременно нельзя рассмотреть больше 300, и рассуждать о миллионах элементов. Задача просмотра даже тысячи значений требует какого-то программного решения. Такие подарки и даром не нужны. С уважением,

С уважением, глупости говорите)
Ясно же, что конкретные значения решения слау в ... 300 или 10000 уравнений никому не нужны. Нужен график этих значений, знаете, намного удобней их так смотреть :) Тем более моя задача состояла в СЛАУ из 10000 уравнений и последующем выводе в график средствами Mathlab, вообще говоря, самое то, чтобы перенаправить поток и сразу вывести график.
Так что, не надо критиковать моя вывод, пожалуйста :)

[hurtsy]

ajaxtpm в сообщении #577440 писал(а):

Тем более моя задача состояла в СЛАУ из 10000 уравнений и последующем выводе в график средствами Mathlab, вообще говоря, самое то, чтобы перенаправить поток и сразу вывести график.
Так что, не надо критиковать моя вывод, пожалуйста

Угу, именно так я должен понимать слова

ajaxtpm в сообщении #574383 писал(а):

Почти что вторую неделю я пытаюсь разобрать, что в моей программе не работает. В начале по заданию не использовал массивы вообще и считал всё по сложным штукам, но потом забил и упростил, стало полегче, но он не химичит всё равно ><

Кроме того что такое

ajaxtpm в сообщении #574383 писал(а):

простенький метод обратной матрицы

и как вам удалось из программы, которая

ajaxtpm в сообщении #574383 писал(а):

не химичит всё равно

сотворить такое, что не боитесь нарушить запрет форума на халяву. Меня интересуют эффективные методы исправления программ, а проходить по вашей ссылке курс уважаемого мною "Интуита" или "давать отлуп вашей программе" (копирайт дед Щукарь) не входит в мои цели. С уважением,

[ajaxtpm]

hurtsy, чувак, это две разных программы, я не исправлял, а писал по новой, без интуитов и интернета вообще)
Решения матлаба X=A/B, даже без очков видно, что он как бэ считает обратную матрицу.
У меня всё не доходили руки взяться и подумать над методом до конца, руки дошли, всё разъяснил, на выходе работающая программа, вывод которой можно заточить под свои цели, уж не буду же я писать все возможные формы вывода, лишь бы ко мне не придрался ни один hurtsy :)

С чего всё началось ?

[hurtsy]

ajaxtpm в сообщении #577817 писал(а):

на выходе работающая программа

(Оффтоп)

Надо отвечать на вопросы, если можете.

Мне приходилось брать работающую программу с тестом из фонда алгоритмов МГУ. Тест пролетал как песня. При других исходных данных программа работала но неправильно. Можно ли называть такое работающей программой. Разумеется нет. Возможно, что вы писали обе свои программы с нуля, хотя и есть сомнения. Больше я ничего по вашему топику писать не собираюсь. Вторая программа ничем не лучше первой и меня удивило, что вы в ней искали. Желаю дальнейших успехов в программировании. С уважением,

[Chifu]

Если бы обратную матрицу было бы легко вычислить, то в других методах решения СЛАУ не было бы нужды. Ну а если бы алгоритмы вели себя так же хорошо на миллионе уравнений как на десяти, то вообще бы никаких проблем не было и не существовало бы Матлаба. Суть Матлаба в том и заключается, что вы пишите обратную матрицу, а как он её вычисляет вы не знаете, вы даёте задание, а уж Матлаб сам решает что делать, чтобы это задание выполнить.

[ajaxtpm]

hurtsy, ерунду несёте, честное слово :)

hurtsy в сообщении #577941 писал(а):

Мне приходилось брать работающую программу с тестом из фонда алгоритмов МГУ. Тест пролетал как песня. При других исходных данных программа работала но неправильно. Можно ли называть такое работающей программой. Разумеется нет. Возможно, что вы писали обе свои программы с нуля, хотя и есть сомнения. Больше я ничего по вашему топику писать не собираюсь. Вторая программа ничем не лучше первой и меня удивило, что вы в ней искали. Желаю дальнейших успехов в программировании. С уважением,

Какое МГУ ? Оно то здесь причём, бы вы лучше чего конкретного сказали, а не воду лили.
У метода есть такое понятие, как "сжимающая матрица", если норма матрицы не меньше еденицы, о 100% точном решении не может быть и речи. А в рамках сжимающих матриц мой алгоритм работает как нужно.
Мне уже кажется, что вы настолько не любите халяву, как входящую, так и исходящую, то готовы уже без всякого обоснования захамить и закритиковать.

Вот честно, от вас я не услышал ни одного нормального замечания, абсолютно ни одного, вам так нравится говорить впустую ?

Chifu, Матлаб (не будем говорить о разных версиях, это глупо) при таком выражении "X=A/B" начинаети добросовестно искать обратную матрицу B, а потом не менее добросовестно умножать. Это факт :)
Конкретное решение системы зависит от вида матрицы А, будь она например сильно разрежена, или какая другая. Я говорю о подходящем численном решении :)