1
$\begingroup$

This is my first question in this forum. I have a relative humidity data as you can see in this excel-sheet:

RH data

I need to find the Di (drying diffusion coefficient) using FindFit, but the final fit is not very accurate as shown below. I will find something for Di but I want the tail of my graph to be fitted well too, as shown in the red marked area... Can anyone help me with this? I am sure I just need to adjust something, like matching with my data. I don't know how. The first 2 columns in the excel file is what I am working on. Please ignore the other columns.

"t" is time, L, y, RHinitial, and RHbound also cannot be changed.

My code

Code:

RHdatafile = Import["RH.xlsx"];
RHdata = RHdatafile[[1, 3 ;; All, All]];
topRH = RHdata[[All, 1 ;; 2]];
midRH = RHdata[[All, 3 ;; 4]];
botRH = RHdata[[All, 5 ;; 6]];
RHinitial = 92.8;
RHbound = 70;
L = 0.102;
y = 0.038;
RHfittop = 
  FindFit[topRH, {RHinitial + 
     RHinitial Erfc[L/Sqrt[
       di t]] + (RHbound - 
        RHinitial) (Erfc[(L - 2 y)/(4 Sqrt[di t])] + 
        Erfc[(3 L + 2 y)/(4 Sqrt[di t])]), 
    0 < di < 0.00003}, {di}, t];

di = di /. RHfittop

Show[Plot[{RHinitial + 
    RHinitial Erfc[L/Sqrt[
      di t]] + (RHbound - 
       RHinitial) (Erfc[(L - 2 0.038)/(4 Sqrt[di t])] + 
       Erfc[(3 L + 2 0.038)/(4 Sqrt[di t])])}, {t, 0, 200}, 
  PlotRange -> {40, 100}, PlotLegends -> "Expressions", 
  FrameLabel -> {"Time (d)", "RH(%)"}], 
 ListPlot[topRH, PlotRange -> All, 
  FrameLabel -> {"Time (d)", "RH(%)"}]]
$\endgroup$
  • $\begingroup$ Please post code, rather than an image of code. $\endgroup$ – Feyre Oct 6 '16 at 19:05
  • $\begingroup$ People here generally like to copy-paste a code so they can play with it, so don't post an image of a code. You can format the code with the {} button; go to help by clicking the ? button for additional details. $\endgroup$ – corey979 Oct 6 '16 at 19:06
  • $\begingroup$ Without going too much in details of your code I would expect that most probable, the model should be improved. $\endgroup$ – Alexei Boulbitch Oct 6 '16 at 19:18
  • 1
    $\begingroup$ with such a single parameter fit you can manually adjust the parameter and look for a fit ( use Manipulate ) and sometimes do better than FindFit or at least find a good initial guess for FindFit. Otherwise I think asking for a better model is outside the scope of this site. $\endgroup$ – george2079 Oct 6 '16 at 19:36
  • $\begingroup$ "t" is time, L, y, RHinitial, and RHbound also cannot be changed. $\endgroup$ – amirtech2005 Oct 6 '16 at 19:43
1
$\begingroup$

Based on the comments, the one-parameter, di, model fitted to the data is not good enough to describe the data adequately. Also, it was stated by the OP that one parameter, RHbound, might not had been measured correctly. So, let's fit a two-parameter model:

RHfittop = 
 FindFit[topRH, {RHinitial + 
    RHinitial Erfc[
      L/Sqrt[di t]] + (RHbound - 
       RHinitial) (Erfc[(L - 2 y)/(4 Sqrt[di t])] + 
       Erfc[(3 L + 2 y)/(4 Sqrt[di t])]), 0 < di < 0.00003}, {di, 
   RHbound}, t]

{di, RHbound} = {di, RHbound} /. RHfittop

{0.0000109214, 68.0971}

Both these parameters don't differ much from RHbound = 70, originally set by the OP, and di = 0.0000154358:

plot1 = Plot[{RHinitial + 
    RHinitial Erfc[
      L/Sqrt[di t]] + (RHbound - 
       RHinitial) (Erfc[(L - 2 0.038)/(4 Sqrt[di t])] + 
       Erfc[(3 L + 2 0.038)/(4 Sqrt[di t])])}, {t, 0, 200}, 
  PlotRange -> {40, 100}, Frame -> True, 
  FrameLabel -> {"Time (d)", "RH(%)"}];

plot2 = ListPlot[topRH, PlotRange -> {All, {40, 100}}, Frame -> True, 
  FrameLabel -> {"Time (d)", "RH(%)"}]

Show[plot1, plot2]

enter image description here

Visually, this result is much better than the original.

$\endgroup$
  • $\begingroup$ Thanks a lot corey979. You are very right about the RHbound being around 68%. In fact I tried 68% and I found a much better fit as well. This means that I should report my measured RHbound as 70+/-2%. You had a very smart answer. $\endgroup$ – amirtech2005 Oct 7 '16 at 0:18
3
$\begingroup$

For what it is worth, even if I let all of the parameters float I can't fit the data with the model supplied by the OP.

The shape of the data looks like a exponential decline superimposed on a line.

RHfittop = 
 FindFit[topRH, a Exp[-b t] + sl t + int,
  {{a, 15}, {b, 0.1}, {sl, -0.03}, {int, 70}}, t]

(* {a -> 15.3197, b -> 0.092728, sl -> -0.0333243, int -> 77.4412} *)

RHfitmid = 
 FindFit[midRH, a Exp[-b t] + sl t + int,
 {{a, 15}, {b, 0.1}, {sl, -0.03}, {int, 70}}, t]

(* {a -> 15.4196, b -> 0.0586449, sl -> -0.0435704, int -> 81.756} *)

RHfitbot = 
 FindFit[botRH, a Exp[-b t] + sl t + int,
 {{a, 15}, {b, 0.1}, {sl, -0.03}, {int, 70}}, t]

(* {a -> 17.7906, b -> 0.0404087, sl -> -0.0406912, int -> 82.0845} *)

Now define three functions based upon the fits

top[t_] := Evaluate[a Exp[-b t] + sl t + int /. RHfittop]

mid[t_] := Evaluate[a Exp[-b t] + sl t + int /. RHfitmid]

bot[t_] := Evaluate[a Exp[-b t] + sl t + int /. RHfitbot]

Plot the data versus the fits

Show[

 Plot[
  {top[t], mid[t], bot[t]},
  {t, 0, 200},
  PlotRange -> {70, 95},
  (*PlotStyle\[Rule]Black,*)
  PlotLegends -> "Expressions",
  Frame -> True,
  FrameLabel -> {"Time (d)", "RH(%)"}
  ],
 ListPlot[
  {topRH, midRH, botRH},
  PlotRange -> All,
  PlotStyle -> Black
  ]
 ]

Mathematica graphics

$\endgroup$
  • $\begingroup$ Thanks Jack LaVigne for your time. I knew the RH function from the beginning but I wanted to find "di" as I mentioned in my question using that expression I have used in Findfit. corey979 nailed it. $\endgroup$ – amirtech2005 Oct 7 '16 at 0:22

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.