Chi-square minimization error

Question

I'm trying to fit some data into a model using the standard procedure of minimizing the chi-square function using this code:

ClearAll["Global`*"]
Needs["ErrorBarPlots`"];
a0 = 1;
b0 = 2;

sol[a_, b_] := 
  NDSolve[{x'[t] == -a*x[t] + Sin[t], y'[t] == a*x[t] - b*y[t]^2, 
    z'[t] == b*y[t] - Cos[t], x[0] == 1, y[0] == 0, z[0] == 0}, {x, y,
     z}, {t, 0, 10}, MaxSteps -> Infinity];

F[t_] := x[t] + y[t]^2 + z[t]^3;

wsol[a_, 
   b_] := (wsol[a, b] = 
    NDSolve[{w'[t] - w[t] - ((F[t] /. sol[a, b])[[1]]) == 0, 
      w[0] == 0}, w, {t, 0, 10}]);

model[t_, a_, b_] := (w[t] /. wsol[a, b])[[1]] // Chop;

SeedRandom[1264645]; 
rangoT = Range[0, 9, 0.5];
ndat = Length[rangoT];

DATA = Table[{t, 
    model[t, a0, b0] + RandomVariate[NormalDistribution[0, .1]], 
    RandomReal[{.05, .1}]}, {t, rangoT}];

chi2[a_, b_] := 
 Sum[((DATA[[i, 2]] - model[DATA[[i, 1]], a, b])/DATA[[i, 3]])^2, {i, 
   1, ndat}]

FindMinimum[chi2[a, b], {w, 0.9}, {b, 1.1}]

But when I run it I get error messages from the FindMinimum command related to NDSolve and ReplaceAll. Everything runs just fine until using that command. I was wondering why's not working and what should be changed.

Answer 1

Three comments:

1. solve for all variables simultaneously instead of solving for $(x,y,z)$ first and $w$ second,
2. use NonlinearModelFit instead of writing your own least-squares solver, and
3. going all the way up to $t_{\text{max}}=10$ seems a bit much and the differential equation seems to diverge.

tmax = 6;
sol = ParametricNDSolveValue[{
        x'[t] == -a*x[t] + Sin[t],
        y'[t] == a*x[t] - b*y[t]^2,
        z'[t] == b*y[t] - Cos[t],
        w'[t] - w[t] == x[t] + y[t]^2 + z[t]^3,
        x[0] == 1, y[0] == 0, z[0] == 0, w[0] == 0},
        {x, y, z, w}, {t, 0, tmax}, {a, b}, 
        MaxSteps -> Infinity];
model[t_?NumericQ, a_?NumericQ, b_?NumericQ] := sol[a, b][[4]][t]

a0 = 1;
b0 = 2;
SeedRandom[1264645];
rangoT = Range[0, tmax, 0.5];
ndat = Length[rangoT];

DATA = Table[{t,
              model[t, a0, b0] + RandomVariate[NormalDistribution[0, .1]],
              RandomReal[{.05, .1}]}, {t, rangoT}];

result = NonlinearModelFit[DATA[[All, 1 ;; 2]], model[t, a, b],
           {{a, a0}, {b, b0}}, t,
           VarianceEstimatorFunction -> (1 &),
           Weights -> 1/DATA[[All, 3]]^2]

result["ParameterTable"]

$$ \begin{array}{l|llll} \text{} & \text{Estimate} & \text{Standard Error} & \text{t-Statistic} & \text{P-Value} \\ \hline a & 1.0013 & 0.00103671 & 965.843 & \text{1.8409735010100922$\grave{ }$*${}^{\wedge}$-28} \\ b & 1.9994 & 0.000490316 & 4077.78 & \text{2.423162678987459$\grave{ }$*${}^{\wedge}$-35} \\ \end{array} $$

update: analytic $x(t)$

Note that the differential equation for $x(t)$ can be solved analytically, which improves the stability and accuracy of the result:

DSolve[{x'[t] == -a*x[t] + Sin[t], x[0] == 1}, x[t], t] // FullSimplify
(*    {{x[t] -> ((2 + a^2) E^(-a t) - Cos[t] + a Sin[t])/(1 + a^2)}}    *)

Therefore we can define the problem more accurately as

x[a_, t_] = ((2 + a^2) E^(-a t) - Cos[t] + a Sin[t])/(1 + a^2);
sol = ParametricNDSolveValue[{
        y'[t] == a*x[a, t] - b*y[t]^2,
        z'[t] == b*y[t] - Cos[t],
        w'[t] - w[t] == x[a, t] + y[t]^2 + z[t]^3,
        y[0] == 0, z[0] == 0, w[0] == 0},
        {y, z, w}, {t, 0, tmax}, {a, b}, MaxSteps -> Infinity];
model[t_?NumericQ, a_?NumericQ, b_?NumericQ] := sol[a, b][[3]][t]

Answer 2

The problem is caused by "NDSolve". It should localize the dummy variable "t", but fails to do so. Consider:

t = 0; 
NDSolve[{f'[t] == f[t], f[0] == 1}, f, {t, 0, 10}]

This gives the error message:

I think this is a bug.