Solve a system of linear ODEs, and plot them

I have a system of Linear ODEs. and I tried the following code, based on an example from matlab.com

clear all
clc
format short e
syms x(t) y(t) v(t) q(t) u(t) s(t) z(t)
A=[-108.15802678, -96035700000, 99.9039643,  100,           -9603.57,  -96.0357,   100;
   -0.00104128,   -96035700000, 0,           0,              0,         0,         0;
    102.8581245,   0,          -100.0960357, 100,            9603.57,   0,         0;
    0.0988755,     0,           0.0960357,  -100,            0,         0,         0;
   -102.957,       0,           100,         0,             -9603.57,   0,         0;
   -5.10111,       0,           0,           0,              0,        -96.0357,   100;
    5.10111,       0,           0,           0,              0,         96.0357,  -100];

B=[ 0.001;         0.001;       0;           0;              0;         0;         0];
tspan = [0, 0.5, 10];
R=[x; y; v; q; u; s; z];
odes= diff(R) == A*R + B

C = R(0) == [10^-7, 10^-7, 0, 1, 0, 0.1, 0];
[xSol(t), ySol(t), vSol(t), qSol(t), uSol(t), sSol(t), zSol(t)] = dsolve(odes,C);

xSol(t) = simplify(xSol(t))
ySol(t) = simplify(ySol(t))
vSol(t) = simplify(vSol(t))
qSol(t) = simplify(qSol(t))
uSol(t) = simplify(uSol(t))
sSol(t) = simplify(sSol(t))
zSol(t) = simplify(zSol(t))

clf
fplot(xSol)
hold on
fplot(ySol)
hold on
fplot(vSol)
hold on
fplot(qSol)
hold on
fplot(uSol)
hold on
fplot(sSol)
hold on
fplot(zSol)
grid on
legend('xSol','ySol','vSol','qSol','uSol','sSol','zSol','Location','best')

Although I get the results, I do not get any plots. It just shows a plot with no graph/lines/points on it. Any opinion on how to solve it?

I really appreciate any help,

[M]

I got the idea for this code from here:

https://www.mathworks.com/help/symbolic/solve-a-system-of-differential-equations.html#buxuujb

Answer 1

[Walter Roberson]

You built R as a column vector but your vector of constants 10^-7 and so on is a row vector. You will accidentally creating conflicting constraints.

syms x(t) y(t) v(t) q(t) u(t) s(t) z(t)
A=[-108.15802678, -96035700000, 99.9039643,  100,           -9603.57,  -96.0357,   100;
   -0.00104128,   -96035700000, 0,           0,              0,         0,         0;
    102.8581245,   0,          -100.0960357, 100,            9603.57,   0,         0;
    0.0988755,     0,           0.0960357,  -100,            0,         0,         0;
   -102.957,       0,           100,         0,             -9603.57,   0,         0;
   -5.10111,       0,           0,           0,              0,        -96.0357,   100;
    5.10111,       0,           0,           0,              0,         96.0357,  -100];

B=[ 0.001;         0.001;       0;           0;              0;         0;         0];
tspan = [0, 0.5, 10];
R=[x; y; v; q; u; s; z];
odes= diff(R) == A*R + B

C = R(0) == [10^-7, 10^-7, 0, 1, 0, 0.1, 0].';   %IMPORTANT CHANGE
sol = dsolve(odes,C);

xSol(t) = sol.x;
ySol(t) = sol.y;
vSol(t) = sol.v;
qSol(t) = sol.q;
uSol(t) = sol.u;
sSol(t) = sol.s;
zSol(t) = sol.z;

fplot([xSol, ySol, vSol, qSol, uSol, sSol, zSol]);

legend({'x(t)','y(t)','v(t)','q(t)','u(t)','s(t)','z(t)'},'Location','best')

However...

You need to restrict the t range to about [0 1e-3] or else you will get NaN due to floating point overflows.

All of the solutions to the equations generate complex numbers. They all include roots of a polynomial of degree 6, and the imaginary components are not canceling out: the imaginary components are in the thousands. And fplot() will not plot imaginary values, so you need to take real() or imag()

In my testing, the solutions appear to be coming out constant in time, at least to the number of digits I processed them to. xSol(1e-10) is the same as xSol(1e-5) is the same as xSol(1e-3) to within the number of digits used.

[M]

Thank you very much for your detailed response. I will consider these and hopefully will find a way through it.

Thanks a lot!

[Walter Roberson]

I seem to be getting contradictions as to whether the results are complex-valued or not.

If you take something like,

Zimg = @(EXPR) mapSymType(EXPR, 'numeric', @(v) piecewise(abs(imag(v))<1e-10, real(v), v))

xSim = simplify(Zimg(vpa(xSol)))

then imaginary parts can disappear.

The purpose of Zimg is to remove small imaginary components that are due to round-off error. This is not the same thing as just taking real() of all of the constants in the expression, because it turns out that the expression has two very significant imaginary components. But it also turns out in my test that those two imaginary components cancel out, leaving only real values.

But when I probe, whether I get imaginary components in practice seems to depend on exactly how I phrase the question, for what should be mathematically identical formulations.

I am, by the way, now seeing some variations in the values of the expression with time, when t is sufficiently small -- below about 1/1000. It seems to reach an asymptope after that.

[M]

Walter,

I really appreciate your time and help. Also, I forgot to say that these euqations may (and in fact they will) reach an equlibrium at some point (probably in a very short period of time). So, getting constant results for x values (or the rest) is not something that worries me at this point. All in all, I really appreciate your detailed response. I need a faster machine to get results as fast as you do (LOL).

Thanks again