ODE45 wont run, just says "error"?

1 Ansicht (letzte 30 Tage)
Nathan Corbissero
Nathan Corbissero am 16 Sep. 2019
Kommentiert: Star Strider am 16 Sep. 2019
Hi all, first post here so sorry If I am leaving out details. I am writing an ODE45 file which contains a nested rate function. I've done this numerous times before, but I am getting an error which does not specify where I should look to debug.
Code:
function res = renalclearance()
t_final = 120; %minutes
[T,M] = ode45(@renal_rate, [0,t_final], 10); %10 = initial concentration (mmol/L)
Y = M(:,1);
W = M(:,2);
figure(1) %dci/dt
plot(T,Y)
figure(2) %dce/dt
plot(T,W)
end
function res = renal_rate(t,Z)
Ci = Z(1);
Ce = Z(2);
kc = 0.77; %L/min
kt = 0.25; %L/min
qf = 0.083; %L/min
gi = 0.1; %mmmol/min
v0 = 40; %L
vr = 1.67;
vt = v0-(t*qf);
vi = ((vt)/(1+(1/vr)));
ve = ((vt)/(1+vr));
vdi = 1/(1.598);
vde = 1/(2.67);
dci_dt = -((vdi+kc)*(Ci/vi))+((kc*Ce + gi)/(vi));
dce_dt = -((vde+kc+kt)*(Ce/ve))+((kc*Ci)/(ve));
res = [dci_dt,dce_dt];
end
And here is the corresponding error message:
Error in renalclearance (line 5)
[T,M] = ode45(@renal_rate, [0,t_final], 10); %10 = initial concentration (mmol/L)

Akzeptierte Antwort

Star Strider
Star Strider am 16 Sep. 2019
You have two problems:
First, your system has two differential equations, so you must have two initial conditions:
ic = [0 10];
[T,M] = ode45(@renal_rate, [0,t_final], ic); %10 = initial concentration (mmol/L)
Second, ‘renal_rate’ has to return a column vector:
res = [dci_dt; dce_dt];
Note the vertical concatenation operator here, (;). With these changes your code runs without error. Change the initial conditions (the ‘ic’ vector) if I guessed wrong as to what they should be.
  2 Kommentare
Nathan Corbissero
Nathan Corbissero am 16 Sep. 2019
That worked perfectly! Thank you so much.
Star Strider
Star Strider am 16 Sep. 2019
As always, my pleasure!

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Kategorien

Mehr zu Programming finden Sie in Help Center und File Exchange

Tags

Produkte

Community Treasure Hunt

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

Start Hunting!

Translated by