Plot coming up blank
1 view (last 30 days)
Show older comments
Hi everyone,
I am trying to plot
function ODE15s20Aug2018
%%VARIABLES
% y(1)
% y(2)
% y(3)
% y(4)
% y(5)
% y(6)
% y(7)
% y(8)
%%EQUATIONS
%(y(1)' = ((1/A) * (B * C– B * y(1))) – (((y(1) * E * F) – (G * ((y(3) * y(8).^2)/(y(8).^2 + H * y(8) + H * I))))/ ((1 / J) + (G / ((1+ ((K * y(5)) / (LL * y(3)))) * M))))
%(y(2)' = ((1/A) * (B * N– B * y(2))) – ((O * (1 + ((K * y(5))/ (P * y(4)))))* (((y(2) * E * F) / (S)) – (((y(4) * y(8).^2) / (y(8).^2 + Q * y(8) + Q * R)))))
%(y(3)) /dt = ( ((y(1) * E * F) – (G * ((y(3) * y(8).^2)/(y(8).^2 + H * y(8) + H * I))))/ ((1 / J) + (G / ((1+ ((K * y(5)) / (LL * y(3)))) * M)))) – (0.162 * (exp(-5153/F)) * (((y(5) * ((y(3) * H * I) / (y(8).^2 + H * y(8) + H * I))) / T) - 1).^3 * ( (U) / (( y(5) * ((y(3) * H * I) / (y(8).^2 + H * y(8) + H * I))) /(T))))
%(y(4)) /dt = ((O* (1 + ((K * y(5))/ (P * y(4)))))* (((y(2) * E * F) / (S)) – (((y(4) * y(8).^2) / (y(8).^2 + Q * y(8) + Q * R))))) + (y(6) * (V* (1 – ((y(5) * ((y(4) * Q * R) / (y(8).^2 + Q * y(8) + Q * R))) / W))))
%(y(5)' = (y(6) * (V* (1 – ((y(5) * ((y(4) * Q * R) / (y(8).^2 + Q * y(8) + Q * R))) / W)))) – (0.162 * (exp(-5153/F)) * (((y(5) * ((y(3) * H * I) / (y(8).^2 + H * y(8) + H * I))) / T) - 1).^3 * ( (U) / (( y(5) * ((y(3) * H * I) / (y(8).^2 + H * y(8) + H * I))) /(T))))
%(y(6)' = -y(6) * (V* (1 – ((y(5) * ((y(4) * Q * R) / (y(8).^2 + Q * y(8) + Q * R))) / W))) * (X / Y)
%(y(7)' = (0.162 * (exp(-5153/F)) * (((y(5) * ((y(3) * H * I) / (y(8).^2 + H * y(8) + H * I))) / T) - 1).^3 * ( (U) / (( y(5) * ((y(3) * H * I) / (y(8).^2 + H * y(8) + H * I))) /(T)))) * (Z/AA)
%y(8) + 2 * y(5) – ((y(3) * H * y(8))/(y(8).^2 + H * y(8) + H * I)) – 2 * (((y(3) * H * I) / (y(8).^2 + H * y(8) + H * I))) – (((y(4) * Q * y(8))/ (y(8).^2 + Q * y(8) + Q * R))) – 2 * (((y(4) * Q * R) / (y(8).^2 + Q * y(8) + Q * R))) – (AB/y(8)) = 0
%%SINGULAR MASS MATRIX
M = diag([1 1 1 1 1 1 1 0 ]);
%%INITIAL VALUES
y0 = zeros(8,1);
y0(1)= 1e-9;
y0(2)= 1e-9;
y0(3)= 1e-9;
y0(4)= 1e-9;
y0(5)= 1e-9;
y0(6)= 4.99e-9;
y0(7)= 1e-9;
y0(8)= 7.413e-6;
tspan = linspace(0,183000,30);
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',1e-4,'Vectorized','off');
[t,y] = ode15s(@f,tspan,y0,options);
plot(t,y)
% --------------------------------------------------------------------------
function out = f(t,y)
% PARAMETERS
A = 1.5e-6;
B = 1.67e-5;
C = 6.51e-2;
E = 8.314;
F = 323.15;
G = 149;
H = 6.24;
I = 5.68e-5;
J = 4.14e-6;
K = 1.39e-9;
LL = 2.89e-9;
M = 8.4e-4;
N = 0;
O = 9.598e-4;
P = 3.53e-9;
Q = 1.7e-3;
R = 6.55e-8;
S = 5.15e3;
T = 1.07e-7;
U = 10;
V = 8.42e-4;
W = 3.2e-1;
X = 100.086;
Y = 2703;
Z = 258.3;
AA = 2540;
AB = 5.3e-8;
out =[ ((1/A) * (B * C - B * y(1,:))) - (((y(1,:) * E * F) - (G * ((y(3,:) * y(8,:).^2)/(y(8,:).^2 + H * y(8,:) + H * I))))/ ((1 / J) + (G / ((1+ ((K * y(5,:)) / (LL * y(3,:)))) * M))))
((1/A) * (B * N - B * y(2,:))) - ((O * (1 + ((K * y(5,:))/ (P * y(4,:)))))* (((y(2,:) * E * F) / (S)) - (((y(4,:) * y(8,:).^2) / (y(8,:).^2 + Q * y(8,:) + Q * R)))))
( ((y(1,:) * E * F) - (G * ((y(3,:) * y(8,:).^2)/(y(8,:).^2 + H * y(8,:) + H * I))))/ ((1 / J) + (G / ((1+ ((K * y(5,:)) / (LL * y(3,:)))) * M)))) - (0.162 * (exp(-5153/F)) * (((y(5,:) * ((y(3,:) * H * I) / (y(8,:).^2 + H * y(8,:) + H * I))) / T) - 1).^3 * ( (U) / (( y(5,:) * ((y(3,:) * H * I) / (y(8,:).^2 + H * y(8,:) + H * I))) /(T))))
((O* (1 + ((K * y(5,:))/ (P * y(4,:)))))* (((y(2,:) * E * F) / (S)) - (((y(4,:) * y(8,:).^2) / (y(8,:).^2 + Q * y(8,:) + Q * R))))) + (y(6,:) * (V* (1 - ((y(5,:) * ((y(4,:) * Q * R) / (y(8,:).^2 + Q * y(8,:) + Q * R))) / W))))
(y(6,:) * (V* (1 - ((y(5,:) * ((y(4,:) * Q * R) / (y(8,:).^2 + Q * y(8,:) + Q * R))) / W)))) - (0.162 * (exp(-5153/F)) * (((y(5,:) * ((y(3,:) * H * I) / (y(8,:).^2 + H * y(8,:) + H * I))) / T) - 1).^3 * ( (U) / (( y(5,:) * ((y(3,:) * H * I) / (y(8,:).^2 + H * y(8,:) + H * I))) /(T))))
-y(6,:) * (V* (1 - ((y(5,:) * ((y(4,:) * Q * R) / (y(8,:).^2 + Q * y(8,:) + Q * R))) / W))) * (X / Y)
(0.162 * (exp(-5153/F)) * (((y(5,:) * ((y(3,:) * H * I) / (y(8,:).^2 + H * y(8,:) + H * I))) / T) - 1).^3 * ( (U) / (( y(5,:) * ((y(3,:) * H * I) / (y(8,:).^2 + H * y(8,:) + H * I))) /(T)))) * (Z/AA)
y(8,:) + 2 * y(5,:) - ((y(3,:) * H * y(8,:))/(y(8,:).^2 + H * y(8,:) + H * I)) - 2 * (((y(3,:) * H * I) / (y(8,:).^2 + H * y(8,:) + H * I))) - (((y(4,:) * Q * y(8,:))/ (y(8,:).^2 + Q * y(8,:) + Q * R))) - 2 * (((y(4,:) * Q * R) / (y(8,:).^2 + Q * y(8,:) + Q * R))) - (AB/y(8,:))];
However, my plot comes out blank. What can I do?
See attached for details.
0 Comments
Answers (1)
Walter Roberson
on 20 Aug 2018
>> ODE15s20Aug2018
Warning: Failure at t=5.450176e-06. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.355253e-20) at time t.
> In ode15s (line 668)
In ODE15s20Aug2018 (line 48)
Your function contains a singularity that is encountered in the very first time step after time 0, so your t is coming out as a single value, and your y is coming out as a single row vector. When you ask to plot() a single point, then because by default it only plots lines, you do not see anything in your plot. If you were to add a marker style in your plot, you would have gotten a plot with a single point drawn.
If you have the symbolic toolbox, I recommend that you look at the documentation for odeFunction() and follow the first example there to see the workflow of converting a symbolic system of equations for use with the numeric ode* routines.
2 Comments
See Also
Categories
Find more on Linear Algebra in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!