Метки

, , , , , ,

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

\dot y = -yt,\quad y(0) = 1

с помощью интегрирования 4-м и 5-м порядками метода.Сначала запишем функцию f=-yt, возвращающую значение производной и сохраним её в файле f.m:
function deriv = f(t, y)

deriv = - y * t;

Далее сам метод Рунге-Кутта, выполняемый в функции runge (файл runge.m):

function y = runge(t0, y0, dt, h, tol)

flag = 1;
t = t0;
y = y0;
while (t < t0 + dt)
    k1 = h * f(t, y);
    k2 = h * f(t + 1/4*h,   y + 1/4*k1);
    k3 = h * f(t + 3/8*h,   y + 3/32*k1 +       9/32*k2);
    k4 = h * f(t + 12/13*h, y + 1932/2197*k1 -  7200/2197*k2 +  7296/2197*k3);
    k5 = h * f(t + h,       y + 439/216*k1 -    8*k2 +          3680/513*k3 - 845/4104*k4);
    k6 = h * f(t + h/2,     y - 8/27*k1 +   2*k2 - 3544/2565*k3 + 1859/4104*k4 - 11/40*k5);
    y = y + 25/216*k1 + 1408/2565*k3 + 2197/4104*k4 - k5/5;
    if (flag)
        z = y + 16/135*k1 + 6656/12825*k3 + 28561/56430*k4 - 9/50*k5 + 2/55*k6;
    end
    if (flag && abs(y - z) > tol)
        t = t0;
        y = y0;
        h = h * (tol * h/2/abs(z-y))^0.25;
    else
        t = t + h;
        flag = 0;
    end
end

Здесь t0 и y0 являются значениями функции в текущем узле сетки, dt — шаг по аргументу к следующему узлу сетки, h — начальный шаг интегрирования, tol — требуемая точность. Алгоритм разбивает интервал интегрирования от t0 до t0+dt на подсетку с шагом h. Пока булевая переменная flag содержит значение «Истина» (1), то мы проверяем, достигнута ли требуемая точность, как только она достигнута, выполнение операций на провекру точности прекращается. На каждой итерации алгоритма вычисляется значение функции в шести точках. Далее в y заносится значение в следующем узле подсетки, на которую разбил интервал интегрирования диффура метод Рунге-Кутта 4-го порядка, а в z — 5-го. Если заданная точность еще не достигнута, тоесть различие при использовани 4-го и 5-го порядков метода больше чем tol, шаг уменьшается. Если заданная точность уже достигнута, то просто идем по подсетке.

И наконец, скрипт, непосредственно выполняющий интегрирование. Идем по интервалу от нуля до трех с шагом, равным 0,1, потом проверяем себя встроенной функцией матлаба и рисуем график:

st=0;
en=3;
x=st:0.1:en;

y = zeros(1,length(x));
y(1) = 1;
for i=1:(length(x)-1)
    y(i+1) = runge(x(i), y(i), 0.1, 0.05, 1e-5);
end

[t,z] = ode45(@f, [st,en], 1);
plot(x,y,t,z);
plottools;
grid on;
xlabel(t);
ylabel(y);
title('y' = -yt');
legend('calc', 'ode');

Получили вот такое:runge
Синим цветом рисовалось интегрирование с помощью вышеприведенных файлов, а зеленым проверка функцией matlab’а ode45, делающей то же самое 🙂

При написании использовался материал статьи http://math.fullerton.edu/mathews/n2003/RungeKuttaFehlbergMod.html

Реклама