Matlab: ode45 и метод Рунге-Кутты 4-го порядка дают разные значения - PullRequest
Новая
       6

Matlab: ode45 и метод Рунге-Кутты 4-го порядка дают разные значения

0 голосов
/ 28 августа 2018

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

kuramoto= @(x,K,N,Omega)Omega+(K/N)*sum(sin(x*ones(1,N)-(ones(N,1)*x')))'
%Kuramoto is a model of N coupled ocilators (such as multiple radiowaves)
%The solution to the model is the phase of each ocilator
%[Kuramoto Equation][1]

theta(:,1) = 2*pi*randn(N,1);
t0 = theta(:,1);
[t,y] = ode45(@(t,y)kuramoto(theta(:,1),K,N,omega),tspan,t0);

%Runge-Kutta method
for j=1:iter
k1=kuramoto(theta(:,j),K,N,omega);
k2=kuramoto(theta(:,j)+0.5*h*k1,K,N,omega);
k3=kuramoto(theta(:,j)+0.5*h*k2,K,N,omega);           
k4=kuramoto(theta(:,j)+h*k3,K,N,omega);
theta(:, j+1)=theta(:,j)+(h/6)*(k1+2*k2+2*k3+k4);
end

Оба метода выводят матрицу с N строками (где каждая строка представляет отдельный осциллятор) и M столбцами (где M представляет решение в данный момент времени). У меня есть ode45, обеспечивающий решения от 0 до 0,5 с интервалом 0,1. Чтобы сравнить методы, я вычитаю матрицу, полученную из Рунге-Кутты, с матрицей, полученной с помощью ode45. В идеале, оба должны иметь одинаковые значения, и результат должен быть нулевой матрицей, но вместо этого я получаю такие значения, как: 

0   -0.0003   -0.0012   -0.0027   -0.0048   -0.0076
0    0.0003    0.0012    0.0027    0.0048    0.0076
%here I have only two oscillators from t = [0.0,0.5]

Существует небольшая разница между двумя матрицами (которая увеличивается с большими интервалами времени). Но необычно, что общее значение, рассчитанное в каждый момент времени (т. Е. Каждый столбец), одинаково. Это согласуется независимо от количества генераторов.

Я не уверен, является ли это математической проблемой или проблемой программирования (возможно, обе), и мне кажется, что я вызываю ode45 неправильно, но я не уверен и не смог выяснить, что не так в течение нескольких дней , Любая помощь будет оценена.

Ответы [ 2 ]

0 голосов
/ 04 сентября 2018

Вы действительно использовали линию 

[t,y] = ode45(@(t,y)kuramoto(theta(:,1),K,N,omega),tspan,t0);

в вашем рабочем коде? Тогда результат здесь определенно неверный. Используйте 

[t,y] = ode45(@(t,u)kuramoto(u,K,N,omega),tspan,t0);

для получения результатов, которые хотя бы связаны с интеграцией RK4. То есть использовать объявленную локальную переменную / параметр анонимной функции при вычислении ее значения. (Использование u вместо y или theta, чтобы не использовать повторно имя переменной, которая также используется в более глобальной области. Вместо нее можно использовать thetalocal, если требуются самодокументируемые имена переменных.)

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

0 голосов
/ 28 августа 2018

Вы должны использовать вывод ode45. Рунге Кутта в том виде, в каком он реализован, в конечном итоге будет нестабильным, если вы выберете слишком большой размер шага. Весь смысл ode45 в том, что он внутренне работает по схемам Runge Kutta 4 и Runge Kutta 5. Если результаты одного шага интеграции отличаются, то ode45 будет сокращать временной шаг до тех пор, пока результаты не будут сопоставимы. Использование необработанного метода, как вы делаете, очевидно, не сделает этого.

Технически, такие вещи, как ode45, называются «встроенными методами Рунге Кутты». Вот один из таких методов: https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta%E2%80%93Fehlberg_method Они эффективны, потому что методы Рунге Кутты различного порядка многократно используют одни и те же оценки функций.

Все это, как говорится, вы должны обнаружить, что, если вы достаточно сократите свой временной шаг, результаты будут практически идентичны. Единственная причина, по которой они отличаются, заключается в том, что ode45 внутренне уточняет шаг по времени, когда обнаруживает, что решение может быть неточным.

...