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