Plot coming up blank
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Réponses (1)
Walter Roberson
le 20 Août 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.
Voir également
Catégories
En savoir plus sur Linear Algebra dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!