loop with changing size and variable values

1 view (last 30 days)
H.O.Z.
H.O.Z. on 17 May 2015
Commented: H.O.Z. on 17 May 2015
Hi, I'm pretty new to MATLAB. I need to build a model to simulate molecular flux within a multi-compartment system using the ode15s function. Below is my main script.
dp0 = 2e-7;
nL0 = 5;
r0=zeros(nL0,1);
V0=zeros(nL0,1);
C0_bulk=zeros(nL0,1);
for i = 1:nL0,
r0(i)=dp0/2/nL0;
V0(i)=pi/6*((dp0-2*r0(i)*(i-1))^3-(2*(nL0-i)*r0(i))^3);
C0_bulk(i)=0.797e6*V0(i)*1e11/338*6.02e23/1e6;
end
c0 = [0;C0_bulk];
options = odeset('RelTol',1e-10);
[t,c] = ode15s(@myfun,[0:120],c0,options,nL0,V0,C0_bulk);
plot(t,40.9*4.065e-14*338*c(:,1),'r.',
t,40.9*4.065e-14*338*(c(:,2)+c(:,3)+c(:,4)+c(:,5)+c(:,6)),'b.');
Below is my function:
function [dcdt, dp] = myfun(t,c,nL0,V0,C0_bulk)
nL=nL0;
V=zeros(nL0,1);
r=zeros(nL0,1);
SA=zeros(nL0,1);
c = zeros(nL0+1,1);
for i = 1:nL0,
V(i) = V0(i)*c(i+1)/C0_bulk(i);
end
dp = (6*sum(V)/pi)^(1/3);
for i = 1:nL0-1,
r(i) = (6*sum(V(i:nL))/pi)^(1/3)/2-(6*sum(V(i+1:nL))/pi)^(1/3)/2;
end
r(nL)=(6*V(nL)/pi)^(1/3)/2;
SA(1)= pi*dp^2;
for i = 2:nL0,
SA(i)=pi*(dp-2*sum(r(1:i-1)))^2;
end
if r(1)<5e-9;
r(1)=r(2)+r(1);
V(1)=V(2)+V(1);
c(2)=c(3)+c(2);
for i = 2:nL0-1,
r(i)=r(i+1);
V(i)=V(i+1);
c(i)=c(i+1);
end
r(nL0)=0;
V(nL0)=0;
c(nL0)=c(nL0+1);
c(nL0+1)=0;
nL = nL-1;
end
dcdt(1) = 1.5e12-c(1); % concentration in the gas (out of the compartment system)
dcdt(2) = -(1.5e12-c(1)) +...
(c(3)/V(2)-c(2)/V(1))*8e-5/pi/(r(1)+r(2))*SA(2);
for i = 3:nL,
dcdt(i) = (-c(i)/V(i-1)+c(i-1)/V(i-2))*8e-5/pi/(r(i-2)+r(i-1))*SA(i-1)+...
(c(i+1)/V(i)-c(i)/V(i-1))*8e-5/pi/(r(i-1)+r(i))*SA(i);
end
dcdt(nL+1) = (-c(nL+1)/V(nL)+c(nL)/V(nL-1))*8e-5/pi/(r(nL-1)+r(nL))*SA(nL);
if nL<nL0,
dcdt(nL0+1) = 0;
end
dcdt = dcdt(:);
In this code, the number of compartments ('nL0' as initial value that decreases upon the size of the first compartment ('nL')) changes. A force drives molecules transport from the surface compartment towards the outside region. If the size of the surface compartment (depending upon the remaining volume of molecules) decreases to certain level, this compartment gets mixed with the next compartment and the two sized add up and becomes the size of the new compartment. Then this could happen to the next compartment on and on.
I don't think I could easily change the matrix size of my output (c). So I tried to maintain the total size, but when the above transfer happens, I let the compartments shift towards the surface and leave the last compartment to zero. I am not sure if my code does this right.
My code runs with errors "Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN." and ends up with "Warning: Failure at t=0.000000e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.905050e-323) at time t. " Please help me solve my problem. Thanks a lot.
  1 Comment
Jan
Jan on 17 May 2015
Please post a copy of the complete error message. Most of all it is timeconsuming to guess, which line causes the error.

Sign in to comment.

Answers (1)

Jan
Jan on 17 May 2015
Edited: Jan on 17 May 2015
This looks like a discontinuity:
if nL<nL0,
dcdt(nL0+1) = 0;
end
Matlab's integrators are not designed to handle non-smooth functions. See: http://www.mathworks.com/matlabcentral/answers/59582#answer_72047 . It is not clear, if this causes the problem you observe currently, but it is at least serious.
Note: The code will look much leaner, when you use vectorized code instead of loops. Compare:
r(1)=r(2)+r(1);
V(1)=V(2)+V(1);
c(2)=c(3)+c(2);
for i = 2:nL0-1,
r(i)=r(i+1);
V(i)=V(i+1);
c(i)=c(i+1);
end
r(nL0)=0;
V(nL0)=0;
c(nL0)=c(nL0+1);
c(nL0+1)=0;
with:
r = [r(1)+r(2), r(3:nL0), 0];
V = [V(1)+V(2), V(3:nL0), 0];
c = [c(2)+c(3), c(3:nL0+1), 0];
  1 Comment
H.O.Z.
H.O.Z. on 17 May 2015
Thanks Jan, I think the discontinuity is the problem. Any idea for solutions?

Sign in to comment.

Categories

Find more on Programming in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!