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

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

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

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

Добавлен: 04.04.2021

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

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

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

> RungeKutt:=proc(f,y_a,a,b,n)
# f — функция
# y_a — начальное условие задачи Коши
# a, b — границы
# n — количество точек разбиения промежутка [a,b]

local i,h,x,y,dy,K1,K2,K3,K4;
y[0]:=y_a;
h:=(b-a)/n;
for i from 0 to n do

# задаем конечное множество точек на отрезке [a,b]

x[i]:=a+h*i;

end do;
for i from 0 to n do

K1[i]:=h*f(x[i],y[i]);
K2[i]:=h*f(x[i]+h/2,y[i]+K1[i]/2);
K3[i]:=h*f(x[i]+h/2,y[i]+K2[i]/2);
K4[i]:=h*f(x[i]+h,y[i]+K3[i]);
dy[i]:=(1/6)*(K1[i]+2*K2[i]+2*K3[i]+K4[i]);
y[i+1]:=y[i]+dy[i];

end do;
for i from 0 to n do

print(x=x[i], y=y[i]);

end do;
end proc:

> f:=(x,y)->y+(1+x)*yˆ 2:
> dy/dx=f(x,y);
Проверим работу процедуры
> RungeKutt(f,-1,1,1.5,5);

x = 1., y = -1

x = 1.100000000, y = -0.9090933148
x = 1.200000000, y = -0.8333367499
x = 1.300000000, y = -0.7692344925
x = 1.400000000, y = -0.7142893912
x = 1.500000000, y = -0.6666701276

6.6.4.

Метод Адамса

> restart;
> Adams:=proc(f,y_a,a,b,n)

66


background image

# f — функция
# y_a — начальное условие задачи Коши
# a, b — границы
# n — количество точек разбиения промежутка [a,b]

local i,h,x,y,dy,K1,K2,K3,K4;
y[0]:=y_a;
h:=(b-a)/n;
for i from 0 to n do

# задаем конечное множество точек на отрезке [a,b]

x[i]:=a+h*i;

end do;
y[1]:=y[0]+h*f(x[0],y[0]);

# вычисление приближенных значений в первых трех точках мето-

дом Эйлера

y[2]:=y[1]+h*f(x[1],y[1]);
y[3]:=y[2]+h*f(x[2],y[2]);
for i from 3 to n do

# вычисление значений методом Адамса

y[i+1]:=y[i]+h/24*(55*f(x[i],y[i])-59*f(x[i-1],y[i-1])+37*f(x[i-2],y[i-2])-9*f(x[i-
3],y[i-3]));

end do;
for i from 0 to n do

print(x=x[i], y=y[i]);

end do;
end proc:

> f:=(x,y)->y+(1+x)*yˆ 2:
> dy/dx=f(x,y);
Проверим работу процедуры
> Adams(f,-1,1,1.5,5);

x = 1., y = -1

x = 1.100000000, y = -0.9000000000
x = 1.200000000, y = -0.8199000000
x = 1.300000000, y = -0.7539980778
x = 1.400000000, y = -0.7031569156
x = 1.500000000, y = -0.6568883609

67


background image

7.

Примеры программ на языке Fortran

7.1.

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

! Программа решения нелинейного уравнения методом

! деления отрезка пополам

program delenie
implicit none
real(4), external :: fun ! уравнение, корни которого нужно найти
real a,b,c,eps ! [a,b] - отрезок, на котором есть один корень,
! с - корень уравнения,
! eps - точность решения
integer i ! счетчик числа итераций
logical l ! вспомогательная переменная

write (*,*) "Введите a,b"
read (*,*) a,b
eps=0.0001
i=1
l=.true.
if (fun(a)*fun(b)<=0) then

do while ((b-a)>2*eps)

if (fun(a).eq.0) then

c=a
l=.false.
exit

end if
if (fun(b).eq.0) then

c=b

l=.false.
exit

end if
i=i+1
c=(a+b)/2
if (fun(c)*fun(a)>0) then

a=c

else

b=c

end if

end do
if (l) c=(a+b)/2

68


background image

write(*,*) "Корень=", c

write(*,*) "Число итераций=", i
write(*,*) "Проверка: f(x)=", fun(c)

else

write (*,*) "Корней на отрезке [",a,",",b,"] нет или более 1"

end if

end program

function fun(x)
real(4) :: fun, x

fun=x-2

end function fun

! Программа решения нелинейного уравнения методом хорд
program horda
implicit none
real(4), external :: fun ! уравнение, корни которого нужно найти
real a,b,eps, xo, xn ! [a,b] - отрезок, на котором есть один корень,
! eps - точность решения
! xo, xn - корень уравнения на предыдущей и текущей итерациях,
integer i ! счетчик числа итераций
logical l ! вспомогательная переменная

write (*,*) "Введите a,b"
read (*,*) a,b
if (fun(a)*fun(b)<=0) then

eps=0.0001
i=1
xo=a-(b-a)*fun(a)/(fun(b)-fun(a)) ! первая итерация
if (fun(a)*fun(xo)>0) then

a=xo

else

b=xo

end if
l=.true.
do while (l)

if (fun(a).eq.0) then

xn=a
exit

end if
if (fun(b).eq.0) then

xn=b

69


background image

exit

end if
i=i+1
xn=a-(b-a)*fun(a)/(fun(b)-fun(a))
if (fun(xn)*fun(a)>0) then

a=xn

else

b=xn

end if
l=abs(xn-xo)>eps
xo=xn

end do
write(*,*) "корень=", xn
write(*,*) "Число итераций=", i
write(*,*) "Проверка: f(x)=", fun(xn)

else

write (*,*) "Корней на отрезке [",a,",",b,"] нет или более 1"

end if

end program

function fun(x)
real(4) :: fun, x

fun=(x-2)**2-4*sin(x)

end function fun

! Программа решения нелинейного уравнения методом Ньютона
program newton
implicit none
real(4), external :: fun, funs ! функция и ее производная
real a,b,eps, xo, xn ! [a,b] - отрезок, на котором есть один корень,
! eps - точность решения
!xo, xn - корень уравнения на предыдущей и текущей итерациях,
integer i ! счетчик числа итераций
logical l ! вспомогательная переменная

write (*,*) "Введите x0"! начальное приближение
read (*,*) xo
eps=0.0001
i=0
l=.true.
do while (l)

if (fun(xo).eq.0) exit

70