Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I have a differential equation with parameters. The objective is to vary the parameters such that the first positive root of the solution is equal to a desired value obj.

As the equation is very involved, I decided that I'd use user1084363's question in Find all roots of an interpolating function (solution to a differential equation), which inspired the first part of my solution.

Currently, my solution is to manually manipulate the parameter k to obtain the desired solution. For example, if I'd like the first root of the equation to be 1.6, I would use the following line of code

{1.09 x''[t] - k* x'[t] + 1.1759 Sin[x[t]] == 0,  x[0] == Pi/3, x'[0] == 0}, x, {t, 0, 50}, 
Method -> {"EventLocator", "Event" -> x[t],"EventAction" :> Sow[t]}]]][[2]]-obj, 
{k, 0.01, 0.1}]

and manually vary the value of k till I obtain a value of zero.

However, this is quite time consuming and I was wondering whether there was a way to automate the process. My first thought was to use Root, with k as the variable, but this doesn't work as of now due to the presence of the differential equation.

share|improve this question
Would ParametricNDSolve help ? – b.gatessucks Jan 28 '13 at 9:04
I don't think it would work directly, but I'm trying to do a parametric sweep using ParametricNDSolve as suggested here – Vincent Tjeng Jan 28 '13 at 9:34
As a clarification, what I mean that it would not work directly is that the Reap function behaves differently with ParametricNDSolve. – Vincent Tjeng Jan 28 '13 at 9:35

I'd do :

auxSol = ParametricNDSolve[{1.09 x''[t] - k  x'[t] + 1.1759 Sin[x[t]] == 0, 
                             x[0] == Pi/3, x'[0] == 0}, x, {t, 0, 50}, {k}]

data = {#, FindRoot[(x[k][#] /. auxSol) == 0, {k, 0.01, 0.1}][[1, 2]]} & /@ Range[1, 2, 0.05]

ListPlot[data, Filling -> Axis, AxesLabel -> {"obj", "k"}]

enter image description here

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.