Mathematica running out of memory when fitting parameters of system of differential equations

Question

I encountered a memory problem in Mathematica when I tried to process my experimental data. I'm using Mathematica to find the optimal parameters for a system of three partial differential equations.

When the e parameter was greater than 0.4, Mathematica consumed a lot of memory. For e < 0.4, the program worked properly.

I have tried using $HistoryLength = 0, and reducing AccuracyGoal and WorkingPrecision with no success.

I'm trying to understand what mistakes I made, and how I may limit the memory usage.

Clear[T, L, e, v, q, C0, R, data, model];
T = 13200; 
L = 0.085; 
e = 0.41; 
v = 0.000557197; 
q = 0.1618; 
C0 = 0.0256; 
R = 0.00075;

data = {{L, 600, 0.141124587}, {L, 1200, 0.254134509}, {L, 1800, 
    0.342888644}, {L, 2400, 0.424476295}, {L, 3600, 0.562844542}, {L, 
    4800, 0.657111356}, {L, 6000, 0.75137817}, 
       {L, 7200, 0.815876516}, {L, 8430, 0.879823594}, {L, 9000, 
    0.900771775}, {L, 13200, 1}};

model[(De_)?NumberQ, (Kf_)?NumberQ, (Y_)?NumberQ] := 
 model[De, Kf, Y] =  yeld /. Last[Last[
     NDSolve[{
       v D[Ci[z, t], z] + D[Ci[z, t], t] == -((
         3 (1 - e) Kf (Ci[z, t] - C0))/(
         R e (1 - (R Kf (1 - R/r[z, t]))/De))),
       D[r[z, t], t] == (R^2 Kf (Ci[z, t] - C0))/(
        q r[z, t]^2 (1 - (R Kf (1 - R/r[z, t]))/De)),
       D[yeld[z, t], t] == Y*(v e Ci[z, t])/(L q (1 - e)),
       r[z, 0] == R,
       Ci[z, 0] == 0,
       Ci[0, t] == 0,
       yeld[z, 0] == 0},
      {r[z, t], Ci[z, t], yeld}, {z, 0, L}, {t, 0, T}]]]

fit = FindFit[
  data, {model[De, Kf, Y][z, t], {0.97  < Y < 1.03, 
    10^-6 < Kf < 10^-4, 10^-13 < De < 10^-9}}, 
     {{De, 10^-12}, {Kf,  10^-6}, {Y, 1}}, {z, t}, Method -> NMinimize]

data = {{600, 0.141124587}, {1200, 0.254134509}, {1800, 
    0.342888644}, {2400, 0.424476295}, {3600, 0.562844542}, {4800, 
    0.657111356}, {6000, 0.75137817}, {7200, 0.815876516}, 
       {8430, 0.879823594}, {9000, 0.900771775}, {13200, 1}}; 

YYY = model[De /. fit[[1]], Kf /. fit[[2]], Y /. fit[[3]]]; 

Show[Plot[Evaluate[YYY[L, t]], {t, 0, T}, PlotRange -> All], 
 ListPlot[data, PlotStyle -> Directive[PointSize[Medium], Red]]]

Link to the .nb file: http://www.4shared.com/folder/249TSjlz/_online.html

Answer 1

I have a sneaking suspicion the reason why it's failing on you is because you're caching the results.

Do you need to store every single solution that NDSolve is producing? I'm somewhat skeptical of whether this is useful or not for Findfit, since I highly doubt it would revisit a past result.

Besides, it's not like you're talking about integers where there's a finite domain you're talking about. You're using reals and even over the range you specify, there's A LOT of different solutions possible. I don't think you want to store each and every one of them.

Rewrite your code so that instead of having:

model[(De_)?NumberQ, (Kf_)?NumberQ, (Y_)?NumberQ] := 
 model[De, Kf, Y] =  yeld /. Last[Last[
     NDSolve[..]

You instead have:

model[(De_)?NumberQ, (Kf_)?NumberQ, (Y_)?NumberQ] := 
 NDSolve[..]

By caching your previous results, you're going to eat up memory like no tomorrow with FindFit. Typically it's useful when you have a recurrence relation, but here I'd seriously advise against it.

---

Some Notes:

After running for 2415 seconds on my machine, Mathematica's memory usage went from 112,475,400 bytes to 1,642,280,320 bytes with the caching.

I'm currently running the code without the caching now.