ODE15s Errors defining variable and ODE functions

조회 수: 6(최근 30일)
Michela Benazzi
Michela Benazzi 2021년 10월 12일
편집: Walter Roberson 2021년 10월 14일
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);

답변(1개)

Walter Roberson
Walter Roberson 2021년 10월 13일
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 = 
  댓글 수: 2
Walter Roberson
Walter Roberson 2021년 10월 14일
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

댓글을 달려면 로그인하십시오.

Community Treasure Hunt

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

Start Hunting!

Translated by