如何在 MATLAB 中使用 ode45 求解具有复数系数的方程？

Question

我正在尝试使用 ode45 求解两个具有复系数的方程。但是我收到一条错误消息“输入必须是浮点数，即单个或双倍。”

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 solvers in MATLAB 5 (R12) and later releases properly handle complex valued systems." 所以复数不是问题。

错误 "Inputs must be floats, namely single or double." 源于您使用符号变量定义 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 微秒以上（我假设没有给定单位）对这个系统进行数值时间积分将花费很长时间，因为您的系数矩阵具有 [= 数量级的虚特征值14=]。这些振荡的极短波长很可能会被 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 ，解析解将产生给定任何密度的任何时间向量的精确解。所以对于长期解决方案，矩阵指数在某种意义上可以看作是一种改进。

如何在 MATLAB 中使用 ode45 求解具有复数系数的方程？

How to solve equations with complex coefficients using ode45 in MATLAB?

matlab

ode

ode45