# Fitting composed model to experimental data

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.

• 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. – user484 Sep 5 '14 at 15:31

## 1 Answer

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}]]]


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.