# FindMinimum problem (edited question)

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: 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[]]], model, {a, b, c}, t]
modelf = Function[{t}, Evaluate[model /. fit]];
Plot[modelf[t], {t, 0, 200},
Epilog -> Map[Point, Rest[Rest[data[]]]]] • Thanks for your kind reply and time. – amirtech2005 Oct 17 '16 at 19:31
• But I feel you have extremely simplified the question I asked. I have to use the function I defined for din. – amirtech2005 Oct 17 '16 at 19:34
• Your question explicitly asked for "something better to fit the tail of my actual data." Try a Log! – bill s Oct 17 '16 at 19:42
• 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. – JimB Oct 17 '16 at 19:45
• 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! – amirtech2005 Oct 17 '16 at 19:50

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]] 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}]] 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"}] 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.

• 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. – amirtech2005 Oct 17 '16 at 20:17
• 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. – JimB Oct 17 '16 at 20:41