The graph is not showing up even though the ODE solver is working...

Question

StillANovice 2018-8-1

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/413163-the-graph-is-not-showing-up-even-though-the-ode-solver-is-working

重新打开： Walter Roberson 2018-12-22

Function:

function dAdT = polyadiabatic(T,A)
dAdT = zeros(2,1);
% Kinetics
Rp = (kp).*rho.*A(1).*((2.*ki.*((rho.*A(1).^3)))./ktc).^0.5;
% Formulae
ki = 2.019.*10.*exp(-13810./A(2)); % m6/kg.s
kp = 1.009.*(10.^5).*exp(-3557./A(2)); % m3/kg.s
%kp = 100900;
ktc = 2.205.*(10.^7).*exp(-844./A(2)).*exp(-2.*((A1.*wp)+(A2.*(wp.^2))+(A3.*(wp.^3)))); % m3/kg.s
A1 = 2.57-((5.05.*0.001).*A(2));
A2 = 9.56-((1.76.*0.01).*A(2));
A3 = -3.032+((7.85.*0.001).*A(2));
rho = 845-(A(2)-353)+((200+(A(2)-353)).*wp); % kg/m3
wp = 1-A(1);
% Pressure
Pv_s = exp(135.23395+((-9226.418).*(A(2).^(-1)))+((-17.99862).*log(A(2)))+(0.01669633.*A(2))); % Vapour pressure, Pa
% Parameters
deltaH = -6.7.*(10.^5); % J/kg
cp = 1.884.*1000; % J/kg.K
% Mass balance
dAdT(1) = -(Rp./rho);
% Energy balance
dAdT(2) = -(Rp.*deltaH)./(cp.*rho);
end

Calling the function:

function [T,A] = call_polyadiabatic()
tspan = [0 452100];
% Initial conditions
A1_0 = 1;
A2_0 = 318;
[T,A] = ode15s(@polyadiabatic,tspan,[A1_0 A2_0]);
yyaxis left
ylabel("Temperature, K")
hold on
grid on
plot(T,A(:,2),'linewidth',1)
yyaxis right
ylabel("Vapour pressure, Pa")
plot(T,Pv_s,'linewidth',1)
hold off
xlabel("Time, s")
title("Adiabatic Temperature and Pressure Profiles")
end

Error:

Undefined function or variable 'Pv_s'.
Error in call_polyadiabatic (line 31)
plot(T,Pv_s,'linewidth',1)

How should I modify my code in order for Pv_s to be plotted?

5 个评论
显示 3更早的评论隐藏 3更早的评论

StillANovice 2018-8-30

Appreciate your consideration :)

Star Strider 2018-8-30

As always, my pleasure.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2018-8-1

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/413163-the-graph-is-not-showing-up-even-though-the-ode-solver-is-working#answer_331247

在 MATLAB Online 中打开

First, I had to rearrange ‘polyadiabatic’ to get it to work, because in the version you posted, several statements were out of order.

This version works:

function dAdT = polyadiabatic(T,A)
dAdT = zeros(2,1);
% Formulae
ki = 2.019.*10.*exp(-13810./A(2)); % m6/kg.s
kp = 1.009.*(10.^5).*exp(-3557./A(2)); % m3/kg.s
%kp = 100900;
A1 = 2.57-((5.05.*0.001).*A(2));
A2 = 9.56-((1.76.*0.01).*A(2));
A3 = -3.032+((7.85.*0.001).*A(2));
wp = 1-A(1);
rho = 845-(A(2)-353)+((200+(A(2)-353)).*wp); % kg/m3
ktc = 2.205.*(10.^7).*exp(-844./A(2)).*exp(-2.*((A1.*wp)+(A2.*(wp.^2))+(A3.*(wp.^3)))); % m3/kg.s
% Kinetics
Rp = (kp).*rho.*A(1).*((2.*ki.*((rho.*A(1).^3)))./ktc).^0.5;
% Pressure
Pv_s = exp(135.23395+((-9226.418).*(A(2).^(-1)))+((-17.99862).*log(A(2)))+(0.01669633.*A(2))); % Vapour pressure, Pa
% Parameters
deltaH = -6.7.*(10.^5); % J/kg
cp = 1.884.*1000; % J/kg.K
% Mass balance
dAdT(1) = -(Rp./rho);
% Energy balance
dAdT(2) = -(Rp.*deltaH)./(cp.*rho);
end

Second, you need to calculate ‘Pv_s’ in ‘call_polyadiabatic’ in order to plot it. Add it (with a minor edit to use the second column of ‘A’) just after your ode15s call:

[T,A] = ode15s(@polyadiabatic,tspan,[A1_0 A2_0]);
Pv_s = exp(135.23395+((-9226.418).*(A(:,2).^(-1)))+((-17.99862).*log(A(:,2)))+(0.01669633.*A(:,2))); % Vapour pressure, Pa

Then the plot works.

Also, it might be more efficient to replace:

exp(135.23395+((-9226.418).*(A(2).^(-1)))

with:

exp(135.23395+((-9226.418)./A(2))

I will let you experiment with that.