Group Abstract Group Abstract

Message Boards Message Boards

0
|
7.2K Views
|
3 Replies
|
0 Total Likes
View groups...
Share
Share this post:

NDSolve inside For loop, reseting conditions after each zero

Posted 12 years ago
POSTED BY: Anup Singh
3 Replies
Posted 12 years ago
POSTED BY: Bill Simpson
Posted 12 years ago

Hey Bill, I am having some trouble still. Sorry for the delay.

r = 0;(*defines the critical radius variable*)
For[i = 1.49999, i <= 1.5, 
   i += 0.000015](*iterates through initial condition values*)s = 
 NDSolve[{x'[
      t] - (2.5*10^-25) ((1.2*10^31 ((9*10^-6) - (1.5*10^-3) t) \
(5*10^5 ((9*10^-6) - (1.5*10^-3) t) - 1000000 x[t] + 
            3) (5*10^5 ((9*10^-6) - (1.5*10^-3) t) + 1000000 x[t] - 
             3)/(1000000000000 ((9*10^-6) - (1.5*10^-3) t)^2 + 
              1000000000000 x[t]^2 - 6000000 x[t] + 9)^(7/
              2)) - ((3 ((9*10^-6) - (1.5*10^-3) t)^3 - 
            12 ((9*10^-6) - (1.5*10^-3) t) x[t]^2))/(x[
              t]^2 + ((9*10^-6) - (1.5*10^-3) t)^2)^(7/2)) == 0, 
   x[0] == i*10^-6}, 
  x, {t, 0, 
   1}]; (*defines the differential equation to be solved with init. \
conditions*)
Print["x[.006]=", x[t] /. s /. t -> .006];
If[x[t] /. s /. t -> .006 <= 10^-10, r = r + 0.000015, 
 Print[Plot[
   Evaluate[x[t] /. s], {t, 
    0, .1}]]];(*checks if zero at specified time-adds .000015 to r if \
true*)
Print[r]

The goal of this code is to be able to solve this differential equations for the given set of initial conditions, and produce an ‘r’ value known as the critical radius. Basically, a particle is moving according to this differential equation, and if it reaches within a certain distance (in this case, 10-10 m) of the trap, it will be captured. The differential equation dictates the position of this particle. The range of positions are from 0 to 1.5 units, and what I want to do is find out what the critical radius is – that is, what is the distance ‘r’ which beyond this distance, the particle will not get trapped. Basically r is set to zero and then add a small increment each time you check a new initial condition – if the particle is captured (d-particle < 10^-10 units away form trap), then you increment farther out. This will continue until this condition is not true, and the graph of this specific instance will be shown. The problem is, it is not doing what I want it to, and I think there is something wrong with my syntax, not logic – I am used to other programming languages so I went about it like those, but I realize Wolfram is different so I’m not sure where I’m going wrong. Thanks for the help and let me know if you have any questions!

POSTED BY: Anup Singh
Posted 12 years ago

From your example code above, I'm guessing you have had some prior experience in another programming language and are trying to use some of that to write Mathematica.

First problem

For[n = 1, n = 3, n++]

does almost nothing and is done.

This does a little more

For[n = 1, n <= 3, n++]

but not much more.

This does a little more

For[n = 1, n <= 3, n++,
 Print[n];
 Print[n^2]
 ]

Note that Mathematica does not use { and } to group multiple statements, it uses semicolons to accomplish that. I am guessing your prior programming experience with other "more normal" programming languages is leading you to expect to group statements with { and }.

r = 0;(*defines the critical radius variable*)
For[i = 0.000015, i <= 0.00015, i += 0.000015,(*iterates through initial condition values*)
 s = NDSolve[{x'[t] - (2.5*10^-25) ((1.2*10^31 ((9*10^-6) - (1.5*10^-3) t) (5*10^5 ((9*10^-6) - (1.5*10^-3) t) -
    1000000 x[t] + 3)(5*10^5 ((9*10^-6) - (1.5*10^-3) t) + 1000000 x[t] - 3)/(1000000000000 ((9*10^-6) -
    (1.5*10^-3) t)^2 + 1000000000000 x[t]^2 -6000000 x[t] + 9)^(7/2)) - ((3 ((9*10^-6) - (1.5*10^-3) t)^3 -
    12 ((9*10^-6) - (1.5*10^-3) t) x[t]^2))/(x[t]^2 + ((9*10^-6) - (1.5*10^-3) t)^2)^(7/2)) == 0, 
     x[0] == i*10^-6}, x, {t, 0, 1}][[1]];
 Print[Plot[Evaluate[x[t] /. s], {t, 0, .1}]];
 Print["x[.006]=", x[t] /. s /. t -> .006];
 If[x[t] /. s /. t -> .006 == 0, r = r + 0.000015](*checks if zero at specified time-adds .000015 to r if true*)
 ];
Print[r]

I'm not sure you really mean to do 100000 plots so I reduced 1.5 to 0.00015 for a first try. There are several other changes in that and you should probably compare that with your original, character by character, word by word, and try to figure out why each change was made. But at least it appears to be correctly iterating and displaying plots for you. If I have misunderstood then please let me know and I will try to fix what I can. Please test this very carefully to make certain that the result is correct before depending on it.

I am concerned that your If statement may never find x[t] to be exactly 0 and thus never True. That may not be good programming.

POSTED BY: Bill Simpson
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard