Как решить уравнения с комплексными коэффициентами, используя ode45 в MATLAB?

Согласно группе поддержки MathWorks, "решатели ODE в MATLAB 5 (R12) и более поздние выпуски должным образом обрабатывают комплексные системы». Так что комплексные числа не проблема.

Ошибка «Входные данные должны быть с плавающей запятой, а именно одинарные или двойные». вытекает из вашего определения f с использованием символических переменных, которые, в отличие от комплексных чисел, не являются числами с плавающей запятой. Самый простой способ обойти это — вообще не использовать Symbolic Toolbox; просто делает Eqns анонимной функцией:

Eqns= @(t,X) [-(X(1)*23788605396486326904946699391889*1i)/38685626227668133590597632 + (X(2)*23788605396486326904946699391889*1i)/38685626227668133590597632; (X(2)*23788605396486326904946699391889*1i)/38685626227668133590597632 + X(1)*(- 2500000 + (5223289665997855453060886952725538686654593059791*1i)/324518553658426726783156020576256)] ;

[t,Xabc]=ode45(Eqns,[0 300*10^-6],[0 1]);

При этом я хотел бы отметить, что численное время интегрирования этой системы в течение 300 микросекунд (я предполагаю, что без указания единиц измерения) займет много времени, поскольку ваша матрица коэффициентов имеет мнимые собственные значения порядка 10E+10. Чрезвычайно короткая длина волны этих колебаний, скорее всего, будет разрешена с помощью адаптивных методов Matlab, и это займет некоторое время, чтобы решить временной интервал, всего на несколько порядков превышающий длину волны.

Поэтому я бы предложил аналитический подход к этой проблеме; если это не ступенька к другой проблеме, которая не может быть решена аналитически.

Системы обыкновенных дифференциальных уравнений вида

,

которая представляет собой линейную однородную систему с постоянной матрицей коэффициентов, имеет общее решение

,

где m-подстрочная экспоненциальная функция является матричной экспоненциальной. Следовательно, аналитическое решение системы может быть вычислено точно, если можно вычислить матричную экспоненту. В Matlab экспоненциальная матрица вычисляется с помощью функции expm. Следующий код вычисляет аналитическое решение и сравнивает его с числовым за короткий промежуток времени:

%   Set-up
A    = [-23788605396486326904946699391889i/38685626227668133590597632,23788605396486326904946699391889i/38685626227668133590597632;...
        -2500000+5223289665997855453060886952725538686654593059791i/324518553658426726783156020576256,23788605396486326904946699391889i/38685626227668133590597632];
Eqns = @(t,X) A*X;
X0   = [0;1];

%   Numerical
options = odeset('RelTol',1E-8,'AbsTol',1E-8);
[t,Xabc]=ode45(Eqns,[0 1E-9],X0,options);

%   Analytical
Xana = cell2mat(arrayfun(@(tk) expm(A*tk)*X0,t,'UniformOutput',false)')';

k = 1;
%   Plots
figure(1);
    subplot(3,1,1)
        plot(t,abs(Xana(:,k)),t,abs(Xabc(:,k)),'--');
        title('Magnitude');
    subplot(3,1,2)
        plot(t,real(Xana(:,k)),t,real(Xabc(:,k)),'--');
        title('Real');
        ylabel('Values');
    subplot(3,1,3)
        plot(t,imag(Xana(:,k)),t,imag(Xabc(:,k)),'--');
        title('Imaginary');
        xlabel('Time');

График сравнения:

Выход ode45 очень хорошо соответствует величине и действительной части решения, но мнимая часть не совпадает по фазе ровно на . Однако, поскольку оценщик ошибок ode45 рассматривает только нормы, разница фаз не замечается, что может привести к проблемам в зависимости от приложения.

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