error in ode45 - must return a column vector

Question

Gustavo 2012-4-25

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/36550-error-in-ode45-must-return-a-column-vector

采纳的回答： Jan

在 MATLAB Online 中打开

Hi everyone,

Im trying to solve a PDE using ode45. But the following error are ocurring: 'EDP must return a column vector'.

I put my code blow for you guys help me to find the error.

------------------------------------------------------------------

global pipe Tin I Tg W T TW;
T = [];
Ts = [];
TW = [];
Tws = [];
% condições iniciais
Iini = IC151(1,2);
Tgini = TA075(1,2);
Wini = FA061(1,2);
Tini = TA072(1,2);
T = Tini;   
pipe.Tref = 250;  
pipe.rho = 903 - 0.672 * pipe.Tref;     
pipe.Cp = 1820 + 3.76 * pipe.Tref;   
pipe.rhow = 7800;      
pipe.Cw = 450;  
pipe.ri = 0.009; 
pipe.ro = 0.013; 
pipe.di = pipe.ri * 2;   
pipe.do = pipe.ro * 2;  
pipe.Ai = pi * pipe.ri^2; 
pipe.Ao = pi * pipe.ro^2; 
pipe.ho = 11;
pipe.hi=800;   
pipe.L=172;   
pipe.S=150; 
                                                       pipe.dX=pipe.L /pipe.S;    
pipe.n0=0.43;  
pipe.G=1.82;  
                                                     pipe.vs=Wini /pipe.Ai;   
pipe.gamma = (pipe.n0 * pipe.G) / (pipe.Ao * pipe.rhow * pipe.Cw);
pipe.Tau1 = (pipe.Ai * pipe.rho * pipe.Cp) / (pi * pipe.di * pipe.hi);
pipe.Tau12 = (pipe.Ao * pipe.rhow * pipe.Cw) / (pi * pipe.di * pipe.hi);
pipe.Tau2 = (pipe.Ao * pipe.rhow * pipe.Cw) / (pi * pipe.do * pipe.ho);
pipe.c = (pipe.vs * pipe.Tau1 * (1 + (pipe.Tau2 / pipe.Tau12)));
TW = (Tini * pipe.Tau2) / (pipe.Tau12 + pipe.Tau2) + (Tgini * pipe.Tau12) / (pipe.Tau12 + pipe.Tau2) + (Iini * pipe.gamma * pipe.Tau12 * pipe.Tau2) / (pipe.Tau12 + pipe.Tau2);
for(i=1:pipe.S)
    T(i+1) = pipe.gamma * pipe.Tau2 * Iini + Tgini + (T(1) - pipe.gamma * pipe.Tau2 * Iini - Tgini) * exp(-(i * pipe.dX) / pipe.c);
    TW(i+1) = (T(i+1) * pipe.Tau2) / (pipe.Tau12 + pipe.Tau2) + (Tgini * pipe.Tau12) / (pipe.Tau12 + pipe.Tau2) + (Iini * pipe.gamma * pipe.Tau12 * pipe.Tau2) / (pipe.Tau12 + pipe.Tau2);
end
Ts = T(1,pipe.S+1);
Tws = TW(1,pipe.S+1);
n = 2;
for k = 1:2481
      I = IC151(n,2);              % from workspace
      Tin = TA072(n,2);            % from workspace
      Tg = TA075(n,2);             % from workspace
      W = FA061(n,2);              % from workspace
      [dT dTW] = ode45(@EDP,[0 1],[T,TW],[]);
      n= n+1; 
end

------------------------------------------------------------------

and the EDP file has:

-------------------------------------------------------------------

function [dT, dTW] = EDP(t,x)
global pipe W I Tg Tin T TW;
T = x(1:pipe.S+1);
T(1) = Tin;
TW = x(pipe.S+1+1:2*(pipe.S+1));
vs = W / pipe.Ai;
dT(1) = 0;
dTW(1) = 0;
for(i=2:pipe.S+1)
    dT(i) = -vs * ((T(i) - T(i-1)) / pipe.dX) + ((TW(i) - T(i)) / pipe.Tau1);
    dTW(i) = pipe.gamma * I + ((Tg - TW(i)) / pipe.Tau2) - ((TW(i) - T(i)) / pipe.Tau12);
end