Skip to main content
Mathematica

Return to Answer

added 815 characters in body
Source Link
Roman
Roman
  • 49.9k
  • 2
  • 57
  • 131

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]

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]
Source Link
Roman
Roman
  • 49.9k
  • 2
  • 57
  • 131

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} $$