# fitting a non-linear curve

I have an equation: $\ddot{x}+(\delta+\epsilon\cos{t})x=0$ known as the Mathieu equation.The $\delta-\epsilon$ parameter space of this equation looks something like

The red lines in this diagram indicate that if a point on the red lines is chosen, and the corresponding value of $(\delta,\epsilon)$ is plugged into the Mathieu equation, then it is guaranteed that one would get a periodic solution. I have managed to generate these curves by something known as the arc-length continuation method, thus I have the exact data to plot these curves with me.

On the other hand, I have an analytic relation approximating these red lines which looks something like: $\delta=\frac{1}{4}+f_1(k)\epsilon+f_2(k)\epsilon^2+f_3(k)\epsilon^3+...$

Here, the functions $f_1(k),f_2(k)...$ are known and are highly nonlinear. This equation is for just one of these red lines. This way I have equations for all the red lines. Also, it is only possible to approximate these curves till an upper bound on $\epsilon$, as I have only one degree of freedom (the parameter k) in analytic $\delta-\epsilon$ relation.

First let me choose one particular red line and the corresponding $\delta-\epsilon$ relation for it.Now, what I want to do is two things: 1) I want to find the optimum $\epsilon$ by which I can fit $\delta-\epsilon$ curve. 2) I want to find the optimum $k$ for that optimum $\epsilon$.

I want to find both the above things by some numerical algorithms which converge. This is necessary as my functions $f_1(k),f_2(k)...$ are highly non-linear in nature.

I will tell you where I have reached so far. I have no idea how to even approach point 1). For point 2), I have managed to implement the mean-square errors approach, by which I tried minimizing the error, but that leaves me to solve a polynomial with powers of k like $k^{3500}$. This polynomial is going to have a lot of local minimas and I have no clue of any numerical algorithms which help me in finding a global minima.

Finally, if I manage to implement 1) and 2), I want to do some analysis (like R-squared coefficient) to see how well I have fit my curve.

Please help me out if you know of any numerical procedures that can achieve the above. Thanks.

• The question is a little too hard to follow. What does it mean to find "the optimum $\epsilon$"? What is $k$? Exactly what is being fitted to what? – Kirill Nov 3 '14 at 2:51

First let me choose one particular red line and the corresponding $\delta$−$\epsilon$ relation for it.Now, what I want to do is two things: 1) I want to find the optimum $\epsilon$ by which I can fit $\delta$−$\epsilon$ curve. 2) I want to find the optimum $k$ for that optimum $\epsilon$.

Regarding point 1): you have to define what you mean by "an optimal $\epsilon$". Usually, this definition suggests an objective function (and potentially a feasible set for that objective function) you can then use as a starting point for a problem formulation.

For point 2), I have managed to implement the mean-square errors approach, by which I tried minimizing the error, but that leaves me to solve a polynomial with powers of $k$ like $k^{3500}$. This polynomial is going to have a lot of local minimas [sic] and I have no clue of any numerical algorithms which help me in finding a global minima [sic].

There are two major approaches to global optimization: deterministic, and nondeterministic. Assuming your feasible set is compact (highly likely), a properly formulated least-squares approach will yield a twice continuously differentiable objective function, so your formulation will be amenable to deterministic methods. The advantage of these methods is that they return a provably good solution, in that the optimal objective function value can be no more than a certain tolerance (user-specified) less than the returned optimal objective function value. The disadvantage is that these methods have a run time that scales exponentially with problem size, so these methods are likely to take a long time if you are fitting more than a handful of parameters. An example of a solver you can try is BARON; if you can somehow pose your optimization using the GAMS or AMPL modeling languages, you can then submit a job on the NEOS server to solve your problem.

Nondeterministic algorithms (e.g., genetic algorithms, particle swarm) do not generate provably good solutions, and merely explore the feasible set of a problem in a random fashion. If you have more than say, 30-300 parameters, you may consider using an algorithm of this type simply to explore your search space in the hopes of calculating a feasible point with a reasonably good objective function value.

You could also try to multistart a local optimization algorithm: use the same approach you've been using, but seed the local optimizer with lots of different initial guesses spread out over your feasible set. You may have to content yourself with selecting a provably good local optimum, even if it might be globally suboptimal.

Finally, if I manage to implement 1) and 2), I want to do some analysis (like R-squared coefficient) to see how well I have fit my curve.

You're likely to get a better answer re: statistical analysis on CrossValidated, a sister Stack Exchange site. As a non-statistician who has taken basic statistics a few times, I can only provide basic advice: at first glance, I would look at a plot of the residual values for your fit to see if there is any systematic bias in the residual (error) values, to see if you are violating any assumptions about normality in the error, which will give you an idea if the model and fit criteria you've selected are appropriate. If you have multiple candidate models and you are looking to select among them, you might also look at model selection statistics, such as the Akaike Information Criterion.

• Basically, the red lines are nothing but a lot of data. Now, my optimum $k$ would depend on how much data I take to formulate an expression for the mean square error in terms of k. So, if I take data till $\epsilon=2$, I will get a different optimum k than when I take data till $\epsilon=4$. Also, after a particular value of $\epsilon$, a good fit would cease to exist (as there wouldn't be a unique value for k anymore). – zynga Nov 3 '14 at 3:38