Файл: Вычислительный эксперимент и методы вычислений.pdf

ВУЗ: Не указан

Категория: Не указан

Дисциплина: Не указана

Добавлен: 04.04.2021

Просмотров: 770

Скачиваний: 1

ВНИМАНИЕ! Если данный файл нарушает Ваши авторские права, то обязательно сообщите нам.
background image

end do

end function

function fun(x)! подынтегральная функция
real*8 fun, x

fun=sin(x)*exp(x)

end function
end ! Расчет определенного интеграла методом трапеций
program trapec
implicit none
real*8 a, b, eps, I1, I2
integer n, dn
logical pr

! ввод a,b - отрезка интегрирования
write(*, "(A,

\

)") "Введите a "; read(*, *) a

write(*, "(A,

\

)") "Введите b "; read(*, *) b

! ввод точности вычислений
write(*, "(A,

\

)") "Введите eps "; read(*, *) eps

! ввод начального числа частей,
! на которые разбивается отрезок [a,b]
write(*, "(A,

\

)") "Введите n0 (>1) "; read(*, *) n

! при каждом вычислении число разбиений
! увеличивается на 50% от исходного
dn=0.5*n
I1=integr(n) ! Значение интеграла с числом разбиений n
write (*, *) "I(",n, ")=", I1
pr=.true.
do while(pr)

n=n+dn ! увеличиваем n
I2=integr(n) ! вычисляем интеграл при увеличенном n
pr=(abs(I2-I1)>eps)
! если разность между двумя
! последовательными вычислениями меньше точности,
! выходим из цикла
write (*, *) "I(",n, ")=", I2
I1=I2

end do
write (*, *) "Integral= ", I2, "n= ", n, "eps= ", eps
! вывод результата на экран

contains

76


background image

! расчет интеграла в зависимости от числа разбиений
! по методу трапеций
function integr(n)
real*8 integr, x, h
integer n, i

integr=0.0
x=a
h=(b-a)/n !длина элементарного интервала
do i=1, n

integr=integr+(fun(x)+fun(x+h))/2.0*h
x=x+h

end do

end function

! подынтегральная функция
function fun(x)
real*8 fun, x

fun=sin(x)

end function
end

! Расчет определенного интеграла методом Симпсона
program simpson
implicit none
real*8 a, b, eps, I1, I2
integer n, dn
logical pr

! ввод a,b - отрезка интегрирования
write(*, "(A,

\

)") "Введите a "; read(*, *) a

write(*, "(A,

\

)") "Введите b "; read(*, *) b

! ввод точности вычислений
write(*, "(A,

\

)") "Введите eps "; read(*, *) eps

! ввод начального числа частей (четное),
! на которые разбивается отрезок [a,b]
pr=.true.
do while(pr)

write(*, "(A,

\

)") "Введите n0-четное "; read(*, *) n

if (mod(n,2)==1) then

write (*,*) "n - нечетное"

else

pr=.false.

77


background image

end if

end do
! при каждом вычислении число разбиений
! увеличивается на 50% от исходного
dn=0.5*n
I1=integr(n) ! Значение интеграла с числом разбиений n
write (*, *) "I(",n, ")=", I1

pr=.true.

do while(pr)

n=n+dn ! увеличиваем n
! Если n оказывается нечетным, добавляем к нему 1
if (mod(n,2)==1) n=n+1
I2=integr(n) ! вычисляем интеграл при увеличенном n
pr=(abs(I2-I1)>eps) ! если разность между двумя
! последовательными вычислениями меньше точности,
! выходим из цикла
write (*, *) "I(",n, ")=", I2
I1=I2

end do
write (*, *) "Интеграл= ", I2, "n= ", n, "eps= ", eps
! вывод результата на экран

contains

! расчет интеграла в зависимости от числа разбиений
! по методу Симпсона
function integr(n)
real*8 integr, x, h
integer n, i

integr=0.0
h=(b-a)/n !длина элементарного интервала
integr=fun(a)+fun(b)
do i=1, n-1, 2

x=a+i*h
integr=integr+4*fun(x)

end do
do i=2, n-2, 2

x=a+i*h
integr=integr+2*fun(x)

end do
integr=integr*h/3.0

78


background image

end function

! подынтегральная функция
function fun(x)
real*8 fun, x

fun=sin(x)*exp(x)

end function
end
! Расчет определенного интеграла методом Гаусса
program gauss_integr
implicit none
real*8 a, b, eps, I1, I2
integer n, dn
logical pr

! ввод a,b - отрезка интегрирования
write(*, "(A,

\

)") "Введите a "; read(*, *) a

write(*, "(A,

\

)") "Введите b "; read(*, *) b

! ввод точности вычислений
write(*, "(A,

\

)") "Введите eps "; read(*, *) eps

! ввод начального числа частей,
! на которые разбивается отрезок [a,b]
write(*, "(A,

\

)") "Введите n0 (>1) "; read(*, *) n

! при каждом вычислении число разбиений
! увеличивается на 50% от исходного
dn=0.5*n
I1=integr(n) ! Значение интеграла с числом разбиений n
write (*, *) "I(",n, ")=", I1
pr=.true.
do while(pr)

n=n+dn ! увеличиваем n
I2=integr(n) ! вычисляем интеграл при увеличенном n
pr=(abs(I2-I1)>eps) ! если разность между двумя
! последовательными вычислениями меньше точности,
! выходим из цикла
write (*, *) "I(",n, ")=", I2
I1=I2

end do
write (*, *) "Интеграл= ", I2, "n= ", n, "eps= ", eps
! вывод результата на экран

contains

79


background image

! расчет интеграла по методу Гаусса
! на каждом элементарном отрезке применяется
! 4-х точечный метод Гаусса
function integr(n)
real*8 integr, x, h, c(4), t(4)
integer n, i, j

! Массив t
t=(/-0.8611363115940492, -0.3399810435848646, &

0.3399810435848646, 0.8611363115940492/)

! Массив c
c=(/2*0.1739274225687284, 2*0.3260725774312716, &

2*0.3260725774312716, 2*0.1739274225687284/)

integr=0.0
x=a
h=(b-a)/n !длина элементарного интервала
do i=1, n

do j=1, 4

integr=integr+c(j)*fun((2*x+h)/2.0+h/2.0*t(j))

end do
x=x+h

end do
integr=integr*h/2.0

end function

! подынтегральная функция
function fun(x)
real*8 fun, x

fun=sin(x)*exp(x)
end function

end

7.4.

Решение систем линейных уравнений

!

решение

систем

линейных

алгебраических

уравнений

методом Гаусса

program gauss_ls
implicit none
integer, parameter :: n=4 !размер системы
integer i,j,k
real*8 a(n,n), b(n), apr(n,n), bpr(n), x(n), s

80