Solve simultaneous differential equations with embedded functions and a parameter estimation

The aim is to solve the below equations and plot $$m$$ with time, i.e. $$\frac{dm}{dt}$$

$$k$$ is unknown and needs to be estimated. For the parameter estimation, the below values in the table for m versus t can be used. This is a dummy dataset, as I would like to get an idea if this problem is solvable. If needed, an initial guesstimate for $$k$$ can be provided $$6\times10^{-5} \mathrm{m^2/(bar \cdot hr)}$$.

I have added an initial idea for the code, but not sure if my implementation of the $$A(H)$$ and $$P(H)$$ functions in the odes is correct. The code as it is doesn't run. Also I need a way to estimate for $$k$$ whilst fitting to the provided $$m$$ vs $$t$$ data.

1. Surface area in $$\mathrm{m}^2$$: $$A(H)=\pi r^2+\frac{2}{r}(1+\frac{\rho_0}{\rho_H}H)\frac{w_0}{\rho_0}$$ where $$A(H)_{t=0}=2\pi rL+2\pi r^2$$. Note that $$H$$ is differentiated later in the equation $$\frac{dH}{dt}$$.
2. Pressure in bar: $$P(H)=5004.6H+150.3$$
3. Volumetric flowrate in $$\frac{\mathrm{m^3}}{\mathrm{hr}}$$: $$\frac{dV}{dt}=\frac{k}{h}A(H)P(H)$$
4. Hydration in $$\mathrm{hr}^{-1}$$ $$\frac{dH}{dt}=\frac{\rho_H}{w_0}\frac{k}{h}A(H)P(H)$$ where $$H=\frac{w_0}{w_h}$$
5. Mass flowrate in $$\frac{\mathrm{kg}}{\mathrm{hr}}$$ $$\frac{dm}{dt}=C_d\frac{dV}{dt}$$

Where

• $$L=0.106 \mathrm{m}$$
• $$D=0.00264 \mathrm{m}$$, $$r=D/2$$
• $$\rho_0=1320\frac{\mathrm{kg}}{\mathrm{m}^3}$$ (material 1)
• $$\rho_h=1000\frac{\mathrm{kg}}{\mathrm{m}^3}$$ (water)
• $$w_0(t=0)=6.48\times10^{-4}$$ (decreases with time)
• $$w_H(t=0)=0$$ (increases with time)
• $$L=0.0015 \mathrm{m}$$
• $$C_d=300 \frac{\mathrm{kg}}{\mathrm{m}^3}$$
Time (hr) Release Rate (kg/hr)
2 .000012
4 .000018
6 .00002
8 .000015
10 .00002
12 .000015
14 .0000075
16 .000001
18 .000002
20 .000001
22 0
24 0
function dydt = diffunTAR(~,y,k, h, A_H, P_H, rho_h, rho_o, wo)
dydt = zeros (2,1);
D = 2.64/1000 %mm/1000 to convert to m
r=D/2
rho_o = 1320 %kg/m3
rho_h = 1000 % kg/m3
wo=648/1000000 %weight component 1 in mg converted to kg
h= 1.5/1000 %1.5mm expressed in m
function A_H= pi*(D/2)^2 + (2/r)*(1+(rho_o/rho_h*y(2)))*(wo/rho_o);
end
function P_H = 5004.6*y(2) + 150.39;
end
%dV/dt
dydt(1)= (k/h)* function A_H* function P_H;
%dH/dt
dydt(2)= (rho_h*k/wo*h)*function A_H*function P_H;
end

tspan = linspace(0,21*24); %time in hrs for 21 days
y0 = [0 6.4800e-04];
[t,y]=ode45(@(t,y) diffunTAR(t,y, k, h, A_H0, rho_h, rho_o, wo), tspan, y0);


Edit I realised that my previous syntax had errors, the below might be improved, but still gives me an error, displayed at the end.

function dydt = diffunTAR(~,y,k,h,rho_h,rho_o,wo,L,r)
dydt = zeros (2,1);
%dV/dt
dydt(1)= k/h* AH * PH;
%dH/dt
dydt(2)= (rho_h*k/wo*h) * AH * PH;
AH
PH
function PH
(5004.6*y(2))+150.3;
end
function AH
pi*(r^2)+(2/r)*(1+(rho_o/rho_h*y(2)))*(wo/rho_o);
setInitialConditions (AH, 2*pi*(r)*L+2*pi*(r)^2);
end
end

D = 2.64/1000 %mm/1000 to convert to m
r=D/2
L = 10.6/100 %cm/100 to convert to m
wo=648/1000000 %weight osmotic tablets in mg converted to kg
k=6e-5 %guess for permeability
h= 1.5/1000 %1.5mm expressed in m
rho_o = 1320 %kg/m3
rho_h = 1000 % kg/m3
tspan = linspace(0,21*24); %time in hrs for 21 days
y0 = [0 6.4800e-04];
[t,y]=ode45(@(t,y) diffunTAR(t,y,k,h,rho_h,rho_o,wo,L,r), tspan, y0);


Error

Error using diffunTAR/AH Too many output arguments.

Error in diffunTAR (line 4) dydt(1)= k/h* AH * PH;

Error in FitODETAR>@(t,y)diffunTAR(t,y,k,h,rho_h,rho_o,wo,L,r)

Error in odearguments (line 90) f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 115) odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

• Reading your question, you seem to have all input (except initial conditions) you need to solve the system of coupled odes. Can you clarify what difficulties you encounter or what is unclear to you?
– Bort
Nov 22 at 10:36
• I submitted an edit to convert your images into formatted text. If there were any errors in my transcription, see if you can use this same format to correct them. Nov 22 at 18:40
• @Tyberius Nice thanks! Nov 22 at 19:39
• I can't figure out from the post what the actual question is. Have you tried your code? What happens when you run it? Are the results reasonable, or nonsensical? In other words, what you have you tried already, and what is your specific question? Nov 22 at 19:42
• Start with the first things. You have a bunch of syntax errors that MATLAB complains about. That should be relatively easy to fix, but this isn't the forum to ask for general programming questions. Once you have a code that at least compiles properly, update the question with what your next problem is. Nov 22 at 21:14