Re: How can I "perturbate" a NDSolve
- To: mathgroup at smc.vnet.net
- Subject: [mg102129] Re: How can I "perturbate" a NDSolve
- From: Albert Retey <awnl at gmx-topmail.de>
- Date: Thu, 30 Jul 2009 05:32:18 -0400 (EDT)
- References: <200907260754.DAA18931@smc.vnet.net> <h4p3ao$iik$1@smc.vnet.net>
Hi,
> BTW, the perturbation (at that probability) makes virtually no difference
> because, for instance, after running the code,
>
> DownValues[a][[All, -1]]
> Length@%
>
> {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
>
> 180
>
> That is, a was called 180 times, and all the results were zero.
>
> When I change the probability to 0.1, there are lots more a values and
> about 10% are non-zero:
>
> DownValues[a][[All, -1]];
> Length@%
> Length@DeleteCases[%%, 0]
>
> 2100
>
> 208
>
> ...but there's no visible difference in the plot! (Interpolations are
> being used, and that smoothes everything.)
>
> There's a critical point around p = 0.17 or 0.18, where NDSolve begins to
> give error messages and the plot changes radically.
I have also played with the problem and I think that the reason that one
doesn't see a difference is the variable step solver: because of the
step in a[t] when it changes the solver decreases the step size so much
that you don't see any effect. The next step a[t] will be back to
a[t]=0, so if the step size is small enough, the effect is negligible.
If one uses higher probabilities for jumps in a the solver will at some
point find that it can't make the step sizes small enough to make the
perturbation negligible and thus will give error messages and stop
because it can't satisfy its precision or stability checks. I think to
achieve what the original poster wants, one needs to switch to a fixed
step solver method, e.g. like this:
ClearAll[a];
a[t_?NumericQ] := a[t] = (
If[RandomReal[] < 0.001, RandomReal[], 0]
);
sol = NDSolve[{x1'[t] == (Sin[t] + a[t])*x2[t],
x2'[t] == (Cos[t] - a[t]*x1[t]^2)*x1[t], x1[0] == 1,
x2[0] == 0}, {x1[t], x2[t]}, {t, 0, 10},
MaxSteps -> 1000000,
Method -> {"FixedStep", Method -> "ExplicitEuler"}
];
Plot[Evaluate[Flatten[{x1[t], x2[t]} /. sol]], {t, 0, 10},
ImageSize -> 800]
probably there are better submethods, but the results looks as what I
believe the OP could have had in mind...
Another idea would be to determine an a[t] in advance, since it doesn't
have a dependence on the differential equations that should be possible.
But then one would need to ensure that the solver uses steps so it
doesn't miss the intervals where a!=0. If one does so, it is also clear
that one needs to define a rate at which a is supposed to update, which
is the solvers step size in the above cases. If that step size is
variable, the update frequency of a is also variable and that seems to
not make too much sense for me...
hth,
albert
- References:
- How can I "perturbate" a NDSolve expresion?
- From: Iván Lazaro <gaminster@gmail.com>
- How can I "perturbate" a NDSolve expresion?