VPF::решение трёхдиагональной матрицы методом прогонки

Muzencia

Дата 5.5.2008, 14:48 (ссылка)

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

Загрузка ...

Новичок

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

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

плиз, у кого есть прога или кто знает, где её можно достать, поделись smile

v1rtu0z

Дата 7.5.2008, 09:34 (ссылка)

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

Загрузка ...

Шустрый

Профиль
Группа: Участник
Сообщений: 61
Регистрация: 21.4.2008
Где: Воронеж

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

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

popovda

Дата 2.6.2008, 15:50 (ссылка)

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

Загрузка ...

Опытный

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

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

Вот рабочая программа решения ОДУ прогонкой. Там прогоночные коэффициенты из дифура берутся. А так - всё то же самое

Код


! Решение ОДУ прогонкой
program Progonka

implicit none;

! Точность
real(8), parameter :: eps = 0.00005;
! Мин. и макс. аргументы
real(8), parameter :: x_min = 0.0, x_max = 0.5;

! Кол-во точек
integer(4), parameter :: N = 11000, N2 = 2*N;

! Переменные
! Параметр задачи
real(8) :: p_alpha;
! Шаг
real(8) :: h;
! Счетчик
integer(4) :: i;

! Файл вывода
integer(4),parameter :: id = 1000;
character(80) :: f_name;
logical :: RungeRule;

! Прогоночные коэффициенты
real(8), allocatable, dimension(:) :: y, alpha, beta; 
real(8), allocatable, dimension(:) :: y2, alpha2, beta2; 
! Вспом. переменные
real(8) :: A, B, C, mu1, mu2, xi1, xi2;


!=======================================================================
! С А М А    П Р О Г Р А М М А
p_alpha = 0.5;

main_loop: do

! Вычисляем шаг 
h = (x_max - x_min)/dble(N);
! Выделяем память под массивы
allocate( y(0:N), alpha(1:N), beta(1:N) );
! Начальные Условия
xi1 = 1.0;
mu1 = -h;

alpha(1) = xi1;
beta(1) = mu1;

write(*,*)
write(*,"(a10,F10.5)") 'Alpha = ', p_alpha;
write(*,"(a10,F10.5)") 'h = ',h;

write(*,"(a10,F10.5)") 'xi1 = ',xi1;
write(*,"(a10,F10.5)") 'mu1 = ',mu1;

! Прямой ход 
do i = 1,N-1

    A = 1.0 + (0.5*h*(p_alpha - 1.0)/(i*h + 1));
    B = 1.0 - (0.5*h*(p_alpha - 1.0)/(i*h + 1));
    C = 2.0 + p_alpha*h*h;
    
    ! Находим прогоночные коэффициенты
    alpha(i+1) = B/( C - A*alpha(i) );
    beta(i+1) = ( A*beta(i) + F(i) )/( C - A*alpha(i) );
end do;


! Обратный ход
mu2 = 2*h*((1.5)**p_alpha)/(p_alpha*h + 1.5);
xi2 = 1.5/(p_alpha*h + 1.5);

write(*,"(a10,F10.5)") 'xi2 = ',xi2;
write(*,"(a10,F10.5)") 'mu2 = ',mu2;

y(N) = (mu2 + xi2*beta(N))/(1 - xi2*alpha(N));

do i = N-1,0,-1
    y(i) = alpha(i+1)*y(i+1) + beta(i+1);
end do;

! Повторяем прогонку с уменьшенным в 2 раза шагом
! Вычисляем шаг 
h = h/2.0;
! Выделяем память под массивы
allocate( y2(0:N2), alpha2(1:N2), beta2(1:N2) );
! Начальные Условия
xi1 = 1.0;
mu1 = -h;

alpha2(1) = xi1;
beta2(1) = mu1;

! Прямой ход 
do i = 1,N2-1

    A = 1.0 + (0.5*h*(p_alpha - 1.0)/(i*h + 1));
    B = 1.0 - (0.5*h*(p_alpha - 1.0)/(i*h + 1));
    C = 2.0 + p_alpha*h*h;
    
    ! Находим прогоночные коэффициенты
    alpha2(i+1) = B/( C - A*alpha2(i) );
    beta2(i+1) = ( A*beta2(i) + F(i) )/( C - A*alpha2(i) );
end do;

! Обратный ход
mu2 = 2*h*((1.5)**p_alpha)/(p_alpha*h + 1.5);
xi2 = 1.5/(p_alpha*h + 1.5);

y2(N2) = (mu2 + xi2*beta2(N2))/(1 - xi2*alpha2(N2));

do i = N2-1,0,-1
    y2(i) = alpha2(i+1)*y2(i+1) + beta2(i+1);
end do;

! Проверка Правила Рунге
RungeRule = .true.;
do i = 0,N
    if( abs(y2(2*i) - y(i)) >= eps) then
            write(*,*) 'Runge Rule is wrong';
            RungeRule = .false.;
            exit;
    end if;
end do;

if( RungeRule) then

! Выводим в файл
write(f_name, "('Out_Alpha=',F3.1,'.csv')") p_alpha;

open(id, file = f_name, status = 'replace');
    write(id,*) 'i; x; y; alpha; beta';    
    
    write(id, "(1x,i4,';',F15.7,';',F15.7,';;')") 0, 0.0, y2(0);
    
    do i = 1,N2
        write(id, "(1x, i4, ';' ,F15.7, ';' ,F15.7, ';' ,F15.7 ,';', F15.7)") i, i*h, y2(i), alpha2(i), beta2(i);
    end do;
close(id, status = 'keep');

end if;

! Освобождаем память
deallocate( y, alpha, beta );
deallocate( y2, alpha2, beta2 );

! Меняем параметр alpha
p_alpha = p_alpha + 0.5;

if(p_alpha > 5.0) exit main_loop;

end do main_loop;

write(*,*) 'Press enter';
read(*,*);

contains

! Правая часть
real(8) function F(i)
    ! На входе
    integer(4), intent(in) :: i;
    
    F = (1.0 + i*h)**p_alpha * h*h;
end function F;


end program Progonka;

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

С уважением, Попов Д.А.

rrrFer

Дата 2.6.2008, 16:25 (ссылка)

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

Загрузка ...

Бывалый

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

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

Muzencia,
писал, когда выч.мат. изучал:

Код


 /*----------------------INCLUDES---------------------*/
#include <stdio.h>
#include <conio.h>
#include <string.h>
#include <stdlib.h>
/*-------------------procedures----------------------*/
void output(float A[5][5],int tn,int tm)
{
    int i,j;
    for(i=0;i<tn;i++)
    {
        printf("|");
        for(j=0;j<tm;j++)
        {
            if(A[i][j]>=0)
                printf(" ");
            if(A[i][j]<0)
                printf("");
            printf("%1.1f |",A[i][j]);
        }
        printf("\n");
    }
}    
void output(float A[5],int tn)
{
    int i;
    for(i=0;i<tn;i++)
    {
        printf("|");
        if(A[i]>=0)
            printf(" ");
        if(A[i]<0)
            printf("");
        printf("%1.3f |",A[i]);
    }
}
void main()
{
    float A[5][5]=
    {{1,-1,0,0,0},{2,10,-5,0,0},{0,-3,8,2,0},{0,0,1,-12,-7},{0,0,0,-5./6,1}},
    Y[5]={4,109,-30,-48,25./6},Z[5],B[5],a[5],K[5],V[5],C[5];
    int n=5,i;
    for(i=0;i<n;i++){
        Z[i]=0;C[i]=0;B[i]=0;a[i]=0;V[i]=0;K[i]=0;
    }
    K[0]=-A[0][1];
    K[n-1]=-A[n-1][n-2];
    V[0]=Y[0];
    V[n-1]=Y[n-1];
    for(i=0;i<n;i++){
        B[i]=A[i][i+1];
        a[i]=A[i][i-1];
        C[i]=-A[i][i];
    }
    for(i=1;i<n-1;i++){
        K[i]=B[i]/(C[i]-(a[i]*K[i-1]));
        V[i]=((a[i]*V[i-1])-Y[i])/(C[i]-(a[i]*K[i-1]));
    }
    Z[n-1]=(V[n-1]+K[n-1]*V[n-2])/(1-(K[n-1]*K[n-2]));
    for(i=n-2;i>=0;i--)
        Z[i]=K[i]*Z[i+1]+V[i];
    float s=0;
    printf("\n\nFind errors^\n");
    output(Z,n);
    for(int j=0;j<n;printf("\n%d: %3.5f",j+1,s-Y[j]),s=0,j++)
        for(i=0;i<n;i++)
            s+=A[j][i]*Z[i];
    getch();
}

не гаранитирую, что правильно работает

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