ODE15s Errors defining variable and ODE functions
2 views (last 30 days)
Show older comments
Hello everyone, I am trying to model changing T, P for a system with a system of stiff ODE's.
I get the following error:
Unrecognized function or variable 'T'.
Error in che561project1 (line 12)
ode1 = diff(T) == 2.0*(T/P)*(r/cv)*(T0-T)-4*10^(-3)*(T/P)*(r^(2)/cv)*(P-P0);
How do I define T?
I am also unsure of what is wrong with line 12. I am definitely going to run into more problems because it is my first time using MatLab, so I was wondering if the code looks good otherwise.
Thank you very much!! [Code below]
clc;
clear all;
close all;
%Units
cp = 29 ; %J/g*K
r = 8.314 ; %J/g*K
cv = cp - r;
tspan = [0 0.001] ; %min
T0 = 300 ; %K
P0 = 20 ; %MPa
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];
[t,T] = ode15s(@(t,T) 2.0*(T/P)*(r/cv)*(T0-T)-4*10^(-3)*(T/P)*(r^(2)/cv)*(P-P0), tspan, T0);
[t,P] = ode15s(@(t,P) 2.0*(P/T)*(r/cv)*(T0-T)-4*10^(-3)*r*(cp/cv)*(P-P0), tspan, P0);
0 Comments
Answers (1)
Walter Roberson
on 13 Oct 2021
Your equations involve components which subtract out T0 and P0. When you initialize the boundaries to T0 and P0 that results in no "force" being available to move the conditions away from T0 and P0. The plot I did here, I increased the initial temperature by 1 so that you could see something
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] ; %min
T0 = 300 ; %K
P0 = 20 ; %MPa
syms T(t) P(t)
ode1 = diff(T) == Q(2.0)*(T/P)*(r/cv)*(T0-T)-Q(4)*Q(10)^(-3)*(T/P)*(r^(2)/cv)*(P-P0);
ode2 = diff(P) == Q(2.0)*(P/T)*(r/cv)*(T0-T)-Q(4)*Q(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')
figure
plot(t, TP(:,2))
title('P')
simplify(f)
subs(f,vars,[T0;P0])
vpa(ans)
2 Comments
Walter Roberson
on 14 Oct 2021
Edited: Walter Roberson
on 14 Oct 2021
Your values are notably out of the range you expected. In particular, your x range is not even remotely close to 30000
Q = @(v) sym(v);
%Units
cp = Q(29) ; %J/g*K
r = Q(8.314) ; %J/g*K
cv = cp - r;
tspan = [0 300] ; %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])
xlim auto; ylim auto
figure
plot(t, TP(:,2))
title('P')
%axis([0 300000 0 20000000])
xlim auto; ylim auto
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!