Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2010

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: FindMinimum and Module

  • To: mathgroup at smc.vnet.net
  • Subject: [mg113746] Re: FindMinimum and Module
  • From: Daniel Lichtblau <danl at wolfram.com>
  • Date: Wed, 10 Nov 2010 06:28:47 -0500 (EST)

Versiane UFOP wrote:
> Dear Mathematica users
>
>
> I would like to request any help with the code below, which I have been
> struggling in the last couple of weeks. I am trying to use FindMinimum
> to find the best k value that fits experimental data (extra=E7=E3o) using
> the code enclosed in a Module structure. This is a heterogeneous
> kinetics problem where the Shrinking Core Model is modified to
> incorporate particle size distribution and k is a rate constant. As many
> of you always express here FindRoot does not produce the expected
> outcome and give the standard message ""The function value...is not a
> real number at {k} == {0.03`}". I would very much appreciate any help.
>
> The code is:
>
>
> mean :== 100;
> std :== 10;
> a == 5.032;
> b == 2.974;
> dmax == 190;
> u == 1;
>
> extra=E7=E3o == {{0, 0.}, {10, 0.068973}, {20, 0.105011}, {30, 0.127376},
> {40, 0.152779}, {50, 0.177607}, {60, 0.186656}, {80, 0.224261}, {100,
> 0.25433}, {120, 0.275795}}
>
> f[D_] :== (1/(b^a*Gamma[a]))*D^(a - 1)*Exp[-D/b]
>
> t :== Table[extra=E7=E3o[[i, 1]], {i, Length[extra=E7=E3o]}]
>
> model == ConstantArray[0, {Length[extra=E7=E3o]}]
>
> conversion[(k_)?NumericQ] :== Module[{dmin, rn, n, d, y, w, unconvt},
>       For[j == 1, j <== Length[t], j++,
>         dmin == (k*t[[j]])^0.5;
>         rn == (dmax - dmin)/(2*u) + 1; n == Round[rn];
>         d == Table[dmin + (2*i - 1)*u, {i, n}];
>     w == Table[y /. FindRoot[1 - 3*(1 - y)^(2/3) + 2*(1 - y) ====
> (k/d[[i]]^2)*t[[j]], {y, 0.3}], {i, n}];
>     unconvt ==     1 - Sum[(1 - w[k][[i]])*N[Integrate[f[D], {D, d[[i]] -
> u, d[[i]] + u}]], {i, 1, n}];
>           model[[j]] == unconvt; ];
>           model1 == {t, model};
>           x == Transpose[model1]; ]
>          fit[(k_)?NumericQ] :==  Sum[(conversion[k][extra=E7=E3o[[i, 1]]] -
> extra=E7=E3o[[i, 2]])^2, {i, 1, Length[extra=E7=E3o]}];
>
> FindMinimum[{fit[k], k > 0}, {k, 0.03}]


This code is a bit of a mess. I took a few guesses as to what it is
being attempted.

mean == 100;
std == 10;
a == 5.032;
b == 2.974;
dmax == 190;
u == 1;

extract == {{0, 0.}, {10, 0.068973}, {20, 0.105011}, {30, 0.127376},
   {40, 0.152779}, {50, 0.177607}, {60, 0.186656},
   {80, 0.224261}, {100, 0.25433}, {120, 0.275795}};

f[D_] :== (1/(b^a*Gamma[a]))*D^(a - 1)*Exp[-D/b]

t == extract[[All, 1]]

conversion[k_?NumericQ] :== Module[
   {dmin, rn, n, d, y, w, unconvt},
       model == Table[
	    dmin == (k*t[[j]])^(1/2);
         rn == (dmax - dmin)/(2*u) + 1;
		n == Round[rn];
         d == Table[dmin + (2*i - 1)*u, {i, n}];
		w == Table[y /.
		  FindRoot[1 - 3*(1 - y)^(2/3) + 2*(1 - y) ==== (k/d[[i]]^2)*t[[j]],
           {y, 0.3}], {i, n}];
         unconvt == 1 - Re[Total[Table[(1 - w[[i]])*
		  NIntegrate[f[D], {D, d[[i]] - u, d[[i]] + u}], {i, 1, n}]]];
         unconvt
		, {j, Length[extract]}];
       model1 == {t, model};
     Transpose[model1]
	]

fit[(k_)?NumericQ] :== With[
   {vec==conversion[k][[All,2]]-extract[[All,2]]}, vec.vec]

FindMinimum[{fit[k], k > 0}, {k, 0.03}]

Out[29]== {0.000716899, {k -> 0.0334664}}

This last will take some time to run.

Daniel Lichtblau
Wolfram Research


  • Prev by Date: Re: pure function
  • Next by Date: Re: command to save as .m file
  • Previous by thread: Re: Find max of a Find root function
  • Next by thread: Repost (Boundary Condition NDSolve)