NDSolve problem with switching equations

[Haibo Min]

Hello, everyone.

I am trying to solve a switching ND problem using NDSolve. Specifically,
suppose the system equation is

y1'[t]==y2;
y2'[t]==1000(1-y1^2)y2-y1;

y1(0)==0;y2(0)==1;

Then every half a second (0.5s), the system equation transforms to

y1'[t]==y2;
y2'[t]==1000(1+y1^2)y2;

and then, after 0.5s, it transforms back again.

How to address this problem?

Thank you!

haibo

[DrMajorBob]

First, we need a square wave of amplitude one, zero up to 0.5, 0 from .05
to 1.0, with period 1.

The following will do well enough:

Plot[UnitStep@Sin[2 Pi t], {t, 0, 2}]

So here's a solution with NDSolve:

Clear[y1, y2]
eqns = {y1'[t] == y2[t],
y2'[t] == 1000 (1 - y1[t]^2) (y2[t] - y1[t] UnitStep@Sin[2 Pi t]),
y1[0] == 0, y2[0] == 1};
errors = Take[eqns, 2] /. Equal -> Subtract;
{y1, y2} = {y1, y2} /. First@
NDSolve[eqns, {y1, y2}, {t, 0, 2}];
Plot[{y1@t, y2@t}, {t, 0, 1}, AxesOrigin -> {0, 2},
PlotStyle -> {Red, Blue}]

Accuracy is good for most of the graph:

Plot[Norm@errors, {t, .015, 2}, PlotRange -> All]

But not all:

Plot[Norm@errors, {t, 0, .015}, PlotRange -> All]

(Very interesting plots!)

Increasing WorkingPrecision works well:

Clear[y1, y2]
eqns = {y1'[t] == y2[t],
y2'[t] == 1000 (1 - y1[t]^2) (y2[t] - y1[t] UnitStep@Sin[2 Pi t]),
y1[0] == 0, y2[0] == 1};
errors = Take[eqns, 2] /. Equal -> Subtract;
{y1, y2} = {y1, y2} /. First@
NDSolve[eqns, {y1, y2}, {t, 0, 1}, WorkingPrecision -> 40];
Plot[{y1@t, y2@t}, {t, 0, 1}, AxesOrigin -> {0, 1},
PlotStyle -> {Red, Blue}]

Do the error plots again, and you'll find they're very slow.

ecause the interpolations involve a lot of points? Because Plot is working
hard to identify features?

I'm not sure.

Bobby

--
DrMaj...@yahoo.com

[Haibo Min]

Thank you very much, Bobby. I think it is really a brilliant solution. As
for the slow plotting in the second case, I think it is due to your second
explanation (Plot is working hard to identify features).

Best regards,

Haibo

[DrMajorBob]

Here are simpler choices for the box wave.

box[x_] = Piecewise@{{1, FractionalPart@x < .5}}
Plot[box@t, {t, 0, 2}]

or

box = If[FractionalPart@# < .5, 1, 0] &;

But, strangely enough, they don't speed things up in NDSolve.

Clear[y1, y2]
eqns = {y1'[t] == y2[t],
y2'[t] == 1000 (1 - y1[t]^2) (y2[t] - y1[t] UnitStep@Sin[2 Pi t]),
y1[0] == 0, y2[0] == 1};
errors = Take[eqns, 2] /. Equal -> Subtract;

First@Timing[{y1, y2} = {y1, y2} /. First@

NDSolve[eqns, {y1, y2}, {t, 0, 2}]]
Plot[{y1@t, y2@t}, {t, 0, 1}, AxesOrigin -> {0, 2},
PlotStyle -> {Red, Blue}]

0.021479

Clear[y1, y2]
box = If[FractionalPart@# < .5, 1, 0] &;

eqns = {y1'[t] == y2[t],

y2'[t] == 1000 (1 - y1[t]^2) (y2[t] - y1[t] box@t), y1[0] == 0,

y2[0] == 1};
errors = Take[eqns, 2] /. Equal -> Subtract;

First@Timing[{y1, y2} = {y1, y2} /. First@

NDSolve[eqns, {y1, y2}, {t, 0, 2}]]
Plot[{y1@t, y2@t}, {t, 0, 1}, AxesOrigin -> {0, 2},
PlotStyle -> {Red, Blue}]

0.020186

Bobby

--
DrMaj...@yahoo.com

[Bill Rowe]

On 1/4/10 at 5:59 AM, btr...@austin.rr.com (DrMajorBob) wrote:

>First, we need a square wave of amplitude one, zero up to 0.5, 0
>from .05 to 1.0, with period 1.

>The following will do well enough:

>Plot[UnitStep@Sin[2 Pi t], {t, 0, 2}]

=46WIW, if you are using vers 7 or later there is the built in
function SquareWave. That is

Plot[.5 + .5 SquareWave[t],{t, 0, 2}]

produces the same graphic as your code. And using .5 + .5
SquareWave[t] in place of UnitStep@Sin[2 Pi t] in the arguments
to NDSolve appears to get the identical solution. At least the
error plots look the same. On my machine, subjectively, there
didn't seem to be any performance difference for either function
for producing a square wave.

[Albert Retey]

Hi,

> Here are simpler choices for the box wave.
>
> box[x_] = Piecewise@{{1, FractionalPart@x < .5}}
> Plot[box@t, {t, 0, 2}]
>
> or
>
> box = If[FractionalPart@# < .5, 1, 0] &;
>
> But, strangely enough, they don't speed things up in NDSolve.
>

I think there is no reason to expect they would. I think the problem is
that a numeric differential equation solver will always use a lot of
sampling points near one of the steps. Actually the resolution it will
use is just limited by the precision it tries to achieve, since it would
need infinitesimaly small steps to resolve the step. So if you want to
see speedup, one pragmatic approach would be to use a "smoothing" of
your step function, e.g. by replacing it with an appropriate ActTan.
Then you can control the resolution of the step function independently
of the desired precision. In most pratical problems you will have an
idea about what a reasonable resolution for such a step function would be...

hth,

albert