Message Boards Message Boards

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

NDSolve inside For loop, reseting conditions after each zero

Anup Singh
Posted 10 years ago

Hi guys, basically what I am trying to do is see whether my graphical solution to an ODE will give a zero at the correct point. For a certain set of initial conditions, this should be the case up until some point where the zeros will vanish (or occur after that specified point). I tried to do this on code but it is not working how I thought - I want to have it see whether x[t-spec] == 0 and if it does, then I will add a certain value (0.000015) to "r" (from the code) and check the next set of initial conditions (and if they satisfy the x[t-spec] == 0 as well, then add a 0.000015 to r). I am interested in finding the "r" value at which beyond it, there will be no more zeros at that specified point. Here is my code below:

    r = 0 (*defines the critical radius variable*)
For[i = 0.000015, 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}]
    Plot[Evaluate[x[t] /. s], {t, 0, .1}]    (*plots sol'n to diff eq*)
    If[{x[.006] /. s} == 0, r = r + 0.000015, r = r] (*checks if zero at specified time - adds .000015 to r if true)
  }
Print[r]

The logic makes sense but it seems like I am having trouble with the syntax. I am most concerned with the last two lines in the loop - where I plot and evaluate to see whether there is a zero at the specified time. Any help with this would be needed! Also, whenever I run it, I keep getting it showing up as r = 0, so I'm guessing it's not adding to "r" - that's the only clue I can see so far. Thanks for the help! Please reply if you have any questions about this so I can clear it up if you're unsure on something I said!
POSTED BY: Anup Singh
3 Replies
Sort By:
Bill Simpson
Bill Simpson
Posted 10 years ago
POSTED BY: Bill Simpson
Anup Singh
Anup Singh
Posted 10 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
Bill Simpson
Bill Simpson
Posted 10 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.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract