1

我正在尝试使用 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])

我怎样才能解决这个问题 ？有人可以帮助我吗？

4

1 回答 1

1

根据 MathWorks 支持团队，“MATLAB 5 (R12) 及更高版本中的 ODE 求解器可以正确处理复值系统。” 所以复数不是问题。

错误“输入必须是浮点数，即单精度或双精度。” 源于您对f使用符号变量的定义，这些变量与复数不同，不是浮点数。解决这个问题的最简单方法是根本不使用符号工具箱；只是做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 的自适应方法来解决，而这将需要一段时间才能解决仅比波长大几个数量级的时间跨度。

因此，我建议对这个问题采取一种分析方法；除非它是一个垫脚石，否则另一个无法通过分析解决的问题。

常微分方程组的形式

常微分方程的线性、齐次、常系数系统的方程,

它是一个具有常系数矩阵的线性齐次系统，具有通解

ODE 系统的一般矩阵指数解,

其中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 相同数量的时间向量元素的成本要高得多，但解析解将为给定的任何密度的任何时间向量产生精确解。因此，对于长时间的解决方案，矩阵指数在某种意义上可以被视为一种改进。

于 2015-07-31T21:02:42.307 回答