Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

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 reference.wolfram.com/mathematica/ref/ParametricNDSolve.html –  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
add comment

1 Answer

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
add comment

Your Answer

 
discard

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.