Рунге Кутта в Фортране

Question

Я пытаюсь реализовать метод Рунге Кутты в Фортране и сталкиваюсь с проблемой сходимости. Я не знаю, сколько кода мне нужно показать, поэтому я опишу проблему подробно и, пожалуйста, подскажите, что мне следует добавить / удалить в / из поста, чтобы сделать его ответственным.

У меня есть 6-мерный вектор положения и скорости шара, и соответствующая система различий. землетрясения. которые описывают уравнения движения, из которых я хочу вычислить траекторию движения мяча, и сравнивать результаты для разных порядков метода РК.

Давайте сосредоточимся на 3-м порядке РК. Модель, которую я использую, реализована следующим образом:

k1 = h * f(vec_old,omega,phi)
k2 = h * f(vec_old + 0.5d0 * k1,omega,phi)
k3 = h * f(vec_old + 2d0 * k2 - k1,omega,phi)
vec = vec_old + (k1 + 4d0 * k2 + k3) / 6d0

куда f это функция, которая составляет уравнения движения (или эквивалентно RHS моей системы diff. eqs). Обратите внимание, что f не зависит от времени, поэтому имеет только 1 аргумент. h берет на себя роль небольшого временного шага.

Если мы хотим рассчитать траекторию мяча за конечное время total_timeи учесть общую ошибку epsilon, тогда мы должны убедиться, что каждый шаг принимает пропорциональную долю ошибки. Для первого шага я сделал следующее:

vec1 = solve(3,vec_old,h,omega,phi)
vec2 = solve(3,vec_old,h/2d0,omega,phi)
do while (maxval((/(abs(vec1(i) - vec2(i)),i=1,6)/)) > eps * h / (tot_time - current_time))
    h = h / 2d0
    vec1 = solve(3,vec_old,h,omega,phi)
    vec2 = solve(3,vec_old,h/2d0,omega,phi)
end do
vec = (8d0/7d0) * vec2 - (1d0/7d0) * vec1

куда solve(3,vec_old,h,omega,phi) является функцией, которая вычисляет один шаг RK, описанный выше. 3 обозначает порядок РК, который мы используем, vec_old текущее состояние вектора положения-скорости, h, h/2d0 оба представляют используемый временной шаг, и omega,phi только некоторые дополнительные параметры для f, Наконец, для первого шага мы устанавливаем current_time = 0d0,

Дело в том, что если мы используем RK 3-го порядка, у нас должна быть ошибка в $O(h^3)$, и поэтому она падает быстрее, чем линейно в h, Поэтому мы должны ожидать, что цикл while будет в конечном итоге остановлен для достаточно малых h,

Моя проблема в том, что цикл не сходится и даже не закрывается - соотношение

maxval(...) / eps * (...)

остается довольно постоянным, вплоть до eps * h / (tot_time - current_time)) становится нулевым из-за конечной точности.

Для полноты, это мое определение f:

function f(vec_old,omega,phi) result(vec)
    real(8),intent(in) :: vec_old(6),omega,phi
    real(8) :: vec(6)
    real(8) :: v,Fv

    v = sqrt(vec_old(4)**2+vec_old(5)**2+vec_old(6)**2)
    Fv = 0.0039d0 + 0.0058d0 / (1d0 + exp((v-35d0)/5d0))
    vec(1) = vec_old(4)
    vec(2) = vec_old(5)
    vec(3) = vec_old(6)
    vec(4) = -Fv * v * vec_old(4) + 4.1d-4 * omega * (vec_old(6)*sin(phi) - vec_old(5)*cos(phi))
    vec(5) = -Fv * v * vec_old(5) + 4.1d-4 * omega * vec_old(4)*cos(phi)
    vec(6) = -Fv * v * vec_old(6) - 4.1d-4 * omega * vec_old(4)*sin(phi) - 9.8d0
end function f

Кто-нибудь имеет представление о том, почему цикл while не сходится? Если что-то еще понадобится (вывод, другие фрагменты кода и т. Д.), Пожалуйста, скажите мне, и я добавлю это. Кроме того, если требуется обрезка, я обрежу все, что будет сочтено ненужным. Спасибо!

Answer 1

Чтобы вычислить ошибку шага, используя метод полшага, вам нужно вычислить аппроксимацию в t+h в обоих случаях, что означает два шага с размером шага h/2, Как это сейчас вы сравниваете приближение в t+h в приближении при t+h/2 который дает вам ошибку размера f(vec(t+h/2))*h/2,

Таким образом, перейдите на 3-х этапную процедуру

vec1 = solve(3,vec_old,h,omega,phi)
vec2 = solve(3,vec_old,h/2d0,omega,phi)
vec2 = solve(3,vec2   ,h/2d0,omega,phi)

в обоих местах разница vec2-vec1 должен быть порядка h^4,