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

Question

Я пытаюсь решить два уравнения с комплексными коэффициентами, используя ode45. Но iam получает сообщение об ошибке «Входы должны быть поплавками, а именно одиночными или double».

X = sym(['[',sprintf('X(%d) ',1:2),']']); 

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

f=@(t,X)[Eqns]; 

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

Как я могу это исправить? Может кто-нибудь может мне помочь?

Answer 1

В «MathWorks Support Team» «Решения 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 -подписью экспоненциальная функция - matrix exponential. Следовательно, аналитическое решение системы может быть рассчитано точно, если предположить, что матричная экспонента может быть вычислена. В 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для того же число временных векторных элементов, аналитическое решение будет производить точное решение для любого вектора времени любой плотности данной к нему. Поэтому для решений в течение длительного времени матричную экспоненту можно рассматривать как улучшение в некотором смысле.