Simultaneous differential equations - derivative of y^3 wrt t

I have been trying to solve this problem for sometime and so far I have received an answer that requires the symbolic toolbox which I dont have - I would appreciate any help!

I need to solve these 2 differential equations simultaneously.

But I dont know how to code the dr*^3/dt. My code is below:

function dydt=odefcnNY_v3(t,y,D,Cs,rho,r0,N,V,Af)
dydt=zeros(2,1);
dydt(1)=(-3*D*Cs/rho*r0^2)*y(1)*(1-y(2));
dydt(2)=(D*4*pi*N*r0*(1-y(2))*y(1)-(Af*y(2)))/V;
end

y(1) = r* and

y(2) = C*

So

dydt(1) = dr*/dt and

dydt(2) = dC*/dt

In my case dydt(1) needs to be replaced with something that would solve dr*^3/dt and not dr*/dt. The rest of the code is below.

D=4e-9;%m2/s
rho=1300; %kg/m3
r0=10.1e-6; %m dv50
Cs=0.0016; %kg/m3
V=1.5e-6;%m3
W=4.5e-8; %kg
N=W/(4/3*pi*r0^3*rho);
Af=0.7e-6/60; %m3/s
tspan=[0 24*3600]; %s in 24 hrs
y0=[r0 0];
[t,y]=ode45(@(t,y) odefcnNY_v3(t,y,D,Cs,rho,r0,Af,N,V), tspan, y0);
plot(t/3600,y(:,1),'-o') %plot time in hr, and r*
xlabel('time, hr')
ylabel('radius,um')
legend('DCU')
plot(t/3600,y(:,2),'-') %plot time in hr, and C*
xlabel('time,hr')
ylabel('C* (C/Cs)')
legend('DCU')

[Nora Rafael]

So I got an idea to simplify the left hand side of the first ode from dr^3/dt to 3*r^2*(dr/dt).

Dividing both sides of the first equation by 3*r^2 gives:

dydt(1)=((-D*Cs)/(rho*r0^2*y(1)))*(1-y(2)); % dr*/dt

It took a long time to run and when it finished I got this error:

"Error using horzcat Requested 2x528701600 (7.9GB) array exceeds maximum array size preference. Creation of arrays greater than this limit may take a long time and cause MATLAB to become unresponsive. See array size limit or preference panel for more information. Error in ode45 (line 484) yout = [yout, zeros(neq,chunk,dataType)]; "

Any idea what this means?

I also attach my files incase anyone has time to look into it.

Thanks!

[David Goodmanson]

Hi Nora,

your original expression for dydt(1) is algebraically incorrect, but you changed it later anyway.

I am speculating that r* = r/r0, in which case your starting value should be y0 = [1 0] and not y0 = [r0 0]. However that is just a guess.

With starting values [1 0] and the new expression for dydt(1) you do get results,but not sensible results for C*. Those get a lot better after fixing things to make the argument order the same in your function as in the ode45 call (Af and V are switched). Still and all, up around 12000 seconds r* collapses to zero and the solution dies.

[Nora Rafael]

Hi David

Thanks very much for looking into it.

I have made improvements since then.

Thank you for pointing out that y0 should be [1 0]!!

I applied an implicit method in view of the near singularity at r=0

then a suitable mollification of 1/r is

r/(eps^2+r^2)

with some small

eps

Function is then:

function dydt=odefcnNY_v10(t,y,D,Cs,rho,r0,N,V,Af)
dydt=zeros(2,1);
dydt(1)=(-D*Cs)/(rho*r0^2)*(1-y(2))*y(1)/(1e-6+y(1)^2); % dr*/dt
dydt(2)=(D*4*pi*N*r0*(1-y(2))*y(1)-Af*y(2))/V; %dC*/dt
end

and the rest of the code was adapted after feedback from the authors of the paper on what parameters they used

MW=224; % molecular weight
D=9.916e-5*(MW^-0.4569)*60/600000 %m2/s - [D(cm2/min)=9.916e-5*(MW^-0.4569)*60] equation provided by authors, divide by 600,000 to convert to m2/s 
rho=1300; %kg/m3
r0=10.1e-6; %m dv50
Cs=1.6*1e6/1e9; %kg/m3 - 1.6ug/m3 converted to kg/m3
V=5*0.3/1e6;%m3 5ml/Kg animal * 0.3Kg animal, divide by 1e6 to convert to m3
W=30*0.3/1000000; %kg; 30mg/Kg animal * 0.3Kg animal, divide by 1e6 to convert to m3
N=W/((4/3)*pi*r0^3*rho); % particle number
Af=0.7e-6/60; %m3/s
tspan=[0 24*3600]; %s in 24 hrs
y0=[1 0];
[t,y]=ode113(@(t,y) odefcnNY_v11(t,y,D,Cs,rho,r0,Af,N,V), tspan, y0);
plot(t/3600,y(:,1),'-o') %plot time in hr, and r*
xlabel('time, hr')
ylabel('r*, (rp/r0)')
legend('DCU')
title ('r*');
plot(t/3600,y(:,1)*r0*1e6); %plot r in microns
xlabel('time, hr');
ylabel('r, microns');
legend('DCU');
title('r');
plot(t/3600,y(:,2),'-') %plot time in hr, and C*
xlabel('time, hr')
ylabel('C* (C/Cs)')
legend('DCU')
title('C*');
plot(t/3600, y(:,2)*Cs) % time in hr, and bulk concentration on y
xlabel('time, hr')
ylabel('C, kg/m3')
legend('Dissolved drug concentration')
title ('C');

Eventually the r* graph makes sense, but not the C* as you mentioned.

And now I am stuck :) I attach the paper I am trying to reproduce the data from, C* in Figure 3, along with my updated Matlab files.

[Nora Rafael]

Hi David

You mentioned that:

>> "Those get a lot better after fixing things to make the argument order the same in your function as in the ode45 call (Af and V are switched). "

Can you tell me a bit more about what you did there?

I also switched to ode113 as ode45 was taking a while to work, maybe it is improved due to the increased step size.

Thanks

[David Goodmanson]

Hi Nora,

you have

[t,y]=ode113(@(t,y) odefcnNY_v11(t,y,D,Cs,rho,r0,Af,N,V), tspan, y0);   % and
      function dydt=odefcnNY_v10(t,y,D,Cs,rho,r0,N,V,Af)

As you can see, the arguments lists don't match. What ode115 has for Af, the function thinks is N, and vice versa. I noticed this because after implementing the

dr*^3/dt = 3r*^2 dr/dt 

fix, C* was changing by an absurdly small amount so something still had to be wrong. After matching up the argument lists, C* changes at a sensible rate.

The differential equations seem inclined to crush r* down to a point after a certain amount of time goes by. I would not expect that from a diffusion model. Could you comment on what process is being modeled?

[Nora Rafael]

Ah so the order of the arguments list have to match.

This fixes the strange numbers with C*. They are now more reasonable, plateaus at around 0.9 however, whereas in the paper they report 0.25.

The process is the dissolution of a microsuspension which ultimately results in concentration changes at the site of dissolution and a reduction in the microparticle radius. I attach some slides which explain first the Nernst-Brunner equation for dissolution, and secondly the derivation of the equations used in the model based on mass balances. I am not sure why the solution collapses, it seems to me like a mathematical rather than a physical issue. If you run the code I attached to 1000 hours (so tspan up to 3600000s) you will see that the particle radius drops to 0 at around 150 hours or 540000s. Then there is no more solid available to be dissolved and C* also drops to zero.

I have received feedback that a suitable mollification of 1/r in the final dr/dt equation is

r/(eps^2+r^2)

with some small eps (I have used 1e-6). This should trap the radius at zero without becoming negative during the numerical process. This is what I have now implemented in the odefcNY_v11 function file.

[David Goodmanson]

There appears to be a mistake on the second slide. For dC*/dt the highlighted equation has a term Af*(C_L), whereas your equation uses Af*(C*). But if you look at the derivation on that slide, Af*(C*) would be correct.

[Nora Rafael]

Hi David, that's right, the slide is wrong - it should be C*. I originally had the equations in terms of CL and later converted to C*, and didn't update this part of the equation on the slides.

[Nora Rafael]

I am now trying to implement a similar problem, but with slightly different parameters and a simplified function with no AfxC* term. This one seems to be taking hours and still no solution, whereas the previous one was taking only a few seconds.

Answer 1

[darova]

Just get

dydt(1)=(-3*D*Cs/rho/r0^2)*y(1)*(1-y(2));
dydt(1)=nthroots(dydt(1),3);