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.

Improve this question
$\endgroup$
3
  • $\begingroup$ Would ParametricNDSolve help ? $\endgroup$
    – b.gates.you.know.what
    Jan 28, 2013 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, 2013 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, 2013 at 9:35

1 Answer 1

Reset to default
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

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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