ODE15s Errors defining variable and ODE functions
2 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Antworten (1)
Walter Roberson
am 13 Okt. 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 Kommentare
Walter Roberson
am 14 Okt. 2021
Bearbeitet: Walter Roberson
am 14 Okt. 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
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Produkte
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!