Рунге-кутта для спаренных ОДУ

Осторожно: в формуле RK4 слишком много h:

k2 = h*dfuns{ [...] +(0.5*h*k1)]);
k3 = h*dfuns{ [...] +(0.5*h*k2]);

должно быть

k2 = h*dfuns{ [...] +(0.5*k1)]);
k3 = h*dfuns{ [...] +(0.5*k2]);

(последний h удален).

Однако это не имеет значения для приведенного вами примера, поскольку там h=1.

Но кроме этой маленькой ошибки, я не думаю, что вы на самом деле делаете что-то не так.

Если я построю график решения, созданного более продвинутым адаптивным RK порядка 4/5, реализованным в ode45:

F{1} = @(e,y) +0.6*y(3);
F{2} = @(e,y) -0.6*y(3) + 0.001407*y(4)*y(3);
F{3} = @(e,y)            -0.001407*y(4)*y(3);

tend  = 24;
steps = 24;
y0    = [0 5 995];
plotN = 2;

sol1 = coupled_ode(0,F, steps, 0,tend, y0, 'euler');
sol2 = coupled_ode(0,F, steps, 0,tend, y0, 'heun');
sol3 = coupled_ode(0,F, steps, 0,tend, y0, 'rk4');

figure(1), clf, hold on
plot(sol1(:,1), sol1(:,plotN+1),...
     sol2(:,1), sol2(:,plotN+1),...
     sol3(:,1), sol3(:,plotN+1));

% New solution, generated by ODE45
opts = odeset('AbsTol', 1e-12, 'RelTol', 1e-12);

fcn = @(t,y) [F{1}(0,[0; y])
              F{2}(0,[0; y])
              F{3}(0,[0; y])];
[t,solN] = ode45(fcn, [0 tend], y0, opts);    
plot(t, solN(:,plotN))

legend('Euler', 'Heun', 'RK4', 'ODE45');
xlabel('t');    

Тогда у нас есть что-то более правдоподобное для сравнения.

Теперь простой и понятный RK4 действительно ужасно работает в этом изолированном случае:

Однако, если я просто переверну знаки последнего члена в последних двух функциях:

%                       ± 
F{2} = @(e,y) +0.6*y(3) - 0.001407*y(4)*y(3);
F{3} = @(e,y)            +0.001407*y(4)*y(3);

Тогда получаем вот что:

Основная причина, по которой RK4 плохо работает в вашем случае, заключается в размере шага. Адаптивный RK4 / 5 (с допуском, установленным на 1 вместо 1e-12, как указано выше) дает среднее значение δt = 0,15. Это означает, что базовый анализ ошибок показал, что для этой конкретной проблемы h = 0.15 - это самый большой шаг, который вы можете предпринять, не допуская недопустимой ошибки.

Но вы принимали h = 1, что действительно дает большую накопленную ошибку.

Тот факт, что Хойн и Эйлер так хорошо работают в вашем случае, ну просто простая удача, как показывает приведенный выше пример с инверсией знаков.

Добро пожаловать в мир числовой математики - никогда не бывает 1 метода, который лучше всего подходил бы для всех задач при любых обстоятельствах :)

Помимо ошибки, описанной в более старом ответе, в реализации действительно есть фундаментальная методологическая ошибка. Во-первых, реализация верна для скалярных дифференциальных уравнений первого порядка. Но в тот момент, когда вы пытаетесь использовать его в связанной системе, разделенная обработка этапов в методе Рунге-Кутта (обратите внимание, что Хойн является просто копией шага Эйлера) сводит их к методу первого порядка.

В частности, начиная с

      k2=h*dfuns{j-1}(E, [sol(i-1,1)+(0.5*h), sol(i-1,2:end)+(0.5*h*k1)]);

добавление 0.5*k1 к sol(i-1,2:end) означает добавление вектора наклонов первой ступени, а не добавление одного и того же значения наклона ко всем компонентам вектора положения.

Учет этого приводит к изменению реализации

  function sol = coupled_ode(E, dfuns, steps, a, b, ini, method)
    range = b-a;
    h=range/steps;  
    rows = steps+1;
    columns = size(dfuns)(2)+1;
    sol= zeros(rows,columns);
    k = ones(4,columns);
    sol(1,1)=a;
    sol(1,2:end)=ini(1:end);
    for i=2:abs(rows)
      sol(i,1)=sol(i-1,1)+h;      
      if strcmp("euler",method)
        for j=2:columns
          sol(i,j)=sol(i-1,j)+h*dfuns{j-1}(E, sol(i-1,1:end));    
        end  
      elseif strcmp("heun",method)
        for j=2:columns
          k(1,j) = h*dfuns{j-1}(E, sol(i-1,1:end));
        end
        for j=2:columns
           sol(i,j)=sol(i-1,j)+h*dfuns{j-1}(E, sol(i-1,1:end)+k(1,1:end)); 
        end         
      elseif strcmp("rk4",method)
        for j=2:columns
          k(1,j)=h*dfuns{j-1}(E, sol(i-1,:));
        end
        for j=2:columns
          k(2,j)=h*dfuns{j-1}(E, sol(i-1,:)+0.5*k(1,:));
        end
        for j=2:columns
          k(3,j)=h*dfuns{j-1}(E, sol(i-1,:)+0.5*k(2,:));
        end
        for j=2:columns
          k(4,j)=h*dfuns{j-1}(E, sol(i-1,:)+k(3,:)); 
        end
        sol(i,2:end)=sol(i-1,2:end)+(1/6)*(k(1,2:end)+(2*k(2,2:end))+(2*k(3,2:end))+k(4,2:end));       
      end
    end
  end 

Как можно видеть, цикл по компонентам вектора часто повторяется. Это можно скрыть, используя полную векторизацию с использованием векторнозначной функции для правой стороны связанной системы ОДУ.

График для второй компоненты решения с этими изменениями дает гораздо более разумный график для размера шага 1.

и с разбиением на 120 интервалов для шага 0,2

где график для RK4 практически не изменился, в то время как два других двигались к нему снизу и сверху.