stiff ODE system plots disappear when I increase axes limits

2 views (last 30 days)
Michela Benazzi
Michela Benazzi on 13 Oct 2021
Answered: Abhiroop Rastogi on 25 Oct 2021
Hello everyone! I recently asked for help to solve my ODE equations and I thought I might have to open another thread for this one.
My plots disappear when I set limits for them. I do see a line (straight line, probably due to the fact that it is only a tiny portion of the graph) before doing so, and nothing after adding limits on lines 24 and 28.
I need to be able to predict the rate of change until 5 min (300 second, with unit set on tspan as 0.001 s).
Thank you!
clc;
clear;
close all;
Q = @(v) sym(v);
%Units
cp = Q(29) ; %J/g*K
r = Q(8.314) ; %J/g*K
cv = cp - r;
tspan = [0 0.001] ; %sec
T0 = 300 ; %K
P0 = 20000000 ; %Pa (20 MPa in Pa due to R being in terms of m3)
syms T(t) P(t)
ode1 = diff(T) == 2.0*(T/P)*(r/cv)*(T0-T)-4*10^(-3)*(T/P)*(r^(2)/cv)*(P-P0);
ode2 = diff(P) == 2.0*(P/T)*(r/cv)*(T0-T)-4*10^(-3)*r*(cp/cv)*(P-P0);
odes = [ode1; ode2];
[eqs,vars] = reduceDifferentialOrder(odes, [T(t), P(t)]);
[M,F] = massMatrixForm(eqs,vars) ;
f = M\F ;
ode = odeFunction(f, vars) ;
[t,TP] = ode15s(ode, tspan, [T0+1, P0]);
plot(t, TP(:,1))
title('T')
axis([0 300000 120 300])
figure
plot(t, TP(:,2))
title('P')
axis([0 300000 0 20000000])
simplify(f)
subs(f,vars,[T0;P0])

Answers (1)

Abhiroop Rastogi
Abhiroop Rastogi on 25 Oct 2021
Hi Michela,
The limits that you have chosen are very much out of bounds what the data is achieving. You can either let the limits be chosen by default or you can set the limits as shown below.
clc;
clear;
close all;
Q = @(v) sym(v);
%Units
cp = Q(29) ; %J/g*K
r = Q(8.314) ; %J/g*K
cv = cp - r;
tspan = [0 0.001] ; %sec
T0 = 300 ; %K
P0 = 20000000 ; %Pa (20 MPa in Pa due to R being in terms of m3)
syms T(t) P(t)
ode1 = diff(T) == 2.0*(T/P)*(r/cv)*(T0-T)-4*10^(-3)*(T/P)*(r^(2)/cv)*(P-P0);
ode2 = diff(P) == 2.0*(P/T)*(r/cv)*(T0-T)-4*10^(-3)*r*(cp/cv)*(P-P0);
odes = [ode1; ode2];
[eqs,vars] = reduceDifferentialOrder(odes, [T(t), P(t)]);
[M,F] = massMatrixForm(eqs,vars) ;
f = M\F ;
ode = odeFunction(f, vars) ;
[t,TP] = ode15s(ode, tspan, [T0+1, P0]);
plot(t, TP(:,1))
title('T')
%axis([0 300000 120 300])
axis([0 1e-3 300.999999992 301]) % modified limits
figure
plot(t, TP(:,2))
title('P')
%axis([0 300000 0 20000000])
axis([0 1e-3 1.999994e7 2e7]) % modified limits
simplify(f)
ans = 
subs(f,vars,[T0;P0])
ans = 

Community Treasure Hunt

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

Start Hunting!

Translated by