Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2007
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2007

[Date Index] [Thread Index] [Author Index]

Search the Archive

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
  • Next by thread: Re: conditional plotstyles in ListPlot [additional note]