MATLAB Answers

Solve a system of linear ODEs, and plot them

5 views (last 30 days)
M
M on 20 Sep 2021
Commented: M on 20 Sep 2021
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,

Accepted Answer

Walter Roberson
Walter Roberson on 20 Sep 2021
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.
  3 Comments
M
M on 20 Sep 2021
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

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!