0
$\begingroup$

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.

$\endgroup$
4
  • $\begingroup$ What error message do you get? Do you still get an answer? $\endgroup$
    – MarcoB
    Dec 3, 2023 at 2:27
  • 1
    $\begingroup$ The evaluation of DATA doesn't work. NDSolve throws copious errors at that point. Is that not the case for you? $\endgroup$
    – MarcoB
    Dec 3, 2023 at 2:33
  • $\begingroup$ @MarcoB that's my case... But everything appears when using the FindMinimum. Errors related to NDSolve and ReplaceAll... If you evaluate chi2 function, everything is ok $\endgroup$
    – Syn1110
    Dec 3, 2023 at 3:47
  • 1
    $\begingroup$ 1. So something is already wrong in the line DATA=..., please correct the description. 2. Use something rather than t in the Table will fix the issue. (Don't forget t is already used in those NDSolves and you haven't localize them! ) 3. {w, 0.9} in last line is obviously wrong. 4. You need _?NumericQ, see this post for more info: mathematica.stackexchange.com/a/26037/1871 $\endgroup$
    – xzczd
    Dec 3, 2023 at 4:28

2 Answers 2

2
$\begingroup$

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]
$\endgroup$
3
  • $\begingroup$ Wow! Thanks... The use of ParametricNDSolveValue and NonlinearModelFit solved those issues. What I was wondering is that if could plot the contour plots for sigma amplitudes using the results of NonlinearModelFit and the definition of chi-square function used above. $\endgroup$
    – Syn1110
    Dec 3, 2023 at 23:50
  • $\begingroup$ @Syn1110 Graphics[result["ParameterConfidenceRegion"], Frame -> True, FrameLabel -> {a, b}] shows the confidence ellipse of the parameters. Please have a look at result["Properties"] to see all the goodies you can get out of the fit result. $\endgroup$
    – Roman
    Dec 4, 2023 at 7:43
  • $\begingroup$ Thanks a lot! That only shows confidence region at 1sigma amplitude right? Is there a way to obtain 2sigma and 3sigma? $\endgroup$
    – Syn1110
    Dec 5, 2023 at 0:30
0
$\begingroup$

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:

enter image description here

I think this is a bug.

$\endgroup$
4
  • $\begingroup$ No, it's not. NDSolve simply doesn't have attribute like HoldAll that adjusts the evaluation order, this is by design. $\endgroup$
    – xzczd
    Dec 3, 2023 at 9:26
  • $\begingroup$ @xzczd Do you see any reason why this should be by design? $\endgroup$ Dec 3, 2023 at 9:30
  • $\begingroup$ Otherwise things like this will happen on Solve/DSolve/NDSolve: opt=PlotStyle->Red; Plot[x,{x,0,1},opt] $\endgroup$
    – xzczd
    Dec 3, 2023 at 9:51
  • $\begingroup$ Actually in early versions, even domain = {x, 0, 1}; Plot[x, domain] is enough to cause trouble: i.stack.imgur.com/LpOl8.png This is all because of the HoldAll. If Solve, etc. behaves similarly, it'll be rather inconvenient. (Just imagine the world that code like eq = {a + b == 0, a - 2 b == 3}; vars = {a, b}; Solve[eq, vars] results in failure! ) $\endgroup$
    – xzczd
    Dec 3, 2023 at 11:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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