Message Boards Message Boards

How do I use NDSolve within a Do loop?

Posted 11 years ago
I am developing a code which uses NDSolve to find solutions to two sets of two simultaneous differential equations, I would like to use a Do loop or a For loop to change some variables and repeat the calculations, to find the point where the sets diverge, however if I use a Do loop or a For loop, even without changing any other information, for instance, using a Do loop to do the calculations only once, NDSolve no longer works.
  m = 9.1*10^-28; (*Mass of Electron in g*)(*Initialising constants, \
  
  starting with universal constants*)
  
  c = 3*10^10;(*Speed of light in cm/s*)
  
  e = 4.8*10^-10;(*Electron charge in StatColoumb*)
  
  \[Lambda] = 10^-4;(*Wavelength of laser in cm*)
 
 r0 = 30*10^-4;(*Spot size in cm*)
 
 n0 = 6*10^17;(*Electron density in cm^-3*)
 
 energy = 10^5;(*Energy of laser pulse in erg*)
 
 FWHM = 100*10^-15;(*Full-width half maximum of power distribution in \
 
 seconds*)
 
 \[Tau] = FWHM/(2 Sqrt[Log[2]]) // N;(*Time Constant in seconds*)
 
 \[Omega]0 = 2*Pi*c/\[Lambda];
 
 \[Omega]p = Sqrt[4*Pi*(e^2)*n0/m];(*Plasma Frequency in rad.s^-1*)
 
 kp = \[Omega]p/c // N;
 
 \[Tau]p = 2*Pi/\[Omega]p;
 
 vp = c (1 - (\[Omega]p/\[Omega]0)^2)^0.5;
 
 kR = r0*kp;
 
 kZ = kp*vp*\[Tau];
 
 P0 = (energy/FWHM)*Sqrt[2*Log[2]/Pi];(*Peak Power in erg/s*)
 
 I0 = 2 P0/(Pi*r0^2);(*Peak Intensity in erg/(cm^2/s)*)
 
 a0 = Sqrt[(2/(c*Pi))*(\[Lambda]*e/(m*c^2))^2*I0];(*Amplitude constant*)
 
 
 
 asquared1[Z_, R_, T_] :=
 
   a0^2*Exp[-(((Z - T)/kZ)^2)]*Exp[-2 ((R/kR)^2)];
 
 
 
 Tmin = -20;
 
 Tmax = 20;
 
 Zmin = 0;
 
 Zmax = 0.1;
 
 Rmin = -0.1;
 
 Rmax = 0.1;
 
 T0 = -20;
 
 
 
 Do[
 
 
 
  s = NDSolve[{D[
 
          u1[Z, R, T], {T, 2}]*(1 + 1.5*D[u1[Z, R, T], T]^2 +
 
           0.5*D[v1[Z, R, T], T]^2) + u1[Z, R, T] -
 
        D[v1[Z, R, T], T]*D[u1[Z, R, T], T]*v1[Z, R, T] == -0.5*
 
        D[asquared1[Z, R, T], R],
 
      D[v1[Z, R, T], {T, 2}]*(1 + 1.5*D[v1[Z, R, T], T]^2 +
 
           0.5*D[u1[Z, R, T], T]^2) + v1[Z, R, T] -
 
        D[u1[Z, R, T], T]*D[v1[Z, R, T], T]*u1[Z, R, T] == -0.5*
 
        D[asquared1[Z, R, T], Z],
 
      (D[u1[Z, R, T], T] /. T -> T0) == 0, v1[Z, R, T0] == 0,
 
      u1[Z, R, T0] == 0, (D[v1[Z, R, T], T] /. T -> T0) == 0}, {u1,
 
      v1}, {Z, Zmin, Zmax}, {R, Rmin, Rmax}, {T, Tmin, Tmax},
 
     MaxSteps -> Infinity, MaxStepSize -> {0.1, 0.1, 0.1},
 
     Method -> "ExplicitRungeKutta"]
 
    

   First[NMaximize[{v1[0, 0, T] /. s[[1]], 10 < T < 20}, T]]



, {1}]
Removing the Do  loop, the code works, but if I add it or a For loop, it gives the errros NMaximize::nnum: and Thread::tdlen: Could anyone suggest where I'm going wrong? Thank you.

Sorry for any mistakes with formatting, this is my first time using this website.
POSTED BY: Henry Banks
6 Replies
One of the confusing things for me when I started using Mathematica was the use of the semicolon.  Normally you're just typing in a cell and leaving a semicolon off the end of a line means the output of that line will appear in its own output cell.  But your typed code is interpreted differently when it appears inside another command, such as  Do  in the original question,  In that case, you need a semicolon or the lines are multiplied together.  The semicolon in that case stands for  CompoundExpression.

As for  Do, it is like the command  Table  but with no output, which Bill Simpson has pointed out.  Most of the time in my work, I need  Table  and only occasionally  Do.
POSTED BY: Michael Rogers
Posted 11 years ago
Wow, thank you, I already had the functionality, but that code is so much neater than my previous set up and with this I can compare the first set of results before the code has finished evaluating.

I'm currently converting the code into ParametricNDSolve, it's actually much simpler than I thought it would be (or possibly my initial impulse is just rather lazy) but with a working NDSolve, it actually requires almost no changes.

Thank you to everyone, this site is incredibly helpful.
POSTED BY: Henry Banks
Posted 11 years ago
If you are filling a Table with results and you are hesitant about changing from NDSolve right now then perhaps something like this might work
Table[ First[NMaximize[{v1[0, 0, T] /. NDSolve[AllYourNDSolveStuff][[1]], 10 < T < 20}, T]], {j, 1, 10}]
That creates the Table directly, rather than calculating a result using NDSolve and NMaximize and then storing the result in the Table, and hopefully eliminates all the confusion about how to get the result stored in a Table.
.
I assume you want to change some parameter with each iteration of your Do loop, but you didn't show that in your original post, so you will need to replace j with the name you want to iterate and choose appropriate bounds and step sizes.

Then you can deal with figuring out ParametricNDSolve another day.
POSTED BY: Bill Simpson
Posted 11 years ago
Yeah, I added a semicolon and it worked, though I don't entirely understand how, is it something to do with Do not producing output, if you know why this occurs, I'd be really interested in finding out. The Do loop was just temporary while I tried to break down the problem to something simpler, in my current code, the values go into a table which then get's printed at the end of the loop. I may try out ParametricNDSolve, it seems like exactly what I want to do (though given how long it's taken me to get to grips with NDSolve, I'm a bit hesitant). 

Thank you both for your help. 
POSTED BY: Henry Banks
I think you need a semicolon after NDSolve inside the do loop.  (See CompoundExpression in the docs).  Also Do never produces output, so you won't see the output of NMaximize.

That said, if it were me, I'd try out Frank Kampas's advice.
POSTED BY: Michael Rogers
If you're going to solve the same problem repeated with NDSolve, only changing a parameter each time, you should use ParametricNDSolve, as that saves setting up the problem each time.  I realize that doesn't answer your question, however.  I suggest you post a simpler example that shows the same problem.
POSTED BY: Frank Kampas
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract