0

I am trying to solve a similar problem to this one: Solving Differential Equations in Matlab

However, this time the scenario is not injection of a drug into the subcutaneous tissue and its subsequent dissolution, but a more simple situation where the suspension is allowed to dissolve in a dissolution bath of volume 900 ml.

function dydt=odefcnNY_v12(t,y,D,Cs,rho,r0,N,V)
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))/V; %dC*/dt
end

i.e. the absorption term from the previous question is removed: 

Absorption term: Af*y(2)

The compound is also different, so the MW, Cs and r0 are different, and the experimental setup is also different so W and V are now changed. To allow for these changes, the ode113 call changes to this:

MW=336.43; % molecular weight
D=9.916e-5*(MW^-0.4569)*60/600000 %m2/s - [D(cm2/min)=9.916e-5*(MW^-0.4569)*60], divide by 600,000 to convert to m2/s 
rho=1300; %kg/m3
r0=9.75e-8; %m dv50
Cs=0.032; %kg/m3 
V=0.0009;%m3 900 mL dissolution bath
W=18e-6; %kg 18mg
N=W/((4/3)*pi*r0^3*rho); % particle number
tspan=[0 200*3600]; %s in 200 hours
y0=[1 0];
[t,y]=ode113(@(t,y) odefcnNY_v12(t,y,D,Cs,rho,r0,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*');

The current problem is that this code has been running for 3 hours, and still not complete. What is different now to the previous question in the link above, that is making it take so long?

Thanks

enter image description here

Engineer_1331
  • 15
  • 4
  • I used the previous code, setting Af to zero in the definition of the new constants. No problem in the integration, using an implicit Radau method. Complete dissolution and stabilization at C*=0.625 inside the first second. You really should use `odeset` to define the absolute and relative tolerance explicitly. If I set them to high, the solver indeed takes a lot of time. The default tolerances should be low enough, but check this or set them. Abstol to 1e-6 and RelTol to 1e-9 or smaller. – Lutz Lehmann Nov 29 '19 at 09:40
  • I changed my code to this `options = odeset('RelTol',1e-6,'AbsTol',1e-9); [t,y]=ode113(@(t,y) odefcnNY_v11(t,y,D,Cs,rho,r0,N,V,Af), tspan, y0, options);` and it still takes a long time to run. I paused it to look at the errors, I attach a screenshot in the question. – Engineer_1331 Nov 29 '19 at 10:46
  • Did you try any of a) reduce the integration interval to 1 sec b) change to a different solver, c) introduce a global counter variable and print out the times for every 100th call to the ODE function and the 10 following (or do a full log) to get some progress reporting and find out where the integrator is stalling, or other debugging measures? – Lutz Lehmann Nov 29 '19 at 11:13
  • We are probably horribly off-topic. Your code is working without errors, it is just that the results are not as expected. Or that the numerical solver stalls for this specific ODE. Which is more to do with the ODE than the solver. Such discussions are better located at scicomp.SE. – Lutz Lehmann Nov 29 '19 at 11:59

1 Answers1

0

I can not really reproduce your problem. I use the "standard" python modules numpy and scipy, copied the block of parameters, 

MW=336.43; # molecular weight
D=9.916e-5*(MW**-0.4569)*60/600000 #m2/s - [D(cm2/min)=9.916e-5*(MW^-0.4569)*60], divide by 600,000 to convert to m2/s 
rho=1300.; #kg/m3
r0=9.75e-8; #m dv50
Cs=0.032; #kg/m3 
V=0.0009;#m3 900 mL dissolution bath
W=18e-6; #kg 18mg
N=W/((4./3)*pi*r0**3*rho); # particle number
Af = 0; # bath is isolated

used the same ODE function like in the previous post (remember Af=0)

def odefcnNY(t,y,D,Cs,rho,r0,N,V,Af):
    r,C = y;
    drdt = (-D*Cs)/(rho*r0**2)*(1-C) * r/(1e-10+r**2); # dr*/dt
    dCdt = (D*4*pi*N*r0*(1-C)*r-(Af*C))/V;            # dC*/dt
    return [ drdt, dCdt ];

and solved the ODE 

tspan=[0, 1.0]; #1 sec
#tspan=[0, 200*3600]; #s in 200 hours
y0=[1.0, 0.0];
method="Radau"

sol=solve_ivp(lambda t,y: odefcnNY(t,y,D,Cs,rho,r0,N,V,Af), tspan, y0, method=method, atol=1e-8, rtol=1e-11);

t = sol.t; r,C = sol.y;
print(sol.message)
print("C*=",C[-1])

This works in a snap, using 235 steps for the first second and 6 further steps to cover the constant behavior in the remaining time of the 200 hours.

plots of the components

I can only mess it up by increasing the tolerances to unreasonably large values like 1e-4, and only if the epsilon used in the mollification is 1e-12. Then the hard turn when the radius reaches zero is too hard, the step size controller falls into a loop. This is more an error of the crude implementation of the step size controller, that should not be the case in the Matlab routines.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51