Prev: mathproject10
Next: JLink error
From: DrMajorBob on 5 Jan 2010 01:46
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(a)Sin[2 Pi t]),
y1[0] == 0, y2[0] == 1};
errors = Take[eqns, 2] /. Equal -> Subtract;
First(a)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(a)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

On Mon, 04 Jan 2010 05:00:04 -0600, Haibo Min <haibo.min(a)gmail.com> wrote:

> 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
>
>
>
>
> On Mon, Jan 4, 2010 at 9:08 AM, DrMajorBob <btreat1(a)austin.rr.com> 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(a)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(a)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(a)errors, {t, .015, 2}, PlotRange -> All]
>>
>> But not all:
>>
>> Plot[Norm(a)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(a)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
>>
>>
>> On Sun, 03 Jan 2010 02:42:43 -0600, Haibo Min <haibo.min(a)gmail.com>
>> wrote:
>>
>> 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(a)yahoo.com
>>
>
>


--
DrMajorBob(a)yahoo.com

From: Albert Retey on 6 Jan 2010 05:59
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

 | 
Pages: 1
Prev: mathproject10
Next: JLink error