loop with changing size and variable values

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.

[Jan]

Please post a copy of the complete error message. Most of all it is timeconsuming to guess, which line causes the error.

Answer 1

[Jan]

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];

[H.O.Z.]

Thanks Jan, I think the discontinuity is the problem. Any idea for solutions?