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

Question

I am trying to solve two equations with complex coefficients using ode45. But iam getting an error message as "Inputs must be floats, namely single or 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])

How can i fix this ? Can somebody can help me ?

Hmmm... Can you separate your diff Eq. into real and imaginary parts and use ode45 on each separately? — anon01
@anon0909 I do not think that the problem is in complex numbers, ode45 should easily handle them. — madbitloman

TroyHaskin TroyHaskin · Accepted Answer · 2015-07-31T21:02:42

Per the MathWorks Support Team, the "ODE solvers in MATLAB 5 (R12) and later releases properly handle complex valued systems." So the complex numbers are the not the issue.

The error "Inputs must be floats, namely single or double." stems from your definition of f using Symbolic Variables that are, unlike complex numbers, not floats. The easiest way to get around this is to not use the Symbolic Toolbox at all; just makes Eqns an anonymous function:

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]);

That being said, I'd like to point out that numerically time integrating this system over 300 microseconds (I assume without units given) will take a long time since your coefficient matrix has imaginary eigenvalues on the order of 10E+10. The extremely short wavelength of those oscillations will more than likely be resolved by Matlab's adaptive methods, and that will take a while to solve for a time span just a few orders greater than the wavelength.

I'd, therefore, suggest an analytical approach to this problem; unless it is a stepping stone another problem that is non-analytically solvable.

Systems of ordinary differential equations of the form

,

which is a linear, homogenous system with a constant coefficient matrix, has the general solution

,

where the m-subscripted exponential function is the matrix exponential. Therefore, the analytical solution to the system can be calculated exactly assuming the matrix exponential can be calculated. In Matlab, the matrix exponential is calculate via the expm function. The following code computes the analytical solution and compares it to the numerical one for a short time span:

%   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');

The comparison plot is:

The output of ode45 matches the magnitude and real parts of the solution very well, but the imaginary portion is out-of-phase by exactly π. However, since ode45's error estimator only looks at norms, the phase difference is not noticed which may lead to problems depending on the application.

It will be noted that while the matrix exponential solution is far more costly than ode45 for the same number of time vector elements, the analytical solution will produce the exact solution for any time vector of any density given to it. So for long time solutions, the matrix exponential can be viewed as an improvement in some sense.

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

1 Answers