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.

Improve this question
$\endgroup$
4
  • $\begingroup$ What error message do you get? Do you still get an answer? $\endgroup$ Commented 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$ Commented 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$ Commented 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$ Commented Dec 3, 2023 at 4:28

2 Answers 2

Reset to default
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]
Improve this answer
$\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$ Commented 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$ Commented 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$ Commented 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.

Improve this answer
$\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$ Commented Dec 3, 2023 at 9:26
  • $\begingroup$ @xzczd Do you see any reason why this should be by design? $\endgroup$ Commented 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$ Commented 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.sstatic.net/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$ Commented Dec 3, 2023 at 11:03

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.