MATLAB: Iteratively solve IVP with ode45

functionodeode45

I have a problem where I'm trying to solve an ODE dh/dt = (4/pi*D_t)*(Qin – Qout). Qout(t) is related to h(t), but cannot be manipulated to give an explicit expression for Qout in terms of h. I created a function file shown below to define the relationship between Qout and h:
function x = flowratea(h)
L = 30;
D_p = 1/12;
rho = 5.374e-4;
mu = 3.367e-7;
g = 1.159e+5;
x = @(Q) h - (1 + (0.25/(log(5.74/(4*rho*Q/(pi*mu*D_p))^0.9)^2))*L/D_p)*(4*Q/(pi*D_p^2))^2/(2*g);
I also created the following function file for dh/ht:
function dh = dheight(t,D_t,Qin,Qout)
dh = (4/pi*D_t)*(Qin - Qout);
Because I need to solve for Qout(t) to evaluate dh/dt, and I need h(t) to solve for Qout(t), I am trying to use a loop to use the initial h to solve Qout(t0), then dh/dt and then h(t1) and so on. I've been messing around with this for a while without any luck. Below is my current code:
for i = 1:length(time)
    t = i;
    Q = fzero(flowratea(h),40);
    dh = dheight(t,D_t,Qin,Q);
    [tn(i),hn(i)] = ode45(dh,time(i),h);
    h = hn(i);
end
The error I'm getting with this is:
??? Error using ==> exist The first input to exist is a string.
Error in ==> odearguments at 79 if (exist(ode)==2) % M-file
Error in ==> ode45 at 173 [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, …
Error in ==> TankDrain at 23 [tn(i),hn(i)] = ode45(dh,time(i),h);
I've attempted to call the dheight function file in the ode45, but then it says I have too many inputs. I've tried a number other things I can't remember exactly, and I'm beginning to wonder if my approach is completely off.

Best Answer

function main
h0=...;
Qout0=...;
x0=[h0 Qout0];
M=[1 0
   0 0];
options=odeset('Mass',M);
[T Y]=ode15s(@flowrates,[0 10],x0,options);
function dy=flowrates(t,y) 
L = 30;
D_p = 1/12;
rho = 5.374e-4;
mu = 3.367e-7;
g = 1.159e+5;
Qin=...;
h=y(1);
Qout=y(2);
dy=zeros(2,1);
dy(1)=(4/pi*D_t)*(Qin - Qout);
dy(2)=h - (1 + (0.25/(log(5.74/(4*rho*Qout/(pi*mu*D_p))^0.9)^2))*L/D_p)*(4*Qout/(pi*D_p^2))^2/(2*g);
Best wishes
Torsten.
Related Question