3

I have a large set of coupled non-linear first order ODEs describing the time evolution of multiple variables a[t]....z[t]. For most of these, I have fixed initial conditions. For three of these, however, say x[t], y[t] and z[t], I want the initial conditions to be parametres I can scan over. These are also the variables I am solving for.

Essentially my code looks something like this:

sol = ParametricNDSolve[{.....x'[t]==...,y'[t]==...,z'[t]==...,...x[0]==x0,
y[0]==y0, z[0]==z0}, {x, y, z}, {t, 0, 10}, {x0, y0, z0}]

Now, when I trial specific initial conditions x0, y0, z0 (using regular NDSolve) I am able to get the full solution of x, y and z evolution at any time, and able to plot them. This works fine.

To some extent I can also use ParametricNDSolve with the variable initial conditions and trial specific ones using xinstance=x[0.2,0.3,0.4] /. sol (say). However, I am having some trouble specifying a range of x0, y0, and z0 values and having ParametricNDSolve scan over that range (producing a family of solutions, which I could plot on one plot). I have tried using test[x0_,y0_,z0_]=NDSolveValue[] and scan across test values but can't seem to be able to get the plot I want. I have also tried putting the initial conditions in a table with n entries and writing sol=Table[NDSolveValue[],{i,1,n}] but am running into trouble with this method too.

Finally, I should mention that there are singularities in the functions at different t values (depending on initial conditions) where NDSolve halts and gives me errors, ideally I would avoid these by only solving/plotting x, y, z from 0 to 1.

Any help would be really appreciated.

SarahThompson
  • 425
  • 2
  • 8
  • You haven't presented any sample (system of ODE's), so it's quite hard to guess. – zhk Jun 06 '18 at 15:09
  • @zhk - I am not sure how to modify the code in the other answer for my problem - I have three variable initial conditions and eventually want a plot of x,y,z vs t. Also, the singularities are causing problems -- not sure how to restrict the range of x, y z to prevent this – SarahThompson Jun 06 '18 at 16:08
  • As for the sample of ODE's, I am not sure it will add much to the discussion - if you could provide a template of how to use NDSolve/ParametricNDSolve to find a family of solutions to the simplest possible 3 ODE's in x[t], y[t], z[t] with variable initial conditions specified by [x0min, x0max, step] (and similarly for others), for 0 < x[t], y[t], z[t] <1, that would already be very helpful – SarahThompson Jun 06 '18 at 16:14

1 Answers1

5

Since you didn't provide your equations, here is a made up example which I think does what you are attempting:

family = ParametricNDSolveValue[
  {x'[t] == x[t], y'[t] == 2 y[t], z'[t] == 3 z[t], x[0] == x0, y[0] == y0, z[0] == z0},
  {x[t], y[t], z[t]}, {t, 0, 1}, {x0, y0, z0}
 ];

initcond = 
  Table[{x0, y0, z0}, {x0, {10, 15}}, {y0, {30, 35}}, {z0, 0, 0.1, 0.05}]~Flatten~2;

interpolatingfunctions = family @@@ initcond;

Plot[interpolatingfunctions, {t, 0, 1}]

parametric function

plots

MarcoB
  • 67,153
  • 18
  • 91
  • 189
  • 1
    Thank you, this is exactly what I was looking for. I understand that you have 'applied' the initial conditions (defined separately as a table) to the function defined via ParametricNDSolve. The code makes sense and I am able to replicate the plot with your function.

    A quick follow up -- I get blank plots for my system most likely because x is singular at some t (depending on the initial conditions). How do I get it to plot until the singularity or, which would seem easier and actually more useful for me, plot the (monotonically increasing) x solution from 0 only until x=1?

     – SarahThompson Jun 06 '18 at 17:45
  • 2
    @SarahThompson I think your best bet would be to set a WhenEvent condition in your NDSolve. Essentially you tell NDSolve to monitor the value of a function or parameter, and to trigger an action (e.g. stop the integration). Something along the lines of WhenEvent[x[t] == 1, "StopIntegration". Take a look at the docs for more information. – MarcoB Jun 06 '18 at 17:57
  • 1
    How would I modify the code if I wanted to use a list of triplets list={{x0,y0,z0},{x1,y1,z1}...} where you have inserted trial x, y, z values? – SarahThompson Jun 10 '18 at 22:49
  • 2
    @SarahThompson If I understand what you mean, you would just leave out initcond and use interpolatingfunctions = family @@@ list. – MarcoB Jun 10 '18 at 23:13
  • 1
    Thanks, it works. I realised that I don't have to flatten it if it is already ordered as {{x0,y0,z0},{x1,y1,z1}...} – SarahThompson Jun 11 '18 at 10:09