Prev: Difficulty with NDSolve (and DSolve)
Next: Wrong ODE solution in Mathematica 7?
From: Haibo Min on 3 Jan 2010 03:39
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

From: DrMajorBob on 4 Jan 2010 06:00
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

From: Haibo Min on 4 Jan 2010 06:01
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
>


From: Bill Rowe on 5 Jan 2010 01:47
On 1/4/10 at 5:59 AM, btreat1(a)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(a)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(a)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.


 | 
Pages: 1
Prev: Difficulty with NDSolve (and DSolve)
Next: Wrong ODE solution in Mathematica 7?