Re: numerical_solution

• To: mathgroup at smc.vnet.net
• Subject: [mg73367] Re: numerical_solution
• From: Paul Abbott <paul at physics.uwa.edu.au>
• Date: Wed, 14 Feb 2007 05:19:21 -0500 (EST)
• Organization: The University of Western Australia
• References: <eqsa2a\$jmb\$1@smc.vnet.net>

```In article <eqsa2a\$jmb\$1 at smc.vnet.net>,
"j. r. campanha" <campanha at rc.unesp.br> wrote:

> How can I know which is the best numerical solution?

Easy -- substitute the solution back into the differential equation!

> w == 6;
> A == 0.1;
> B == 0.17;
> F == 0.501;

There were a number of spurious "=" characters in your post. And it is
preferable to use lower case letters for variable names:

w = 6;
a = 0.1;
b = 0.17;
f = 0.501;

> solution1 = NDSolve[{x'[t] == v[t], v'[t] == -(B)*(v[t]) - (A)*(((x[t]
> ))/(1 - (x[t])^2)) + (F)*(Cos[w*t]),x[0] == 0.4, v[0] == 0.3}, {x, v}, =
> {t, 260, 270}];
>
> ParametricPlot[Evaluate[{x[t], v[t]} /. solution1], {t, 260, 270}, Frame ->
>  True]

Solve the differential equations, here for 0 < t < 270,

solution1 = NDSolve[{x'[t] == v[t], v'[t] == -b v[t] -
a x[t]/(1 - x[t]^2) + f Cos[w t], x[0] == 0.4, v[0] == 0.3},
{x, v}, {t, 0, 270}];

and substitute the solution back into (the difference of the right and
left-hand sides of) each differential equation:

ParametricPlot[Evaluate[{x'[t] - v[t],
v'[t] + b v[t] + a x[t]/(1 - x[t]^2) - f Cos[w t]} /. solution1],
{t, 0, 270}, Frame -> True]

As you will see, if you plot over the full range of t, the solution is
approaching a simple limit cycle.

If you try this with Method -> StiffnessSwitching, the numerical
solution obtained is not very good:

solution2 = NDSolve[{x'[t] == v[t], v'[t] == -b v[t] -
a x[t]/(1 - x[t]^2) + f Cos[w t], x[0] == 0.4, v[0] == 0.3},
{x, v}, {t, 0, 270}, Method -> StiffnessSwitching];

ParametricPlot[Evaluate[{x'[t] - v[t],
v'[t] + b v[t] + a x[t]/(1 - x[t]^2) - f Cos[w t]} /. solution2],
{t, 0, 270}, Frame -> True]

Finally, note that there is no need to convert second order differential
equations into systems of first order equations:

solution3 = NDSolve[{x''[t] + b x'[t] + a x[t]/(1 - x[t]^2) ==
f Cos[w t], x[0] == 0.4, x'[0] == 0.3}, x, {t, 0, 270}];

Plot[Evaluate[x[t] /. solution3], {t, 0, 270}]

ParametricPlot[Evaluate[{x[t], x'[t]}/. solution3],
{t, 0, 270}, Frame -> True, PlotRange -> All, PlotPoints -> 300]

Plot[Evaluate[x''[t] + b v[t] + a x[t]/(1 - x[t]^2) -
f Cos[w t] /. solution1], {t, 0, 270}]

Cheers,
Paul

_______________________________________________________________________
Paul Abbott                                      Phone:  61 8 6488 2734
School of Physics, M013                            Fax: +61 8 6488 1014
The University of Western Australia         (CRICOS Provider No 00126G)
AUSTRALIA                               http://physics.uwa.edu.au/~paul

```

• Prev by Date: Re: RandomArray from user defined distribution?
• Next by Date: Re: Dotted Lines at Discontinuities
• Previous by thread: Re: numerical_solution