1
$\begingroup$

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

obj=1.6;
Manipulate[Flatten[Reap[NDSolve[
{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.

$\endgroup$
  • $\begingroup$ Would ParametricNDSolve help ? $\endgroup$ – b.gates.you.know.what Jan 28 '13 at 9:04
  • $\begingroup$ I don't think it would work directly, but I'm trying to do a parametric sweep using ParametricNDSolve as suggested here reference.wolfram.com/mathematica/ref/ParametricNDSolve.html $\endgroup$ – Vincent Tjeng Jan 28 '13 at 9:34
  • $\begingroup$ As a clarification, what I mean that it would not work directly is that the Reap function behaves differently with ParametricNDSolve. $\endgroup$ – Vincent Tjeng Jan 28 '13 at 9:35
2
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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