MATLAB Answers

ODE Solver Function File Help - Must return a column vector

2 views (last 30 days)
Julia Carey
Julia Carey on 5 May 2021
Answered: Cris LaPierre on 5 May 2021
I keep trying to run an ode23tb but I keep getting a function file error saying my function file must return a column vector
function h=Julia_Prob1a_fun(t,y)
global T0 cA0 cB0 cC0 Pe ar PeT Le Nu mc mT N M dsi deta A B
k0=700;
Ea=111e3; %J/mol
vA=1;
vB=0;
vC=-0.36;
vD=-0.09;
nv=(N-2)*(M-2);
cA=y(1:nv);
cB=y(nv+1:2*nv);
cC=y(2*nv+1:3*nv);
T=y(3*nv+1:4*nv);
bA=zeros(nv,nv);
bB=bA;
bC=bA;
bT=bA;
for i=1:nv
lnK=5693.5/T(i)+1.077*log(T(i))+5.44e-4*T(i)-1.125e-7*T(i)^2-49170/T(i)^2-13.148;
K=exp(lnK);
RT=0.082*T(i)*1e-3;
r=k0*exp(-Ea/(8.314*T(i)))*(cA(i)*RT)^vA*(cB(i)*RT)^vB*(cC(i)*RT)^vC*(cC(i)*RT)^vD*(1-cC(i)*cC(i)/(cA(i)*cB(i)*K));
bA(i)=-mc*r;
bB(i)=-mc*r;
bC(i)=mc*r;
bT(i)=mT*r;
alfa=Pe/(2*dsi)+1/dsi^2;
gama=-Pe/(2*dsi)+1/dsi^2;
alfaT=Le*(PeT/(2*dsi)+1/dsi^2);
gamaT=Le*(-PeT/(2*dsi)+1/dsi^2);
end
for j=1:(M-2)
eps=ar^2*(1/(2*j*deta^2)+1/deta^2);
thta=eps*-1;
epsT=Le*eps;
thtaT=Le*thta;
for i=1:(N-2)
I=(N-2)*(j-1)+i;
if i==1
bA(I)=bA(I)+alfa(j)*cA0;
bB(I)=bB(I)+alfa(j)*cB0;
bC(I)=bC(I)+alfa(j)*cC0;
bT(I)=bT(I)+alfaT(j)*T0;
end
if i==(N-2)
bA(I)=bA(I)+gama(j)*cA(I);
bB(I)=bB(I)+gama(j)*cB(I);
bC(I)=bC(I)+alfa(j)*cC0;
bT(I)=bT(I)+alfaT(j)*T0;
end
if i==(N-2)
bA(I)=bA(I)+gama(j)*cA(I);
bB(I)=bB(I)+gama(j)*cB(I);
bC(I)=bC(I)+gama(j)*cC(I);
bT(I)=bT(I)+gamaT(j)*T(I);
end
if j==1
bA(I)=bA(I)+thta*cA(I);
bB(I)=bB(I)+thta*cB(I);
bC(I)=bC(I)+thta*cC(I);
bT(I)=bT(I)+thtaT*T(I);
end
if j==(M-2)
bA(I)=bA(I)+eps*cA(I);
bB(I)=bB(I)+eps*cB(I);
bC(I)=bC(I)+eps*cC(I);
bT(I)=bT(I)+epsT*(T(I)+Nu*deta*T0)/(1+Nu*deta);
end
end
end
hA=A*cA+bA;
hB=A*cB+bB;
hC=A*cC+bC;
hT=B*T+bT;
h=[hA;hB;hC;hT];
end
Main File
clc;
clear all;
global T0 cA0 cB0 cC0 tf Pe ar PeT Le Nu mc mT N M dsi deta A B
vm=0.066; %m/s
D=3.19e-7; %m2/s
P=1.013; % bar
rhob=21.6e3; %gm/m3
rhog=0.37e3; %gm/m3
kT=38.2/3600; %J/(m.s.K)
Cp=1; %J/(gm.K)
delH=-41.1e3; %J/mol;
UT=180/3600; %J/(m2.s.K)
L=70e-3; %m
R=6.35e-3; %m
ra=linspace(0,R); %various distances from the center of the pipe to R
vz=vm.*(1.-ra.^2./R.^2);
T0=773; %K
cT=P/(0.082*T0)*1e3 ; %mol/m3
cA0=0.27*cT; %CO;
cB0=0.72*cT; %H2O;
cC0=0; %CO2;
ar=L/R;
Pe=vz.*L/D;
PeT=rhog*Cp*vz.*L/kT;
Le=kT/(rhog*Cp*D);
Nu=UT*R/kT;
mc=L^2*rhob/D;
mT=Le*(-delH)*rhob*L^2/kT;
N=21; %z direction
M=21; %r direction
dsi=1/(N-1); %step in z direction
deta=1/(M-1); %step size in r direction
si=meshgrid(0:dsi:1); %dimensionless z
eta=meshgrid(0:deta:1)';%dimensionless eta
nv=(N-2)*(M-2);
alfa=Pe/(2*dsi)+1/dsi^2;
beta=-2/dsi^2-2*ar^2/deta^2;
gama=-Pe/(2*dsi)+1/dsi^2;
alfaT=Le*(PeT/(2*dsi)+1/dsi^2);
betaT=Le*(-2/dsi^2-2*ar^2/deta^2);
gamaT=Le*(-PeT/(2*dsi)+1/dsi^2);
A=zeros(nv,nv);
B=A;
for j=1:(M-2)
eps=ar^2*(1/(2*j*deta^2)+1/deta^2);
thta=ar^2*(-1/(2*j*deta^2)+1/deta^2);
epsT=Le*ar^2*(1/(2*j*deta^2)+1/deta^2);
thtaT=Le*ar^2*(-1/(2*j*deta^2)+1/deta^2);
for i=1:(N-2)
I=(N-2)*(j-1)+i;
A(I,I)=beta;
B(I,I)=betaT;
if i>1
A(I,I-1)=alfa(i);
B(I,I-1)=alfaT(i);
end
if i<(N-2)
A(I,I+1)=gama(i);
B(I,I+1)=gamaT(i);
end
if j>1
A(I,I-(N-2))=thta;
B(I,I-(N-2))=thtaT;
end
if j<(M-2)
A(I,I+(N-2))=eps;
B(I,I+(N-2))=epsT;
end
end
end
nv=(N-2)*(M-2);
cAi=cA0*ones(nv,1);
cBi=cB0*ones(nv,1);
cCi=1e-7*ones(nv,1);
Ti=T0*ones(nv,1);
yin=[cAi;cBi;cCi;Ti];
options = odeset('AbsTol',1e-8,'RelTol',1e-5,'MaxStep',0.1);
tf=1e-4;
tau=linspace(0,tf,100);
[t,y]=ode23tb(@Julia_Prob1a_fun,tau,yin',options);
cA=y(:,1:nv);
cB=y(:,nv+1:2*nv);
cC=y(:,2*nv+1:3*nv);
T=y(:,3*nv+1:4*nv);
cAf=zeros(M,N,length(t));
cBf=cAf;
cCf=cAf;
Tf=cAf;
for j=2:M-1
for i=2:N-1
I=(N-2)*(j-2)+(i-1);
for k=1:length(t)
cAf(j,i,k)=cA(k,I);
cBf(j,i,k)=cB(k,I);
cCf(j,i,k)=cC(k,I);
Tf(j,i,k)=T(k,I);
end
end
end
cAf(1,:,:)=cAf(2,:,:);
cAf(M,:,:)=cAf(M-1,:,:);
cAf(:,N,:)=cAf(:,N-1,:);
cAf(:,1,:)=cA0;
cAL=mean(squeeze(cAf(:,N,:)),1);
Xf=100.*(1-cAf./cA0);
XL=mean(squeeze(Xf(:,N,:)),1);
cBf(1,:,:)=cBf(2,:,:);
cBf(M,:,:)=cBf(M-1,:,:);
cBf(:,N,:)=cBf(:,N-1,:);
cBf(:,1,:)=cB0;
cBL=mean(squeeze(cBf(:,N,:)),1);
cCf(1,:,:)=cCf(2,:,:);
cCf(M,:,:)=cCf(M-1,:,:);
cCf(:,N,:)=cCf(:,N-1,:);
cCf(:,1,:)=cC0;
cCL=mean(squeeze(cCf(:,N,:)),1);
Tf(1,:,:)=Tf(2,:,:);
Tf(M,:,:)=(Tf(M-1,:,:)+Nu*deta*T0)./(1+Nu*deta);
Tf(:,N,:)=Tf(:,N-1,:);
Tf(:,1,:)=T0;
TL=mean(squeeze(Tf(:,N,:)),1);

Answers (1)

Cris LaPierre
Cris LaPierre on 5 May 2021
Check the dimensions of h, your output from Julia_Prob1a_fun. The ode solver expects a column vector, which has dimensions mx1. I suspect h is a matrix.

Community Treasure Hunt

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

Start Hunting!

Translated by