2
$\begingroup$

I have mathematical model based on equation:

$y(x) = 1 - \sum_{b=1}^{\inf} 2c_b \sin(c_b x/L) \exp(ax) $ (1),

where $c_b$ are the positive roots of equation:

$c_b \cot c_b + aL = 0$ (2).

$L$ is a known constant, and $a$ is parameter. I would like to fit the model to experimental data and obtain value of $a$.

How to implement model based on system of such equations to Mathematica's FindFit? In particular, my problem lies in the second equation, which has more than one solution, and every positive one must be accounted in eq.1. Any help appreciated.

$\endgroup$
  • $\begingroup$ Wait, equation (2) has arbitrarily large roots for $c_b$, so your model contains massively oscillating terms $c_b\sin(c_b x/L)$? That doesn't sound right. $\endgroup$ – Rahul Sep 5 '14 at 15:31
1
$\begingroup$

Using only the first cbMax amount of $c_b$s, we get:

getmodel[a_, cbMax_] := With[{cb = Table[FindRoot[cb Cot[cb] + a L, {cb, Pi (n + 1/2), 
n Pi, (n + 1) Pi}][[1, 2]], {n, 0, cbMax}]},
1 - Sum[2 cb[[i]] Sin[cb[[i]] x/L] E^(a x), {i, cbMax + 1}]]

And now notice that for the rest of the $c_b$s, the difference becomes constant (Pi) very fast:

ListPlot[Differences[Table[FindRoot[cb Cot[cb] + 1, {cb, Pi (n + 1/2),
n Pi, (n + 1) Pi}][[1, 2]], {n, 0, 100}]]]

enter image description here

Hence a good approximation to account for the rest of the $c_b$s could be

FullSimplify[Sum[2 (cb[[-1]] + i Pi) Sin[(cb[[-1]] + i Pi) x/L] E^(a x), {i, \[Infinity]}]]
(*(2 E^((a + (I Pi)/L) x) ((Pi + cb[[-1]]) Sin[(x cb[[-1]])/L] - 
cb[[-1]] Sin[(x (Pi + cb[[-1]]))/L]))/(-1 + E^((I Pi x)/L))^2*)

So the final expression for the model function becomes:

getmodel[a_, cbMax_] := With[{cb = Table[FindRoot[cb Cot[cb] + a L,
{cb, Pi (n + 1/2), n Pi, (n + 1) Pi}][[1, 2]], {n, 0, cbMax}]},
1 - Sum[2 cb[[i]] Sin[cb[[i]] x/L] E^(a x), {i, cbMax + 1}] -
((2 E^((a + (I Pi)/L) x) ((Pi + cb[[-1]]) Sin[(x cb[[-1]])/L] - 
cb[[-1]] Sin[(x (Pi + cb[[-1]]))/L]))/(-1 + E^((I Pi x)/L))^2)]

Depending on the x-values of your data, a and L choose some cbMax that is sufficiently high for the model to converge and then do:

error[a_?NumericQ] := (
y[x_] = getmodel[a, cbMax]; Total[(dataY - Map[y, dataX])^2])

FindArgMin[error[a], {a, 1, 2}]

And there is a.

$\endgroup$

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.