0
$\begingroup$

I am doing the following analysis to find a RHfit function from "din function" that I have defined in the beginning of the code. To verify if RHfit is great or not to use in my future work, I check if mass loss fit is good or not.

I appreciate if someone can help... I have attached the required mass-loss data:

https://www.dropbox.com/s/d5ywe17ye7300s5/Masslossdata.xlsx?dl=0

Code:

datafile = Import["Masslossdata.xlsx"];
massdata = datafile[[1, 2 ;; All, All]];
RHmax = 96.6438;
RHbound = 70;
L = .5/39.4;
l = 1;
b = 2.4/39.4;
ϕ = .27;
ρ = 1000;
V = l L b;
Sfit2[RH_, y_] = 
    (1 - Log[0.0103473 RH]/0.712045)^(-1.00096 (1 - 78.8 y * -0.110095));

α0 = 0.3; n = 4; He = 30;

din[di0_, RH_] = di0/100000  
  (α0 + (1 - α0)/(1 + ((100 - RH[x, t])/(100 - He))^n))
tdata = massdata[[All, 1]];
RHboundNL = RHmax - (RHmax - RHbound ) Tanh[t*1000];
eq1[di0_] = D[RH[x, t], t] - D[din[di0, RH] D[RH[x, t], x], x] == 0
IC1 = RH[x, 0] == RHmax; 
BC1 = RH[L/2, t] == RHboundNL; 
BC2 = (D[RH[x, t], x] /. x -> -L/2) == 0;
MaxTime = massdata[[Length[massdata], 1]];
fun[data_, deq_] := Module[{sol2, residuals},
   sol2 = NDSolve[{deq, IC1, BC1, BC2},RH, {t, 0, MaxTime}, {x, -L/2, L/2}];
   RHfit[t_, x_] = Evaluate[RH[x, t] /. sol2];
   RHavg[t_] := 1/L Integrate[RHfit[t, x], {x, -L/2, L/2}];
   Savg[t_] := 1/L Integrate[Sfit2[RHavg[t], y], {y, -L/2, L/2}];
   Massloss[t_] := -(Savg[t] - Savg[0.001]) ρ ϕ V*1000;
   modeldata = Map[Evaluate[Massloss], tdata][[All, 1]];
   residuals = N[data[[All, 2]] - modeldata];
   residualssquare = Total[residuals^2]];
fun1[di0_?NumberQ] := fun[massdata, eq1[di0]]
{res, dii} = 
 Reap[FindMinimum[{fun1[di0], .1 < di0 < 2}, {di0, .25}, 
   MaxIterations -> 11, AccuracyGoal -> 4, PrecisionGoal -> 3, 
   WorkingPrecision -> 3, Method -> InteriorPoint, 
   StepMonitor :> Sow[di0]]]
di0 = Last[Last[dii]]
sol2 = NDSolve[{eq1[di0], IC1, BC1, BC2},RH, {t, 0, MaxTime}, {x, -L/2, L/2}];
RHfit[t_, x_] = Evaluate[RH[x, t] /. sol2];
RHavg[t_] := 1/L Integrate[RHfit[t, x], {x, -L/2, L/2}];
Savg[t_] := 1/L Integrate[Sfit2[RHavg[t], y], {y, -L/2, L/2}];
Massloss[t_] := -(Savg[t] - Savg[0.001]) ρ ϕ V*1000;
modeldata = Map[Evaluate[Massloss], tdata][[All, 1]];
residuals = N[massdata[[All, 2]] - modeldata];
residualssquare = Total[residuals^2]
Plot[Savg[t], {t, 0, MaxTime}, PlotRange -> {0, 1}]
Show[
  Plot[Massloss[t], {t, 0, MaxTime}, 
    PlotRange -> {0, 200},
    Frame -> True, 
    FrameLabel -> {"Time (d)", "Mass loss (g)"}],
  ListPlot[massdata]]

I think the problem should be with one of these, which I couldn't fix yet:

1- I have tried many ranges for di0 and the initial guess like the one I showed in the code:

FindMinimum[{fun1[di0], .1 < di0 < 2}, {di0, .25}

2- I have also changed these parameters, but I think they might be ok:

α0 = 0.3; n = 4; He = 30;

The fit I got for different runs is not better than this, but I need something better to fit the tail of my actual mass data as well; after that I can trust my RHfit function for my future use:

mass loss

Thanks in advance!

$\endgroup$
1
$\begingroup$

When trying to fit data, it often helps to remain flexible as to the functional form. Using your data and a Log model gives a better fit at the larger values.

model = a Log[t + b] + c;
fit = FindFit[Rest[Rest[data[[1]]]], model, {a, b, c}, t]
modelf = Function[{t}, Evaluate[model /. fit]];
Plot[modelf[t], {t, 0, 200}, 
 Epilog -> Map[Point, Rest[Rest[data[[1]]]]]]

enter image description here

$\endgroup$
  • $\begingroup$ Thanks for your kind reply and time. $\endgroup$ – amirtech2005 Oct 17 '16 at 19:31
  • $\begingroup$ But I feel you have extremely simplified the question I asked. I have to use the function I defined for din. $\endgroup$ – amirtech2005 Oct 17 '16 at 19:34
  • $\begingroup$ Your question explicitly asked for "something better to fit the tail of my actual data." Try a Log! $\endgroup$ – bill s Oct 17 '16 at 19:42
  • $\begingroup$ I feel your function just can't fit the data. (This is similar to a previous question of yours.) Real data has a habit of getting in the way of theory. The lack of fit of your original function should be convincing to whoever is telling you to keep that function and they should change their minds. $\endgroup$ – JimB Oct 17 '16 at 19:45
  • $\begingroup$ Yes, I understand. I have been trying to use another model for "din", but I have to do minor changes... because, as I said in my question, I am trying to find RHfit function and check if RHfit is ok through the mass loss fit. I just don't know how to do the minor changes in din! $\endgroup$ – amirtech2005 Oct 17 '16 at 19:50
1
$\begingroup$

To follow @bills 's suggestion, you might want to take the logs of both the dependent and independent variable (but getting rid of that first data point as {0,0} is many times not a real data point) and add in the square of the log of the independent variable.

data = {{1, 21.56667}, {3, 30.56667}, {4, 33.03333}, {6, 
    41.36667}, {7, 43.26667}, {8, 45.73333}, {9, 46.83333}, {10, 
    47.83333}, {12, 50.56667}, {15, 53.36667}, {18, 55.86667}, {22, 
    58.83333}, {26, 61.4}, {37, 66.66667}, {45, 69.26667}, {57, 
    71.63333}, {71, 75.2}, {80, 77.56667}, {90, 79.06667}, {113, 
    83.4}, {146, 87.1}, {176, 89.7}};
nlm = NonlinearModelFit[Log[data], a + b logx + c logx^2, {a, b, c}, logx];
Show[ListPlot[data], Plot[Exp[nlm[Log[x]]], {x, 1, 200}, PlotRange -> All]]

Fit on original scale

Sometimes looking at the fit on the log scale highlights both goodness of fit and lack of fit:

Show[ListLogLogPlot[data], LogLogPlot[Exp[nlm[Log[x]]], {x, 1, 200}]]

Fit on log-log scale

Looking at the residuals shows there's still something more going on than what the model can handle:

ListPlot[Transpose[{nlm["PredictedResponse"], data[[All, 2]] - Exp[nlm["PredictedResponse"]]}],
 Frame -> True, FrameLabel -> {"Predicted", "Residual"}]

Residuals

In short, your data is not well described by your original function and models with one or two more parameters still don't provide a fit with desirable residuals.

$\endgroup$
  • $\begingroup$ Thanks for your time. But the problem is somewhere in din... fitting mass loss itself is not my problem, I only use that fit to double check the goodness of RH function which is found through RHfit. $\endgroup$ – amirtech2005 Oct 17 '16 at 20:17
  • $\begingroup$ OK. You might consider rewriting your question to emphasize and clarify that. One of the reasons I simply added one more term to @bills 's answer was that I couldn't make heads or tails out of your code. $\endgroup$ – JimB Oct 17 '16 at 20:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.