ВУЗ: Не указан
Категория: Не указан
Дисциплина: Не указана
Добавлен: 04.04.2021
Просмотров: 738
Скачиваний: 1
i=i+1
xn=xo-fun(xo)/funs(xo)
l=abs(xn-xo)>eps
xo=xn
end do
write(*,*) "Корень=", xo
write(*,*) "Число итераций=", i
write(*,*) "Проверка: f(x)=", fun(xo)
end program
function fun(x)
real(4) :: fun, x
fun=sin(x)+x**2
end function fun
function funs(x) !производная функции
real :: funs, x
funs=cos(x)+2*x
end function
7.2.
Интерполяция
! Программа строит интерполяционный многочлен Лагранжа Lagr(x)
! для функции f(x) на отрезке [a,b] по n+1 узлам интерполяции
program Lagrang
implicit none
real*8 a, b, h, x, chag
real*8, allocatable :: xmas(:), fun(:)
integer i, n, ch
! Вводим a и b - границы отрезка, на котором ищется многочлен
Лагранжа
write (*, "(A,
\
)") "Введите a "; read (*, *) a
write (*, "(A,
\
)") "Введите b "; read (*, *) b
! Вводим n - число отрезков, на которые разбивается
! интервал [a,b] (Всего будет n+1 узел интерполяции)
write (*, "(A,
\
)") "Введите n "; read (*, *) n
! Шаг интерполяции - возможно неравномерный
h=(b-a)/n
! Вводим массивы узлов интерполяции xmas
! и значения функции f(x) в этих узлах
71
! Если шаг неравномерный - ввод данных
! следует осуществлять по-другому
allocate (xmas(0:n), fun(0:n))
do i=0, n
xmas(i)=a+i*h
fun(i)=f(xmas(i))
end do
! Вывод данных в файл для построения графика в Gnuplot
! Число точек для построения графика ch
ch=200
chag=(b-a)/ch
x=a
! открываем файл graph.dat в текущей директории
open (1, file="graph.dat")
! записываем форматированные результаты в файл
do i=0, ch
write (1, 10) x, f(x), Lagr(x)
x=x+chag
end do
10 format (3(d12.4, 3x)) ! формат вывода данных
! освобождаем память под массивы
deallocate (xmas, fun)
! внутренние функции f(x) -задана и Lagr(x)-строим по формуле
contains
function f(x)
implicit none
real*8 f, x
f=sin(x)*cos(2*x)
end function
function Lagr(x)
implicit none
real*8 Lagr, x, s
integer i, j
do i=0, n
s=1
do j=0, n
if(i.ne.j) then
s=s*(x-xmas(j))/(xmas(i)-xmas(j))
end if
end do
72
Lagr=Lagr+fun(i)*s
end do
end function
end program Lagrang
! Программа строит интерполяционный многочлен Ньютона Newt(x)
! для функции f(x) на отрезке [a,b] по n+1 узлам интерполяции
program Newton
implicit none
real*8 a, b, h, x, chag
real*8, allocatable :: raznost(:), temp(:)
integer i,k,n,ch
! Вводим a и b - границы отрезка,
! на котором ищется многочлен Ньютона
write (*, "(A,
\
)") "Введите a "; read (*, *) a
write (*, "(A,
\
)") "Введите b "; read (*, *) b
! Вводим n - число отрезков, на которые разбивается
! интервал [a,b] (Всего будет n+1 узел интерполяции)
write (*, "(A,
\
)") "Введите n "; read (*, *) n
! Шаг интерполяции
h=(b-a)/n
! Вводим вспомогательный массив temp
allocate (temp(0:n))
do i=0, n
x=a+i*h
temp(i)=f(x) ! заносим в массив temp значения функции
! в узлах интерполяции
end do
! Вводим массив, где будем хранить только первые
! конечные разности всех порядков
allocate (raznost(n))
! заполняем массив raznost
do k=1, n
do i=0, n-k
temp(i)=temp(i+1)-temp(i)
end do
raznost(k)=temp(0) ! для дальнейших вычислений
! нужна только первая разность
end do
! Вывод данных в файл для построения графика в Gnuplot
! Число точек для построения графика ch
73
ch=200
chag=(b-a)/ch
x=a
! открываем файл graph.dat в текущей директории
open (1, file="graph.dat")
! записываем форматированные результаты в файл
do i=0, ch
write (1, 10) x, f(x), Newt(x)
x=x+chag
end do
10 format (3(d12.4, 3x)) ! формат вывода данных
! освобождаем память под массивы
deallocate (raznost, temp) ! внутренние функции f(x) -задана и
Newt(x)-строим по формуле
contains
function f(x)
implicit none
real*8 f, x
f=sin(x)*cos(2*x)*x
end function function Newt(x)
implicit none
real*8 Newt, x, prod, xk
Newt=f(a)
prod=1
xk=a
do k=1, n
prod=prod*(x-xk)/(h*k)
Newt=Newt+raznost(k)*prod
xk=xk+h
end do
end function
end program Newton
7.3.
Численное интегрирование
! Расчет определенного интеграла методом прямоугольников
program pramoug
implicit none
real*8 a, b, eps, I1, I2
74
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
! расчет интеграла в зависимости от числа разбиений по
! методу прямоугольников (по центру)
function integr(n)
real*8 integr, x, h
integer n, i
integr=0.0
h=(b-a)/n !длина элементарного интервала
x=a+h/2.0
do i=1, n
integr=integr+fun(x)*h
x=x+h
75