VPF::Работа итерационного метода - Форум программистов

PandaRus

Дата 5.1.2008, 17:17 (ссылка)

(нет голосов)

Загрузка ...

Новичок

Профиль
Группа: Участник
Сообщений: 28
Регистрация: 27.12.2006

Репутация: нет
Всего: нет

Суть программы - решение большой СЛАУ. Первоначально использовал метод Гаусса, но при увеличении размерности СЛАУ (у меня она определяется матрицей G и B) в квадратичной зависимости рос размер памяти под переменные.

Т.к. матрица G сильно разреженная, начал использовать компактное хранение матрицы G и итерационный метод решения СЛАУ методом Гаусса-Зейделя. Реализацию данного метода взял из книжки, переведя с языка Фортран.

Но вот он то и дает сбой (в коде он начинается с комментария "Решение СЛАУ итерационным методом").

В начале программы я задаю параметры моделирования: h, s и r. Размерность СЛАУ = n, где n=h*s*r.
Если введу, например, h=5, s=5, r=5 - программа считает правильно,
h=8, s=9, r=5 - программа выводит значения всех переменных, равным -1.#I0
h=10, s=9, r=9 - программа вылетает.
То есть при превышении какого-то порога n происходит сбой.

Код


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

int main()
{
    int I,J,IAA,IAB,IT,IEND,N,S,H,L,i,j,m,k,s,r,n,ii,si,rj,h,hm,hmj,RS;
    double U,EPS,F;

    double Ps[16][16][31],Pr[16][16][31],Ph[16][16][31],Kg;
    double** G=new double*[2000];
    for(i=0;i<2000;i++)
        G[i]=new double[2000];

    double* Fi1=new double[2000];
    double* Fi2=new double[2000];
    double* B=new double[2000];
    double* S1=new double[2000];

    double* X=new double[2000];
    double* AN=new double[2000];
    double* AD=new double[2000];
    int* IA=new int[10000];
    int* JA=new int[10000];

    printf("\t     Programma rascheta toka\n\n\n");

    FILE *outfile;
    outfile=fopen("1.txt","w");

//Ввод параметров на моделирование-----------------------------
    h=9;                            //размеры тела
    s=5;                            
    r=9;
    k=s*r*h;            N=k;            

    for(i=1; i<=k; i++)                //места подключения электродов
        B[i]=0;
    B[1]=100000; B[9]=-100000;


    for(m=0; m<=h; m++)
        for(i=0; i<=s; i++)
            for(j=0; j<=r; j++)
            {
                Ps[m][i][j]=0;
                Pr[m][i][j]=0;
                Ph[m][i][j]=0;
            }

    for(m=1; m<=h; m++)                //задание характеристики тела
        for(i=1; i<=s; i++)
            for(j=1; j<=r-1; j++)
                Ps[m][i][j]=10;
    for(m=1; m<=h; m++)
        for(i=1; i<=s-1; i++)
            for(j=1; j<=r; j++)
                Pr[m][i][j]=10;
    for(m=1; m<=h-1; m++)
        for(i=1; i<=s; i++)
            for(j=1; j<=r; j++)
                Ph[m][i][j]=10;

//==============================================================
//==============================================================
//Формирование матрицы проводимостей G--------------------------
    printf("\t     Formirovanie matrici G\n");
    for(i=1; i<=k; i++)
        for(j=1; j<=k; j++)
            G[i][j]=0;

    RS=r*s;
    for(i=1; i<=k; i++)
        for(j=1; j<=k; j++)
        {
//Определение уровня, строки и столбца--------------------------
            if(i%RS==0)
                    hm=i/RS;
                else
                    if(i<RS)
                        hm=1;
                    else
                        hm=(i-i%(RS))/(RS)+1;
            if(j%RS==0)
                    hmj=j/RS;
                else
                    if(j<RS)
                        hmj=1;
                    else
                        hmj=(j-j%(RS))/(RS)+1;
            ii=i%(RS);
            if (ii==0)
                ii=RS;
            
            rj=ii%r;
            if(rj==0)
                rj=r;
                
            if(ii==rj)
                si=1;
            else
                si=(ii-rj)/r+1;
//--------------------------------------------------------------        
        
            if(i==j)            //Формирование главной диагонали    
                G[i][j]=Pr[hm][si-1][rj]+Pr[hm][si][rj]+
Ps[hm][si][rj-1]+Ps[hm][si][rj]+Ph[hm][si][rj]+Ph[hm-1][si][rj];
            //    AD[i]=Pr[hm][si-1][rj]+Pr[hm][si][rj]+
//Ps[hm][si][rj-1]+Ps[hm][si][rj]+Ph[hm][si][rj]+Ph[hm-1][si][rj];
            else
            {

                if(hm==hmj)
                {
                    if(j==i+r)
                        G[i][j]=-Pr[hm][si][rj];
                    if(j==i-r)
                        G[i][j]=-Pr[hm][si-1][rj];
                    if(j==i+1)
                        if(rj<=r-1)
                            G[i][j]=-Ps[hm][si][rj];
                    if(j==i-1)
                        if(rj>=2)
                            G[i][j]=-Ps[hm][si][rj-1];
                }
                if(j==i+RS)
                    G[i][j]=-Ph[hm][si][rj];
                if(j==i-RS)
                    G[i][j]=-Ph[hm-1][si][rj];
            }

        }

//    Формирование компактной формы хранения матрицы G



    for(I=1;I<=N;I++)
        AD[I]=G[I][I];
    S=0;
    H=0;
    for(I=1;I<=N;I++)
    {
        L=0;
        for(J=1;J<=N;J++)
            if(G[I][J]!=0 && I!=J)
            {
                L=1;
                S++;
                JA[S]=J;
                AN[S]=G[I][J];
                if(I>H)
                {
                    H++;
                    IA[H]=S;
                }
            }
        if(L==0)
        {
            H++;
            IA[H]=S+1;
        }
    }

    H++;
    IA[H]=S+1;


//    Решение СЛАУ итерационным методом
    
    F=1;

printf("\t    Rechenie\n");
EPS=0.1;
    for(I=1;I<=N;I++)
        X[I]=B[I]/AD[I];
    IT=0;
A:    IT++;
    IEND=0;
    for(I=1;I<=N;I++)
    {
        IAA=IA[I];
        IAB=IA[I+1]-1;
        if(IAB<IAA)
            goto B;
        U=B[I];
        for(J=IAA;J<=IAB;J++)
            U=U-AN[J]*X[JA[J]];
        U=U/AD[I]-X[I];
        if(fabs(U)>EPS)
            IEND=1;
        X[I]=X[I]+F*U;
    }
B:    if(IEND==1)
        goto A;
//------------------------------------------------------------------------------------------        

    for(I=1; I<=N; I++)
        {
                printf(" % 5.3f",X[I]);
                printf("\n\n");
        }

                printf("IT= %d    ",IT);
printf("U= % 10.9f",fabs(U));


    for(i=0;i<2000;i++)
    {
        delete [] G[i];
    }
    delete [] G;


    getch();
    return 0;
}

Kuvaldis

Дата 7.1.2008, 04:37 (ссылка)

(нет голосов)

Загрузка ...

механик-вредитель

Профиль
Группа: Участник Клуба
Сообщений: 1189
Регистрация: 16.6.2006
Где: Минск

Репутация: нет
Всего: 61

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

Я сам на лабах обжигался в свое время. В общем, советую посмотеть теорию для начала (гуд - самарский ) smile

Добавлено через 1 минуту и 1 секунду
можешь еще поставить long double
Скорее всего размерность считаемая увеличится, но тоже будет пороговое значение

--------------------

Помни - когда ты спишь, враг не дремлет
Спи чаще и дольше, изматывай врага бессоницей

0 Пользователей читают эту тему (0 Гостей и 0 Скрытых Пользователей)
0 Пользователей:
« Предыдущая тема \| C/C++: Для новичков \| Следующая тема »