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$
1
  • $\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$
    – user484
    Sep 5, 2014 at 15:31

1 Answer 1

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 and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.