ODE15s Errors defining variable and ODE functions

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);

Answer 1

[Walter Roberson]

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)

Q = function_handle with value:
    @(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)])

eqs = 

vars = 

[M,F] = massMatrixForm(eqs,vars)

M = 

F = 

f = M\F

f = 

ode = odeFunction(f, vars)

ode = function_handle with value:
    @(t,in2)[((in2(1,:).*-7.508314e+7+in2(2,:).*in2(1,:).*4.157e+3+in2(1,:).^2.*2.5e+5).*(-3.215314705597989e-6))./in2(2,:);((in2(2,:).*1.5e+5+in2(1,:).*5.8e+2-in2(2,:).*in2(1,:).*5.29e+2).*1.607657352798994e-3)./in2(1,:)]

[t,TP] = ode15s(ode, tspan, [T0+1, P0]);
plot(t, TP(:,1))
title('T')

figure
plot(t, TP(:,2))
title('P')

simplify(f)

ans = 

subs(f,vars,[T0;P0])

ans = 

vpa(ans)

ans = 

[Michela Benazzi]

Thank you so much for helping me plot the graphs, you are a lifesaver!

I tried extending the range of the graphs (I need to be able to predict up to 5 min) and the graphs disappear. I was wondering if you get the same when running the following:

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])
vpa(ans)

[Walter Roberson]

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