Решение системы дифференциальных уравнений методом Адамса-Бэшфорда 3 порядка - Fortra

[Turb_man]

Здравствуйте уважаемые форумчани.
Есть следующая задача: решить методом Адамса-Бэшфорда 3 порядка следующие уравнение
dv/dt = - (sx + v)/tau, где x=x(t), v=v(t) - координата и скорость движения частицы, t- время движения, s, tau - произвольные параметры.
Начальные условия x(0)=x0, v(0)=0.
Так как x = dv/dt, то получается уравнение x"+x/tau+s/tau*x=0.
Точное решение этого уравнения x = x0/sqrt(1-1/(4*s*tau))*exp(-(i-1)*dt/(2*tau))*sin(sqrt(s/tau-1/(4*tau**2)) * (i-1)*dt + atan(sqrt(4*s*tau-1))).
Вроде бы я точное решение нашёл верно.
Первые три значения по скорости я задаю с помощью метода Эйлера, а по координате с помощью точного решения( хотя можно с помощью метода Эйлера или Рунге-Кутта, но это не принципиально). Визуальна оба точное решение и численное решение полученное с помощью метода Адамса-Бэшфорда совпадают, но когдая хотел удостоверится и проверить, что действительно метод 3 порядка, то получилось, что при уменьшении шага ошибка между точным и численным решением меняется не пропорционально кубу шага по времени. Код я писал сам, скорее всего где в самом методе Адамса-Бэшфорда намудрил, хотя это и сложно, но возможно, и метод Эйлера в этом коде отлично работает.
Буду рад, если скажете где именно намудрил)))))

[Turb_man]

Вот код на Fortran
program interpolation
integer i
real*8 x0, tau, x(1000), s, dt, y(1000000,2), z(1000000,2)
character(10) file_name
c Запись точного решения в файл
file_name = 'myfile.dat'
c Задание начальных условий
x0 = 5.0
tau = 3.0
s = 5.0
dt = 0.02
open(unit = 1, file = file_name)
do i = 1, 1000
c Аналитичсекое решение уравнения
x(i) = x0/sqrt(1-1/(4*s*tau))*exp(-(i-1)*dt/(2*tau))*
&sin(sqrt(s/tau-1/(4*tau**2)) * (i-1)*dt + atan(sqrt(4*s*tau-1)))
enddo
do i = 1, 1000
write(1,*) (i-1)*dt, x(i)
enddo
close(1)
c Вызов функций интерполяции
! call eler(dt,z)
! call adams(dt,y)
call discrepancy
end
c--------------------------------------------------------------------------
c f = (u(x(t),t) - v)/tau
c Метод Адомса-Бэшфорта 3 порядка
subroutine adams(dt,y)
c Использование внешней функции
external F
real*8 k1(2), k2(2), k3(2), F
real*8 y(1000,2), dt, c1, c2, c3
integer i, j
character(10) file_name
c Имя файла для записи данных
file_name = 'adams3.dat'
c Шаг по времени
y = 0.0
c1 = 23.0/12.0
c2 = -16.0/12.0
c3 = 5.0/12.0
c Открытие файла для записи данных
open(unit = 2, file = file_name)
y(1,1) = 5.0
y(1,2) = 0.0
y(2,1) = 4.9983371232164568
y(2,2) = -0.16666666666666669
y(3,1) = 4.9933643379359083
y(3,2) = -0.33222222222222220
c Алгоритм Адомса-Бэшфорта 3 порядка
do i = 1, 997
k1(1) = y(i+2,2)* dt
k1(2) = F(y(i+2,1), y(i+2,2))*dt
k2(1) = y(i+1,2) * dt
k2(2) = F(y(i+1,1), y(i+1,2)) * dt
k3(1) = y(i,2) * dt
k3(2) = F(y(i,1), y(i,2)) * dt
do j = 1,2
y(i+3,j)=y(i+2,j) + (23.0*k1(j) - 16.0*k2(j) + 5.0*k3(j))/12.0
enddo
enddo
c Запись данных в файл adams3
do i = 1, 1000
write(2,*) (i-1)*dt, y(i,1)
enddo
close(2)
end
c-----------------------------------------------------------------------
function F(x,v)
real*8 x,v,f,s,tau
tau = 3.0
s = 5.0
f = 0.0
f = (-s*x-v)/tau
return
end
c-----------------------------------------------------------------------------
c Эйлера
subroutine eler(dt,z)
c Использование внешней функции
external F
real*8 dt, F
real*8 z(1000,2)
integer i
character(15) file_name
c Имя файла для записи данных
file_name = 'eler_file.dat'
c dt = 0.02
z(1,1) = 5.0
z(1,2) = 0.0
c Открытие файла для записи данных
open(unit = 1, file = file_name)
c Алгоритм Эйлера
do i = 1, 999
z(i+1,1) = z(i,1) + z(i,2) * dt
z(i+1,2) = z(i,2) + F(z(i,1), z(i,2)) * dt
enddo
c Запись данных в файл eler_file
do i = 1, 1000
write(1,*) (i-1)*dt, z(i,1), z(i,2)
enddo
close(1)
end
c-----------------------------------------------------------------------
! Проверка порядка сходимости метода, вычисление ошибки метода
subroutine discrepancy
integer i, j
real*8 x0, tau, x(1000000), s, dt, y(1000000,2), ead(10)
real*8 z(1000000,2), eel(10)
character(15) ff_name
ff_name = 'errofile.dat'
x0 = 5.0
tau = 3.0
s = 5.0
dt = 0.02
x(1) = x0
do j = 1, 10
dt = dt / 2
do i = 1, 999996
x(i) = x0/sqrt(1-1/(4*s*tau))*exp(-(i-1)*dt/(2*tau))*
&sin(sqrt(s/tau-1/(4*tau**2)) * (i-1)*dt + atan(sqrt(4*s*tau-1)))
enddo
call adams(dt,y)
call eler(dt,z)
ead(j) = abs(y(20*2**j,1) - x(20*2**j))
write(*,*) y(20*2**j,1), x(20*2**j), 20*2**j*dt
eel(j) = abs(z(20*2**j,1) - x(20*2**j))
enddo
dt = 0.02
open(unit = 5, file = ff_name)
do i = 1, 10
write(5,*) dt/(2**i), ead(i), eel(i)
enddo
close(5)
end
c-----------------------------------------------------------------------