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