Форум программистов [Powered by Invision Power Board]

Код

! Решить урмат: U_t = U_xx + (exp(x^3)-1) в "стакане"
! U(0,x) = 0,
! U(t,0) = U(t,1)=0
! используюя следующий шаблон
!        o
!        |
!      o-o-o, схема явная, устойчива при соотношении 
!шагов h_t~h_x^2.
! Точность eps = 1e-4;

program Ex_2
implicit none;

! Размеры сетки по x и t - Mx подобран минимально путем нескольких запусков. 
! Динамически решить задачу не удалось, т.к. сетка очень большая и 
! превышает размер стека
integer(4), parameter :: Mx=100;
integer(4), parameter :: Mt=5*Mx*Mx;

! Точность
real(8), parameter :: eps = 0.0001_8; ! Что именно 8 байт в типе.

! Дескриптор файла вывода 
integer(4), parameter :: f_out = 1001;
! Счетчики
integer(4) :: m,n;
! Сетки для пересчета и невязка
real(8) , allocatable, dimension(:,:) :: U1, U2, R;
! И переменная максимума невязки
real(8) :: max_R;




allocate(U1(0:Mt,0:Mx), U2(0:2*Mt,0:2*Mx), R(1:Mt,1:Mx-1));

! Первый раз вычисляем ДУ
U1 = DU(Mt,Mx);

! Проверим точность
U2 = DU(2*Mt, 2*Mx);
    
forall( m = 1:Mx-1, n = 1:Mt) 
    R(n,m) = abs(U2(2*n,2*m) - U1(n,m));
end forall;
    
! Т.к. схема 1-го порядка,то 
! max(U2 - U1)<(2*r - 1)*eps == eps
max_R = maxval(R);
print *,' max(R) = ',max_R;
if( max_R < eps) then
    print *,'Answer in file Result.csv';
    
    ! Вычислили, теперь выводим результат
    open(f_out,file='Result.csv',status='replace');
    do n = 0,Mt
        write(f_out,"( <2*Mx+1>(f15.8,';') )") U2(n,:);
    end do;
    close(f_out,status='keep');
else
    print *,'Grid so small';
end if;


deallocate(U1, U2, R);

write(*,*) 'press enter';
read(*,*);
contains
! Конечно-разностная реализация
function DU(Mt,Mx) result(U)
    ! На вход 
    ! Размеры сетки по t и x
    integer(4), intent(in) :: Mt, Mx;
    
    ! На выход
    ! Сама сетка - в строке - t, в столбце x!
    real(8), dimension(0:Mt,0:Mx) :: U;
    
    ! Внутренние переменные --------------------------------
    ! Вспом. сетка для f(x) = exp(x**3+1)
    real(8), dimension(0:Mx) :: f;
    ! Счетчики
    integer(4) :: m,n;
    ! Шаги
    real(8) :: ht; 
    real(8) :: hx; 
    
    
    !------------ Тело функции ------------------------------
    ! Вычисляем шаги
    ht = (1.0-0.0)/dble(Mt);
    hx = (1.0-0.0)/dble(Mx);
    
    ! К О Д      П Р О Г Р А М М Ы
    U = dble(0.0); 
    ! Т.к. На границах у нас везде нули, то мы задали ГУ.
    ! Вспом. функция, чтобы не вычислять постоянно
    forall(m=0:Mx) f(m) = exp((m*hx)**3 - 1.0);

    !write(*,*) 'Size of grid (T,X) = (',N ,',', M, ')';

    do n = 0,Mt-1
        do m = 1,Mx-1 ! Только по внутренним точкам сетки, чтобы не портить ГУ
            if(m==0 .or. m==Mx) print *,'Error!!!';
            
            U(n+1,m) = ht/(hx**2)*( U(n,m+1) - 2*U(n,m) + U(n,m-1) ) + &
                   ht*f(m) + U(n,m);
        end do; ! m
    end do ! n
end function DU;
end program Ex_2;