Is there any in-built function or a recommended package that enables one to find, say, 10 Parameters, which are used as boundary conditions in NDSolve, that minimize $\chi^2$?



Mexp = {171.7, 6.19*10^-1, 1.27*10^-3, 2.89, 5.5*10^-2, 2.9*10^-3, 
   1.74624, 1.0272*10^-1, 4.8657*10^-4};
Merrors = {0.05*Mexp[[1]], 0.05*Mexp[[2]], 0.05*Mexp[[3]], 
   0.05*Mexp[[4]], 0.05*Mexp[[5]], 0.05*Mexp[[6]], 0.05*Mexp[[7]], 
   0.05*Mexp[[8]], 0.05*Mexp[[9]]};

chisquare = Sum[((Mexp[[i]] - m[[i]])/Merrors[[i]])^2, {i, 9}];

Any ideas would be much appreciated!


ParametricNDSolve will return a numerical ODE solution with any number of free parameters. This parametric solution can then be fed into NonLinearModelFit (or whatever home-brew chi-squared algorithm you want to cook up) to find the best-fit values for the parameters.

As an example, suppose we want have the ODE $y''(x) = - y(x)$, with initial conditions $y(0) = 0$ and $y'(0) = a$. We want to find the value of $a$ for which the solution best fits the points $(1,1)$, $(2, 1.5)$, and $(3, 0)$. This is accomplished by the following code:

data = {{1, 1}, {2, 1.5}, {3, 0}};
soln = ParametricNDSolve[{y''[x] == -y[x], y[0] == 0, y'[0] == a}, y, {x, 0, 5}, {a}];
model[a_, x_] := (y /. soln)[a][x]
fit = NonlinearModelFit[data, model[a, x], {a}, x]

(* FittedModel[InterpolatingFunction[Domain: {{0.,5.}} Output: scalar][x]] *)

Show[Plot[fit[x], {x, 0, 5}], ListPlot[data]]

enter image description here

Of course, in this case, we know what the solutions will be ($y(x) = a \sin x$), and they end up yielding a linear model in $a$; so you could probably simplify the code or use different algorithms for this "toy case". But these nice features will not be the case for a general ODE, and so I wrote up this more general method.


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.